OGS 6.2.1-97-g73d1aeda3
PhaseFieldProcess.cpp
Go to the documentation of this file.
1 
12 #include "PhaseFieldProcess.h"
13 
14 #include <cassert>
15 
17 
18 #include "ProcessLib/Process.h"
20 
21 #include "PhaseFieldFEM.h"
22 
23 namespace ProcessLib
24 {
25 namespace PhaseField
26 {
27 template <int DisplacementDim>
29  std::string name,
30  MeshLib::Mesh& mesh,
31  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
32  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
33  unsigned const integration_order,
34  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
35  process_variables,
37  SecondaryVariableCollection&& secondary_variables,
38  NumLib::NamedFunctionCaller&& named_function_caller,
39  bool const use_monolithic_scheme)
40  : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
41  integration_order, std::move(process_variables),
42  std::move(secondary_variables), std::move(named_function_caller),
43  use_monolithic_scheme),
44  _process_data(std::move(process_data))
45 {
46  if (use_monolithic_scheme)
47  {
48  OGS_FATAL(
49  "Monolithic scheme is not implemented for the PhaseField process.");
50  }
51 
52  _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
53  mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
54 }
55 
56 template <int DisplacementDim>
58 {
59  return false;
60 }
61 
62 template <int DisplacementDim>
65  const int process_id) const
66 {
67  // For the M process (deformation) in the staggered scheme.
68  if (process_id == 0)
69  {
70  auto const& l = *_local_to_global_index_map;
71  return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
72  &l.getGhostIndices(), &this->_sparsity_pattern};
73  }
74 
75  // For staggered scheme and phase field process.
77  return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
78  &l.getGhostIndices(), &_sparsity_pattern_with_single_component};
79 }
80 
81 template <int DisplacementDim>
84 {
85  // For the M process (deformation) in the staggered scheme.
86  if (process_id == 0)
87  {
89  }
90 
91  // For the equation of phasefield
93 }
94 
95 template <int DisplacementDim>
97 {
98  // For displacement equation.
99  const int mechanics_process_id = 0;
101 
102  // TODO move the two data members somewhere else.
103  // for extrapolation of secondary variables of stress or strain
104  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
107  std::make_unique<NumLib::LocalToGlobalIndexMap>(
108  std::move(all_mesh_subsets_single_component),
109  // by location order is needed for output
111 
113 
114  // For phase field equation.
117 }
118 
119 template <int DisplacementDim>
121  NumLib::LocalToGlobalIndexMap const& dof_table,
122  MeshLib::Mesh const& mesh,
123  unsigned const integration_order)
124 {
126  DisplacementDim, PhaseFieldLocalAssembler>(
127  mesh.getElements(), dof_table, _local_assemblers,
128  mesh.isAxiallySymmetric(), integration_order, _process_data);
129 
131  "sigma",
133  DisplacementDim>::RowsAtCompileTime,
134  getExtrapolator(), _local_assemblers,
136 
138  "epsilon",
140  DisplacementDim>::RowsAtCompileTime,
141  getExtrapolator(), _local_assemblers,
143 
144  // Initialize local assemblers after all variables have been set.
146  &LocalAssemblerInterface::initialize, _local_assemblers,
148 }
149 
150 template <int DisplacementDim>
152 {
153  // Staggered scheme:
154  // for the equations of deformation.
155  const int mechanical_process_id = 0;
157  *_local_to_global_index_map, mechanical_process_id);
158  // for the phase field
159  const int phasefield_process_id = 1;
161  *_local_to_global_index_map_single_component, phasefield_process_id);
162 }
163 
164 template <int DisplacementDim>
166  const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
167  GlobalVector& b)
168 {
169  DBUG("Assemble PhaseFieldProcess.");
170 
171  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
172  dof_tables;
173 
174  // For the staggered scheme
175  if (_coupled_solutions->process_id == 1)
176  {
177  DBUG(
178  "Assemble the equations of phase field in "
179  "PhaseFieldProcess for the staggered scheme.");
180  }
181  else
182  {
183  DBUG(
184  "Assemble the equations of deformation in "
185  "PhaseFieldProcess for the staggered scheme.");
186  }
187  dof_tables.emplace_back(*_local_to_global_index_map_single_component);
188  dof_tables.emplace_back(*_local_to_global_index_map);
189 
190  const int process_id = _coupled_solutions->process_id;
191  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
192 
193  // Call global assembler for each local assembly item.
196  pv.getActiveElementIDs(), dof_tables, t, x, M, K, b,
198 }
199 
200 template <int DisplacementDim>
202  const double t, GlobalVector const& x, GlobalVector const& xdot,
203  const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
204  GlobalVector& b, GlobalMatrix& Jac)
205 {
206  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
207  dof_tables;
208 
209  // For the staggered scheme
210  if (_coupled_solutions->process_id == 1)
211  {
212  DBUG(
213  "Assemble the Jacobian equations of phase field in "
214  "PhaseFieldProcess for the staggered scheme.");
215  }
216  else
217  {
218  DBUG(
219  "Assemble the Jacobian equations of deformation in "
220  "PhaseFieldProcess for the staggered scheme.");
221  }
222  dof_tables.emplace_back(*_local_to_global_index_map);
223  dof_tables.emplace_back(*_local_to_global_index_map_single_component);
224 
225  const int process_id = _coupled_solutions->process_id;
226  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
227 
228  // Call global assembler for each local assembly item.
231  _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, x, xdot,
232  dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
233 
234  if (_coupled_solutions->process_id == 0)
235  {
236  b.copyValues(*_nodal_forces);
237  std::transform(_nodal_forces->begin(), _nodal_forces->end(),
238  _nodal_forces->begin(), [](double val) { return -val; });
239  }
240 }
241 
242 template <int DisplacementDim>
244  GlobalVector const& x, double const t, double const dt,
245  const int process_id)
246 {
247  DBUG("PreTimestep PhaseFieldProcess %d.", process_id);
248 
249  _process_data.dt = dt;
250  _process_data.t = t;
251  _process_data.injected_volume = _process_data.t;
252 
253  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
254 
257  pv.getActiveElementIDs(), getDOFTable(process_id), x, t, dt);
258 }
259 
260 template <int DisplacementDim>
262  GlobalVector const& x, const double /*t*/, const double /*delta_t*/,
263  int const process_id)
264 {
265  if (isPhaseFieldProcess(process_id))
266  {
267  DBUG("PostTimestep PhaseFieldProcess.");
268 
269  _process_data.elastic_energy = 0.0;
270  _process_data.surface_energy = 0.0;
271  _process_data.pressure_work = 0.0;
272 
273  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
274  dof_tables;
275 
276  dof_tables.emplace_back(*_local_to_global_index_map);
277  dof_tables.emplace_back(*_local_to_global_index_map_single_component);
278 
279  ProcessLib::ProcessVariable const& pv =
280  getProcessVariables(process_id)[0];
281 
284  pv.getActiveElementIDs(), dof_tables, x, _process_data.t,
285  _process_data.elastic_energy, _process_data.surface_energy,
286  _process_data.pressure_work, _coupled_solutions);
287 
288  INFO("Elastic energy: %g Surface energy: %g Pressure work: %g ",
289  _process_data.elastic_energy, _process_data.surface_energy,
290  _process_data.pressure_work);
291  }
292 }
293 
294 template <int DisplacementDim>
296  GlobalVector const& x, const double t, const int process_id)
297 {
298  _process_data.crack_volume = 0.0;
299 
300  if (!isPhaseFieldProcess(process_id))
301  {
302  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
303  dof_tables;
304 
305  dof_tables.emplace_back(*_local_to_global_index_map);
306  dof_tables.emplace_back(*_local_to_global_index_map_single_component);
307 
308  DBUG("PostNonLinearSolver crack volume computation.");
309 
310  ProcessLib::ProcessVariable const& pv =
311  getProcessVariables(process_id)[0];
314  pv.getActiveElementIDs(), dof_tables, x, t,
315  _process_data.crack_volume, _coupled_solutions);
316 
317  INFO("Integral of crack: %g", _process_data.crack_volume);
318 
319  if (_process_data.propagating_crack)
320  {
321  _process_data.pressure_old = _process_data.pressure;
322  _process_data.pressure =
323  _process_data.injected_volume / _process_data.crack_volume;
324  _process_data.pressure_error =
325  std::fabs(_process_data.pressure_old - _process_data.pressure) /
326  _process_data.pressure;
327  INFO("Internal pressure: %g and Pressure error: %.4e",
328  _process_data.pressure, _process_data.pressure_error);
329  auto& u = _coupled_solutions->coupled_xs[0].get();
330  MathLib::LinAlg::scale(const_cast<GlobalVector&>(u),
331  _process_data.pressure);
332  }
333  }
334  else
335  {
336  if (_process_data.propagating_crack)
337  {
338  auto& u = _coupled_solutions->coupled_xs[0].get();
339  MathLib::LinAlg::scale(const_cast<GlobalVector&>(u),
340  1 / _process_data.pressure);
341  }
342  }
343 }
344 
345 template <int DisplacementDim>
347  int const process_id) const
348 {
349  return process_id == 1;
350 }
351 
352 template class PhaseFieldProcess<2>;
353 template class PhaseFieldProcess<3>;
354 
355 } // namespace PhaseField
356 } // namespace ProcessLib
MeshLib::Mesh & _mesh
Definition: Process.h:287
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
virtual void computeCrackIntegral(std::size_t mesh_item_id, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, GlobalVector const &x, double const t, double &crack_volume, CoupledSolutionsForStaggeredScheme const *const cpl_xs)=0
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
Definition: Process.h:288
void initializeProcessBoundaryConditionsAndSourceTerms(const NumLib::LocalToGlobalIndexMap &dof_table, const int process_id)
Definition: Process.cpp:67
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:125
MeshLib::PropertyVector< double > * _nodal_forces
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map_single_component
static void executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)
A class to simulate mechanical fracturing process using phase-field approach in solids described by...
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, double const t, GlobalVector const &x, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, CoupledSolutionsForStaggeredScheme const *const cpl_xs)
void createLocalAssemblers(std::vector< MeshLib::Element *> const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface >> &local_assemblers, ExtraCtorArgs &&... extra_ctor_args)
void postTimestepConcreteProcess(GlobalVector const &x, const double t, const double delta_t, int const process_id) override
NumLib::Extrapolator & getExtrapolator() const
Definition: Process.h:152
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, const double t, GlobalVector const &x, GlobalVector const &xdot, const double dxdot_dx, const double dx_dx, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac, CoupledSolutionsForStaggeredScheme const *const cpl_xs)
bool isAxiallySymmetric() const
Definition: Mesh.h:137
virtual void initialize(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:290
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
std::vector< std::reference_wrapper< GlobalVector const > > const & coupled_xs
References to the current solutions of the coupled processes.
Builds expression trees of named functions dynamically at runtime.
virtual void preTimestep(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
bool isPhaseFieldProcess(int const process_id) const
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
SecondaryVariableCollection _secondary_variables
Definition: Process.h:292
void constructDofTableOfSpecifiedProsessStaggerdScheme(const int specified_prosess_id)
Definition: Process.cpp:281
void assembleWithJacobianConcreteProcess(const double t, GlobalVector const &x, GlobalVector const &xdot, const double dxdot_dx, const double dx_dx, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
std::vector< std::unique_ptr< LocalAssemblerInterface > > _local_assemblers
virtual std::vector< double > const & getIntPtSigma(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
virtual void computeEnergy(std::size_t mesh_item_id, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, GlobalVector const &x, double const t, double &elastic_energy, double &surface_energy, double &pressure_work, CoupledSolutionsForStaggeredScheme const *const cpl_xs)=0
GlobalSparsityPattern computeSparsityPattern(LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
Computes a sparsity pattern for the given inputs.
GlobalSparsityPattern _sparsity_pattern_with_single_component
virtual std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
PhaseFieldProcess(std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable >>> &&process_variables, PhaseFieldProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, NumLib::NamedFunctionCaller &&named_function_caller, bool const use_monolithic_scheme)
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:108
void postNonLinearSolverConcreteProcess(GlobalVector const &x, const double t, int const process_id) override
PhaseFieldProcessData< DisplacementDim > _process_data
CoupledSolutionsForStaggeredScheme * _coupled_solutions
Definition: Process.h:305
void preTimestepConcreteProcess(GlobalVector const &x, double const t, double const dt, const int process_id) override
Handles configuration of several secondary variables from the project file.
Ordering data by spatial location.
void assembleConcreteProcess(const double t, GlobalVector const &x, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
void addSecondaryVariable(std::string const &internal_name, SecondaryVariableFunctions &&fcts)
GlobalSparsityPattern _sparsity_pattern
Definition: Process.h:318
VectorMatrixAssembler _global_assembler
Definition: Process.h:299
void scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:43
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
std::vector< std::size_t > & getActiveElementIDs() const