OGS
RichardsMechanicsProcess.cpp
Go to the documentation of this file.
1 
12 
13 #include <cassert>
14 
15 #include "MeshLib/Elements/Utils.h"
20 #include "RichardsMechanicsFEM.h"
22 
23 namespace ProcessLib
24 {
25 namespace RichardsMechanics
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  bool const use_monolithic_scheme)
39  : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
40  integration_order, std::move(process_variables),
41  std::move(secondary_variables), use_monolithic_scheme),
42  _process_data(std::move(process_data))
43 {
44  _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
45  mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
46 
47  _hydraulic_flow = MeshLib::getOrCreateMeshProperty<double>(
48  mesh, "MassFlowRate", MeshLib::MeshItemType::Node, 1);
49 
50  // TODO (naumov) remove ip suffix. Probably needs modification of the mesh
51  // properties, s.t. there is no "overlapping" with cell/point data.
52  // See getOrCreateMeshProperty.
53  _integration_point_writer.emplace_back(
54  std::make_unique<IntegrationPointWriter>(
55  "sigma_ip",
56  static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
57  integration_order, _local_assemblers, &LocalAssemblerIF::getSigma));
58 
59  _integration_point_writer.emplace_back(
60  std::make_unique<IntegrationPointWriter>(
61  "saturation_ip", 1 /*n components*/, integration_order,
63 
64  _integration_point_writer.emplace_back(
65  std::make_unique<IntegrationPointWriter>(
66  "porosity_ip", 1 /*n components*/, integration_order,
68 
69  _integration_point_writer.emplace_back(
70  std::make_unique<IntegrationPointWriter>(
71  "transport_porosity_ip", 1 /*n components*/, integration_order,
73 
74  _integration_point_writer.emplace_back(
75  std::make_unique<IntegrationPointWriter>(
76  "swelling_stress_ip",
77  static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
78  integration_order, _local_assemblers,
80 
81  _integration_point_writer.emplace_back(
82  std::make_unique<IntegrationPointWriter>(
83  "epsilon_ip",
84  static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
85  integration_order, _local_assemblers,
87 }
88 
89 template <int DisplacementDim>
91 {
92  return false;
93 }
94 
95 template <int DisplacementDim>
98  const int process_id) const
99 {
100  // For the monolithic scheme or the M process (deformation) in the staggered
101  // scheme.
102  if (_use_monolithic_scheme || process_id == 1)
103  {
104  auto const& l = *_local_to_global_index_map;
105  return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
106  &l.getGhostIndices(), &this->_sparsity_pattern};
107  }
108 
109  // For staggered scheme and H process (pressure).
110  auto const& l = *_local_to_global_index_map_with_base_nodes;
111  return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
112  &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
113 }
114 
115 template <int DisplacementDim>
117 {
118  // Create single component dof in every of the mesh's nodes.
119  _mesh_subset_all_nodes =
120  std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
121  // Create single component dof in the mesh's base nodes.
122  _base_nodes = MeshLib::getBaseNodes(_mesh.getElements());
123  _mesh_subset_base_nodes =
124  std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
125 
126  // TODO move the two data members somewhere else.
127  // for extrapolation of secondary variables of stress or strain
128  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
129  *_mesh_subset_all_nodes};
130  _local_to_global_index_map_single_component =
131  std::make_unique<NumLib::LocalToGlobalIndexMap>(
132  std::move(all_mesh_subsets_single_component),
133  // by location order is needed for output
135 
136  if (_use_monolithic_scheme)
137  {
138  // For pressure, which is the first
139  std::vector<MeshLib::MeshSubset> all_mesh_subsets{
140  *_mesh_subset_base_nodes};
141 
142  // For displacement.
143  const int monolithic_process_id = 0;
144  std::generate_n(std::back_inserter(all_mesh_subsets),
145  getProcessVariables(monolithic_process_id)[1]
146  .get()
147  .getNumberOfGlobalComponents(),
148  [&]() { return *_mesh_subset_all_nodes; });
149 
150  std::vector<int> const vec_n_components{1, DisplacementDim};
151  _local_to_global_index_map =
152  std::make_unique<NumLib::LocalToGlobalIndexMap>(
153  std::move(all_mesh_subsets), vec_n_components,
155  assert(_local_to_global_index_map);
156  }
157  else
158  {
159  // For displacement equation.
160  const int process_id = 1;
161  std::vector<MeshLib::MeshSubset> all_mesh_subsets;
162  std::generate_n(std::back_inserter(all_mesh_subsets),
163  getProcessVariables(process_id)[0]
164  .get()
165  .getNumberOfGlobalComponents(),
166  [&]() { return *_mesh_subset_all_nodes; });
167 
168  std::vector<int> const vec_n_components{DisplacementDim};
169  _local_to_global_index_map =
170  std::make_unique<NumLib::LocalToGlobalIndexMap>(
171  std::move(all_mesh_subsets), vec_n_components,
173 
174  // For pressure equation.
175  // Collect the mesh subsets with base nodes in a vector.
176  std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
177  *_mesh_subset_base_nodes};
178  _local_to_global_index_map_with_base_nodes =
179  std::make_unique<NumLib::LocalToGlobalIndexMap>(
180  std::move(all_mesh_subsets_base_nodes),
181  // by location order is needed for output
183 
184  _sparsity_pattern_with_linear_element = NumLib::computeSparsityPattern(
185  *_local_to_global_index_map_with_base_nodes, _mesh);
186 
187  assert(_local_to_global_index_map);
188  assert(_local_to_global_index_map_with_base_nodes);
189  }
190 }
191 
192 template <int DisplacementDim>
194  NumLib::LocalToGlobalIndexMap const& dof_table,
195  MeshLib::Mesh const& mesh,
196  unsigned const integration_order)
197 {
198  using nlohmann::json;
199 
200  ProcessLib::createLocalAssemblersHM<DisplacementDim,
202  mesh.getElements(), dof_table, _local_assemblers,
203  NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
204  _process_data);
205 
206  auto add_secondary_variable = [&](std::string const& name,
207  int const num_components,
208  auto get_ip_values_function)
209  {
210  _secondary_variables.addSecondaryVariable(
211  name,
212  makeExtrapolator(num_components, getExtrapolator(),
213  _local_assemblers,
214  std::move(get_ip_values_function)));
215  };
216 
217  add_secondary_variable("sigma",
219  DisplacementDim>::RowsAtCompileTime,
220  &LocalAssemblerIF::getIntPtSigma);
221 
222  add_secondary_variable("swelling_stress",
224  DisplacementDim>::RowsAtCompileTime,
225  &LocalAssemblerIF::getIntPtSwellingStress);
226 
227  add_secondary_variable("epsilon",
229  DisplacementDim>::RowsAtCompileTime,
230  &LocalAssemblerIF::getIntPtEpsilon);
231 
232  add_secondary_variable("velocity", DisplacementDim,
233  &LocalAssemblerIF::getIntPtDarcyVelocity);
234 
235  add_secondary_variable("saturation", 1,
236  &LocalAssemblerIF::getIntPtSaturation);
237 
238  add_secondary_variable("micro_saturation", 1,
239  &LocalAssemblerIF::getIntPtMicroSaturation);
240 
241  add_secondary_variable("micro_pressure", 1,
242  &LocalAssemblerIF::getIntPtMicroPressure);
243 
244  add_secondary_variable("porosity", 1, &LocalAssemblerIF::getIntPtPorosity);
245 
246  add_secondary_variable("transport_porosity", 1,
247  &LocalAssemblerIF::getIntPtTransportPorosity);
248 
249  add_secondary_variable("dry_density_solid", 1,
250  &LocalAssemblerIF::getIntPtDryDensitySolid);
251 
252  //
253  // enable output of internal variables defined by material models
254  //
256  LocalAssemblerIF>(_process_data.solid_materials,
257  add_secondary_variable);
258 
259  // Assume all materials have same internal variables.
262  _process_data.solid_materials, _local_assemblers,
263  _integration_point_writer, integration_order);
264 
265  _process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
266  const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
268 
269  _process_data.element_porosity = MeshLib::getOrCreateMeshProperty<double>(
270  const_cast<MeshLib::Mesh&>(mesh), "porosity_avg",
272 
273  _process_data.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
274  const_cast<MeshLib::Mesh&>(mesh), "stress_avg",
277  DisplacementDim>::RowsAtCompileTime);
278 
279  _process_data.pressure_interpolated =
280  MeshLib::getOrCreateMeshProperty<double>(
281  const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
283 
284  setIPDataInitialConditions(_integration_point_writer, mesh.getProperties(),
285  _local_assemblers);
286 
287  // Initialize local assemblers after all variables have been set.
288  GlobalExecutor::executeMemberOnDereferenced(&LocalAssemblerIF::initialize,
289  _local_assemblers,
290  *_local_to_global_index_map);
291 }
292 
293 template <int DisplacementDim>
295 {
296  if (_use_monolithic_scheme)
297  {
298  const int monolithic_process_id = 0;
299  initializeProcessBoundaryConditionsAndSourceTerms(
300  *_local_to_global_index_map, monolithic_process_id);
301  return;
302  }
303 
304  // Staggered scheme:
305  // for the equations of pressure
306  const int hydraulic_process_id = 0;
307  initializeProcessBoundaryConditionsAndSourceTerms(
308  *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
309 
310  // for the equations of deformation.
311  const int mechanical_process_id = 1;
312  initializeProcessBoundaryConditionsAndSourceTerms(
313  *_local_to_global_index_map, mechanical_process_id);
314 }
315 
316 template <int DisplacementDim>
318  setInitialConditionsConcreteProcess(std::vector<GlobalVector*>& x,
319  double const t,
320  int const process_id)
321 {
322  if (process_id != 0)
323  {
324  return;
325  }
326 
327  DBUG("SetInitialConditions RichardsMechanicsProcess.");
328 
329  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
330 
332  &LocalAssemblerIF::setInitialConditions, _local_assemblers,
333  pv.getActiveElementIDs(), getDOFTable(process_id), *x[process_id], t,
334  _use_monolithic_scheme, process_id);
335 }
336 
337 template <int DisplacementDim>
339  const double t, double const dt, std::vector<GlobalVector*> const& x,
340  std::vector<GlobalVector*> const& xdot, int const process_id,
342 {
343  DBUG("Assemble the equations for RichardsMechanics");
344 
345  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
346  dof_table = {std::ref(*_local_to_global_index_map)};
347  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
348 
349  // Call global assembler for each local assembly item.
351  _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
352  pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
353  b);
354 }
355 
356 template <int DisplacementDim>
358  assembleWithJacobianConcreteProcess(const double t, double const dt,
359  std::vector<GlobalVector*> const& x,
360  std::vector<GlobalVector*> const& xdot,
361  int const process_id, GlobalMatrix& M,
362  GlobalMatrix& K, GlobalVector& b,
363  GlobalMatrix& Jac)
364 {
365  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
366  dof_tables;
367  // For the monolithic scheme
368  if (_use_monolithic_scheme)
369  {
370  DBUG(
371  "Assemble the Jacobian of RichardsMechanics for the monolithic"
372  " scheme.");
373  dof_tables.emplace_back(*_local_to_global_index_map);
374  }
375  else
376  {
377  // For the staggered scheme
378  if (process_id == 0)
379  {
380  DBUG(
381  "Assemble the Jacobian equations of liquid fluid process in "
382  "RichardsMechanics for the staggered scheme.");
383  }
384  else
385  {
386  DBUG(
387  "Assemble the Jacobian equations of mechanical process in "
388  "RichardsMechanics for the staggered scheme.");
389  }
390  dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
391  dof_tables.emplace_back(*_local_to_global_index_map);
392  }
393 
394  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
395 
398  _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
399  process_id, M, K, b, Jac);
400 
401  auto copyRhs = [&](int const variable_id, auto& output_vector)
402  {
403  if (_use_monolithic_scheme)
404  {
405  transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
406  output_vector,
407  std::negate<double>());
408  }
409  else
410  {
411  transformVariableFromGlobalVector(b, 0, dof_tables[process_id],
412  output_vector,
413  std::negate<double>());
414  }
415  };
416  if (_use_monolithic_scheme || process_id == 0)
417  {
418  copyRhs(0, *_hydraulic_flow);
419  }
420  if (_use_monolithic_scheme || process_id == 1)
421  {
422  copyRhs(1, *_nodal_forces);
423  }
424 }
425 
426 template <int DisplacementDim>
428  std::vector<GlobalVector*> const& x, double const t, double const dt,
429  const int process_id)
430 {
431  if (hasMechanicalProcess(process_id))
432  {
433  DBUG("PostTimestep RichardsMechanicsProcess.");
434  auto const dof_tables = getDOFTables(x.size());
435 
436  ProcessLib::ProcessVariable const& pv =
437  getProcessVariables(process_id)[0];
439  &LocalAssemblerIF::postTimestep, _local_assemblers,
440  pv.getActiveElementIDs(), dof_tables, x, t, dt);
441  }
442 }
443 
444 template <int DisplacementDim>
446  computeSecondaryVariableConcrete(const double t, const double dt,
447  std::vector<GlobalVector*> const& x,
448  GlobalVector const& x_dot,
449  int const process_id)
450 {
451  if (process_id != 0)
452  {
453  return;
454  }
455 
456  DBUG("Compute the secondary variables for RichardsMechanicsProcess.");
457  auto const dof_tables = getDOFTables(x.size());
458 
459  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
461  &LocalAssemblerIF::computeSecondaryVariable, _local_assemblers,
462  pv.getActiveElementIDs(), dof_tables, t, dt, x, x_dot, process_id);
463 }
464 
465 template <int DisplacementDim>
466 std::tuple<NumLib::LocalToGlobalIndexMap*, bool> RichardsMechanicsProcess<
467  DisplacementDim>::getDOFTableForExtrapolatorData() const
468 {
469  const bool manage_storage = false;
470  return std::make_tuple(_local_to_global_index_map_single_component.get(),
471  manage_storage);
472 }
473 
474 template <int DisplacementDim>
477  const int process_id) const
478 {
479  if (hasMechanicalProcess(process_id))
480  {
481  return *_local_to_global_index_map;
482  }
483 
484  // For the equation of pressure
485  return *_local_to_global_index_map_with_base_nodes;
486 }
487 
488 template <int DisplacementDim>
489 std::vector<NumLib::LocalToGlobalIndexMap const*>
491  int const number_of_processes) const
492 {
493  std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
494  dof_tables.reserve(number_of_processes);
495  std::generate_n(std::back_inserter(dof_tables), number_of_processes,
496  [&]() { return &getDOFTable(dof_tables.size()); });
497  return dof_tables;
498 }
499 
500 template class RichardsMechanicsProcess<2>;
501 template class RichardsMechanicsProcess<3>;
502 
503 } // namespace RichardsMechanics
504 } // namespace ProcessLib
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:29
Global vector based on Eigen vector.
Definition: EigenVector.h:28
bool isAxiallySymmetric() const
Definition: Mesh.h:131
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:76
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:103
Properties & getProperties()
Definition: Mesh.h:128
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::unique_ptr< IntegrationPointWriter > > _integration_point_writer
Definition: Process.h:354
std::vector< std::unique_ptr< LocalAssemblerIF > > _local_assemblers
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, int const process_id) override
std::vector< NumLib::LocalToGlobalIndexMap const * > getDOFTables(const int number_of_processes) const
RichardsMechanicsProcess(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, RichardsMechanicsProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &x, double const t, int const process_id) override
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
Handles configuration of several secondary variables from the project file.
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, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac)
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, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:48
static const double t
std::vector< Node * > getBaseNodes(std::vector< Element * > const &elements)
Definition: Utils.h:26
@ BY_LOCATION
Ordering data by spatial location.
GlobalSparsityPattern computeSparsityPattern(LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
Computes a sparsity pattern for the given inputs.
void transformVariableFromGlobalVector(GlobalVector const &input_vector, int const variable_id, NumLib::LocalToGlobalIndexMap const &local_to_global_index_map, MeshLib::PropertyVector< double > &output_vector, Functor mapFunction)
Definition: DOFTableUtil.h:59
void solidMaterialInternalToSecondaryVariables(std::map< int, std::unique_ptr< MaterialLib::Solids::MechanicsBase< DisplacementDim >>> const &solid_materials, AddSecondaryVariableCallback const &add_secondary_variable)
void solidMaterialInternalVariablesToIntegrationPointWriter(std::map< int, std::unique_ptr< MaterialLib::Solids::MechanicsBase< DisplacementDim >>> const &solid_materials, std::vector< std::unique_ptr< LocalAssemblerInterface >> const &local_assemblers, std::vector< std::unique_ptr< IntegrationPointWriter >> &integration_point_writer, int const integration_order)
void setIPDataInitialConditions(std::vector< std::unique_ptr< IntegrationPointWriter >> const &_integration_point_writer, MeshLib::Properties const &mesh_properties, LocalAssemblersVector &local_assemblers)
std::pair< std::vector< GlobalVector * >, std::vector< GlobalVector * > > setInitialConditions(double const t0, std::vector< std::unique_ptr< ProcessData >> const &per_process_data)
Definition: TimeLoop.cpp:171
void createLocalAssemblersHM(std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface >> &local_assemblers, NumLib::IntegrationOrder const integration_order, ExtraCtorArgs &&... extra_ctor_args)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)
virtual std::vector< double > getSwellingStress() const =0
virtual std::vector< double > getSaturation() const =0
virtual std::vector< double > getSigma() const =0
virtual std::vector< double > getPorosity() const =0
virtual std::vector< double > getTransportPorosity() const =0
virtual std::vector< double > getEpsilon() const =0