OGS 6.2.0-97-g4a610c866
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  MeshLib::Mesh& mesh,
30  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
31  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
32  unsigned const integration_order,
33  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
34  process_variables,
36  SecondaryVariableCollection&& secondary_variables,
37  NumLib::NamedFunctionCaller&& named_function_caller,
38  bool const use_monolithic_scheme)
39  : Process(mesh, std::move(jacobian_assembler), parameters,
40  integration_order, std::move(process_variables),
41  std::move(secondary_variables), std::move(named_function_caller),
42  use_monolithic_scheme),
43  _process_data(std::move(process_data))
44 {
45  if (use_monolithic_scheme)
46  {
47  OGS_FATAL(
48  "Monolithic scheme is not implemented for the PhaseField process.");
49  }
50 
51  _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
52  mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
53 }
54 
55 template <int DisplacementDim>
57 {
58  return false;
59 }
60 
61 template <int DisplacementDim>
64  const int process_id) const
65 {
66  // For the M process (deformation) in the staggered scheme.
67  if (process_id == 0)
68  {
69  auto const& l = *_local_to_global_index_map;
70  return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
71  &l.getGhostIndices(), &this->_sparsity_pattern};
72  }
73 
74  // For staggered scheme and phase field process.
76  return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
77  &l.getGhostIndices(), &_sparsity_pattern_with_single_component};
78 }
79 
80 template <int DisplacementDim>
83 {
84  // For the M process (deformation) in the staggered scheme.
85  if (process_id == 0)
86  {
88  }
89 
90  // For the equation of phasefield
92 }
93 
94 template <int DisplacementDim>
96 {
97  // Create single component dof in every of the mesh's nodes.
99  std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
100 
101  // TODO move the two data members somewhere else.
102  // for extrapolation of secondary variables of stress or strain
103  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
106  std::make_unique<NumLib::LocalToGlobalIndexMap>(
107  std::move(all_mesh_subsets_single_component),
108  // by location order is needed for output
110 
112 
113  // For displacement equation.
114  const int process_id = 0;
115  std::vector<MeshLib::MeshSubset> all_mesh_subsets;
116  std::generate_n(
117  std::back_inserter(all_mesh_subsets),
118  getProcessVariables(process_id)[0].get().getNumberOfComponents(),
119  [&]() { return *_mesh_subset_all_nodes; });
120 
121  std::vector<int> const vec_n_components{DisplacementDim};
123  std::make_unique<NumLib::LocalToGlobalIndexMap>(
124  std::move(all_mesh_subsets), vec_n_components,
126 
127  // For phase field equation.
130 }
131 
132 template <int DisplacementDim>
134  NumLib::LocalToGlobalIndexMap const& dof_table,
135  MeshLib::Mesh const& mesh,
136  unsigned const integration_order)
137 {
139  DisplacementDim, PhaseFieldLocalAssembler>(
140  mesh.getElements(), dof_table, _local_assemblers,
141  mesh.isAxiallySymmetric(), integration_order, _process_data);
142 
144  "sigma",
146  DisplacementDim>::RowsAtCompileTime,
147  getExtrapolator(), _local_assemblers,
149 
151  "epsilon",
153  DisplacementDim>::RowsAtCompileTime,
154  getExtrapolator(), _local_assemblers,
156 }
157 
158 template <int DisplacementDim>
160 {
161  // Staggered scheme:
162  // for the equations of deformation.
163  const int mechanical_process_id = 0;
165  *_local_to_global_index_map, mechanical_process_id);
166  // for the phase field
167  const int phasefield_process_id = 1;
169  *_local_to_global_index_map_single_component, phasefield_process_id);
170 }
171 
172 template <int DisplacementDim>
174  const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
175  GlobalVector& b)
176 {
177  DBUG("Assemble PhaseFieldProcess.");
178 
179  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
180  dof_tables;
181 
182  // For the staggered scheme
183  if (_coupled_solutions->process_id == 1)
184  {
185  DBUG(
186  "Assemble the equations of phase field in "
187  "PhaseFieldProcess for the staggered scheme.");
188  }
189  else
190  {
191  DBUG(
192  "Assemble the equations of deformation in "
193  "PhaseFieldProcess for the staggered scheme.");
194  }
195  dof_tables.emplace_back(*_local_to_global_index_map_single_component);
196  dof_tables.emplace_back(*_local_to_global_index_map);
197 
198  const int process_id = _coupled_solutions->process_id;
199  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
200 
201  // Call global assembler for each local assembly item.
204  pv.getActiveElementIDs(), dof_tables, t, x, M, K, b,
206 }
207 
208 template <int DisplacementDim>
210  const double t, GlobalVector const& x, GlobalVector const& xdot,
211  const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
212  GlobalVector& b, GlobalMatrix& Jac)
213 {
214  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
215  dof_tables;
216 
217  // For the staggered scheme
218  if (_coupled_solutions->process_id == 1)
219  {
220  DBUG(
221  "Assemble the Jacobian equations of phase field in "
222  "PhaseFieldProcess for the staggered scheme.");
223  }
224  else
225  {
226  DBUG(
227  "Assemble the Jacobian equations of deformation in "
228  "PhaseFieldProcess for the staggered scheme.");
229  }
230  dof_tables.emplace_back(*_local_to_global_index_map);
231  dof_tables.emplace_back(*_local_to_global_index_map_single_component);
232 
233  const int process_id = _coupled_solutions->process_id;
234  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
235 
236  // Call global assembler for each local assembly item.
239  _local_assemblers, pv.getActiveElementIDs(), dof_tables, t,
240  x, xdot, dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
241 
242  if (_coupled_solutions->process_id == 0)
243  {
244  b.copyValues(*_nodal_forces);
245  std::transform(_nodal_forces->begin(), _nodal_forces->end(),
246  _nodal_forces->begin(), [](double val) { return -val; });
247  }
248 }
249 
250 template <int DisplacementDim>
252  GlobalVector const& x, double const t, double const dt,
253  const int process_id)
254 {
255  DBUG("PreTimestep PhaseFieldProcess %d.", process_id);
256 
257  _process_data.dt = dt;
258  _process_data.t = t;
259  _process_data.injected_volume = _process_data.t;
260 
261  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
262 
265  pv.getActiveElementIDs(), getDOFTable(process_id), x, t, dt);
266 }
267 
268 template <int DisplacementDim>
270  GlobalVector const& x, const double /*t*/, const double /*delta_t*/,
271  int const process_id)
272 {
273  if (isPhaseFieldProcess(process_id))
274  {
275  DBUG("PostTimestep PhaseFieldProcess.");
276 
277  _process_data.elastic_energy = 0.0;
278  _process_data.surface_energy = 0.0;
279  _process_data.pressure_work = 0.0;
280 
281  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
282  dof_tables;
283 
284  dof_tables.emplace_back(*_local_to_global_index_map);
285  dof_tables.emplace_back(*_local_to_global_index_map_single_component);
286 
288  = getProcessVariables(process_id)[0];
289 
292  pv.getActiveElementIDs(), dof_tables, x, _process_data.t,
293  _process_data.elastic_energy, _process_data.surface_energy,
294  _process_data.pressure_work, _coupled_solutions);
295 
296  INFO("Elastic energy: %g Surface energy: %g Pressure work: %g ",
297  _process_data.elastic_energy, _process_data.surface_energy,
298  _process_data.pressure_work);
299  }
300 }
301 
302 template <int DisplacementDim>
304  GlobalVector const& x, const double t, const int process_id)
305 {
306  _process_data.crack_volume = 0.0;
307 
308  if (!isPhaseFieldProcess(process_id))
309  {
310  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
311  dof_tables;
312 
313  dof_tables.emplace_back(*_local_to_global_index_map);
314  dof_tables.emplace_back(*_local_to_global_index_map_single_component);
315 
316  DBUG("PostNonLinearSolver crack volume computation.");
317 
319  = getProcessVariables(process_id)[0];
322  pv.getActiveElementIDs(), dof_tables, x, t,
323  _process_data.crack_volume, _coupled_solutions);
324 
325  INFO("Integral of crack: %g", _process_data.crack_volume);
326 
327  if (_process_data.propagating_crack)
328  {
329  _process_data.pressure_old = _process_data.pressure;
330  _process_data.pressure =
331  _process_data.injected_volume / _process_data.crack_volume;
332  _process_data.pressure_error =
333  std::fabs(_process_data.pressure_old - _process_data.pressure) /
334  _process_data.pressure;
335  INFO("Internal pressure: %g and Pressure error: %.4e",
336  _process_data.pressure, _process_data.pressure_error);
337  auto& u = _coupled_solutions->coupled_xs[0].get();
338  MathLib::LinAlg::scale(const_cast<GlobalVector&>(u),
339  _process_data.pressure);
340  }
341  }
342  else
343  {
344  if (_process_data.propagating_crack)
345  {
346  auto& u = _coupled_solutions->coupled_xs[0].get();
347  MathLib::LinAlg::scale(const_cast<GlobalVector&>(u),
348  1 / _process_data.pressure);
349  }
350  }
351 }
352 
353 template <int DisplacementDim>
355  int const process_id) const
356 {
357  return process_id == 1;
358 }
359 
360 template class PhaseFieldProcess<2>;
361 template class PhaseFieldProcess<3>;
362 
363 } // namespace PhaseField
364 } // namespace ProcessLib
MeshLib::Mesh & _mesh
Definition: Process.h:262
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
PhaseFieldProcess(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::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
Definition: Process.h:263
void initializeProcessBoundaryConditionsAndSourceTerms(const NumLib::LocalToGlobalIndexMap &dof_table, const int process_id)
Definition: Process.cpp:65
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:124
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
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:105
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:151
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
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:265
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:267
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
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:280
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:293
VectorMatrixAssembler _global_assembler
Definition: Process.h:274
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