OGS
HydroMechanicsProcess.cpp
Go to the documentation of this file.
1 
11 #include "HydroMechanicsProcess.h"
12 
13 #include "LocalAssembler/CreateLocalAssemblers.h"
18 #include "MeshLib/ElementStatus.h"
19 #include "MeshLib/Elements/Utils.h"
20 #include "MeshLib/Mesh.h"
22 #include "MeshLib/Properties.h"
29 
30 namespace ProcessLib
31 {
32 namespace LIE
33 {
34 namespace HydroMechanics
35 {
36 template <int GlobalDim>
38  std::string name,
39  MeshLib::Mesh& mesh,
40  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
41  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
42  unsigned const integration_order,
43  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
44  process_variables,
46  SecondaryVariableCollection&& secondary_variables,
47  bool const use_monolithic_scheme)
48  : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
49  integration_order, std::move(process_variables),
50  std::move(secondary_variables), use_monolithic_scheme),
51  _process_data(std::move(process_data))
52 {
53  INFO("[LIE/HM] looking for fracture elements in the given mesh");
54  std::vector<int> vec_fracture_mat_IDs;
55  std::vector<std::vector<MeshLib::Element*>> vec_vec_fracture_elements;
56  std::vector<std::vector<MeshLib::Element*>>
57  vec_vec_fracture_matrix_elements;
58  std::vector<std::vector<MeshLib::Node*>> vec_vec_fracture_nodes;
59  std::vector<std::pair<std::size_t, std::vector<int>>>
60  vec_branch_nodeID_matIDs;
61  std::vector<std::pair<std::size_t, std::vector<int>>>
62  vec_junction_nodeID_matIDs;
64  mesh, _vec_matrix_elements, vec_fracture_mat_IDs,
65  vec_vec_fracture_elements, vec_vec_fracture_matrix_elements,
66  vec_vec_fracture_nodes, vec_branch_nodeID_matIDs,
67  vec_junction_nodeID_matIDs);
69  vec_vec_fracture_elements[0].begin(),
70  vec_vec_fracture_elements[0].end());
73  vec_vec_fracture_matrix_elements[0].begin(),
74  vec_vec_fracture_matrix_elements[0].end());
76  vec_vec_fracture_nodes[0].begin(),
77  vec_vec_fracture_nodes[0].end());
78 
79  if (!_vec_fracture_elements.empty())
80  {
81  // set fracture property assuming a fracture forms a straight line
82  setFractureProperty(GlobalDim,
84  *_process_data.fracture_property);
85  }
86 
87  //
88  // If Neumann BCs for the displacement_jump variable are required they need
89  // special treatment because of the levelset function. The implementation
90  // exists in the version 6.1.0 (e54815cc07ee89c81f953a4955b1c788595dd725)
91  // and was removed due to lack of applications.
92  //
93 
94  if (!_process_data.deactivate_matrix_in_flow)
95  {
96  _process_data.p_element_status =
97  std::make_unique<MeshLib::ElementStatus>(&mesh);
98  }
99  else
100  {
101  auto const range =
103  if (!range)
104  {
105  OGS_FATAL(
106  "Could not get minimum/maximum ranges values for the "
107  "MaterialIDs property in the mesh '{:s}'.",
108  mesh.getName());
109  }
110 
111  std::vector<int> vec_p_inactive_matIDs;
112  for (int matID = range->first; matID <= range->second; matID++)
113  {
114  if (std::find(vec_fracture_mat_IDs.begin(),
115  vec_fracture_mat_IDs.end(),
116  matID) == vec_fracture_mat_IDs.end())
117  {
118  vec_p_inactive_matIDs.push_back(matID);
119  }
120  }
121  _process_data.p_element_status =
122  std::make_unique<MeshLib::ElementStatus>(&mesh,
123  vec_p_inactive_matIDs);
124 
125  const int monolithic_process_id = 0;
126  ProcessVariable const& pv_p =
127  getProcessVariables(monolithic_process_id)[0];
128  _process_data.p0 = &pv_p.getInitialCondition();
129  }
130 }
131 
132 template <int GlobalDim>
134 {
135  //------------------------------------------------------------
136  // prepare mesh subsets to define DoFs
137  //------------------------------------------------------------
138  // for extrapolation
139  _mesh_subset_all_nodes =
140  std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
141  // pressure
142  _mesh_nodes_p = MeshLib::getBaseNodes(
143  _process_data.p_element_status->getActiveElements());
144  _mesh_subset_nodes_p =
145  std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh_nodes_p);
146  // regular u
147  _mesh_subset_matrix_nodes =
148  std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
149  if (!_vec_fracture_nodes.empty())
150  {
151  // u jump
152  _mesh_subset_fracture_nodes =
153  std::make_unique<MeshLib::MeshSubset>(_mesh, _vec_fracture_nodes);
154  }
155 
156  // Collect the mesh subsets in a vector.
157  std::vector<MeshLib::MeshSubset> all_mesh_subsets;
158  std::vector<int> vec_n_components;
159  std::vector<std::vector<MeshLib::Element*> const*> vec_var_elements;
160  // pressure
161  vec_n_components.push_back(1);
162  all_mesh_subsets.emplace_back(*_mesh_subset_nodes_p);
163  if (!_process_data.deactivate_matrix_in_flow)
164  {
165  vec_var_elements.push_back(&_mesh.getElements());
166  }
167  else
168  {
169  // TODO set elements including active nodes for pressure.
170  // cannot use ElementStatus
171  vec_var_elements.push_back(&_vec_fracture_matrix_elements);
172  }
173  // regular displacement
174  vec_n_components.push_back(GlobalDim);
175  std::generate_n(std::back_inserter(all_mesh_subsets), GlobalDim,
176  [&]() { return *_mesh_subset_matrix_nodes; });
177  vec_var_elements.push_back(&_vec_matrix_elements);
178  if (!_vec_fracture_nodes.empty())
179  {
180  // displacement jump
181  vec_n_components.push_back(GlobalDim);
182  std::generate_n(std::back_inserter(all_mesh_subsets), GlobalDim,
183  [&]() { return *_mesh_subset_fracture_nodes; });
184  vec_var_elements.push_back(&_vec_fracture_matrix_elements);
185  }
186 
187  INFO("[LIE/HM] creating a DoF table");
188  _local_to_global_index_map =
189  std::make_unique<NumLib::LocalToGlobalIndexMap>(
190  std::move(all_mesh_subsets),
191  vec_n_components,
192  vec_var_elements,
194 
195  DBUG("created {:d} DoF", _local_to_global_index_map->size());
196 }
197 
198 template <int GlobalDim>
200  NumLib::LocalToGlobalIndexMap const& dof_table,
201  MeshLib::Mesh const& mesh,
202  unsigned const integration_order)
203 {
204  assert(mesh.getDimension() == GlobalDim);
205  INFO("[LIE/HM] creating local assemblers");
210  mesh.getElements(), dof_table, _local_assemblers,
211  NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
212  _process_data);
213 
214  auto mesh_prop_sigma_xx = MeshLib::getOrCreateMeshProperty<double>(
215  const_cast<MeshLib::Mesh&>(mesh), "stress_xx",
217  mesh_prop_sigma_xx->resize(mesh.getNumberOfElements());
218  _process_data.mesh_prop_stress_xx = mesh_prop_sigma_xx;
219 
220  auto mesh_prop_sigma_yy = MeshLib::getOrCreateMeshProperty<double>(
221  const_cast<MeshLib::Mesh&>(mesh), "stress_yy",
223  mesh_prop_sigma_yy->resize(mesh.getNumberOfElements());
224  _process_data.mesh_prop_stress_yy = mesh_prop_sigma_yy;
225 
226  auto mesh_prop_sigma_zz = MeshLib::getOrCreateMeshProperty<double>(
227  const_cast<MeshLib::Mesh&>(mesh), "stress_zz",
229  mesh_prop_sigma_zz->resize(mesh.getNumberOfElements());
230  _process_data.mesh_prop_stress_zz = mesh_prop_sigma_zz;
231 
232  auto mesh_prop_sigma_xy = MeshLib::getOrCreateMeshProperty<double>(
233  const_cast<MeshLib::Mesh&>(mesh), "stress_xy",
235  mesh_prop_sigma_xy->resize(mesh.getNumberOfElements());
236  _process_data.mesh_prop_stress_xy = mesh_prop_sigma_xy;
237 
238  if (GlobalDim == 3)
239  {
240  auto mesh_prop_sigma_xz = MeshLib::getOrCreateMeshProperty<double>(
241  const_cast<MeshLib::Mesh&>(mesh), "stress_xz",
243  mesh_prop_sigma_xz->resize(mesh.getNumberOfElements());
244  _process_data.mesh_prop_stress_xz = mesh_prop_sigma_xz;
245 
246  auto mesh_prop_sigma_yz = MeshLib::getOrCreateMeshProperty<double>(
247  const_cast<MeshLib::Mesh&>(mesh), "stress_yz",
249  mesh_prop_sigma_yz->resize(mesh.getNumberOfElements());
250  _process_data.mesh_prop_stress_yz = mesh_prop_sigma_yz;
251  }
252 
253  auto mesh_prop_epsilon_xx = MeshLib::getOrCreateMeshProperty<double>(
254  const_cast<MeshLib::Mesh&>(mesh), "strain_xx",
256  mesh_prop_epsilon_xx->resize(mesh.getNumberOfElements());
257  _process_data.mesh_prop_strain_xx = mesh_prop_epsilon_xx;
258 
259  auto mesh_prop_epsilon_yy = MeshLib::getOrCreateMeshProperty<double>(
260  const_cast<MeshLib::Mesh&>(mesh), "strain_yy",
262  mesh_prop_epsilon_yy->resize(mesh.getNumberOfElements());
263  _process_data.mesh_prop_strain_yy = mesh_prop_epsilon_yy;
264 
265  auto mesh_prop_epsilon_zz = MeshLib::getOrCreateMeshProperty<double>(
266  const_cast<MeshLib::Mesh&>(mesh), "strain_zz",
268  mesh_prop_epsilon_zz->resize(mesh.getNumberOfElements());
269  _process_data.mesh_prop_strain_zz = mesh_prop_epsilon_zz;
270 
271  auto mesh_prop_epsilon_xy = MeshLib::getOrCreateMeshProperty<double>(
272  const_cast<MeshLib::Mesh&>(mesh), "strain_xy",
274  mesh_prop_epsilon_xy->resize(mesh.getNumberOfElements());
275  _process_data.mesh_prop_strain_xy = mesh_prop_epsilon_xy;
276 
277  if (GlobalDim == 3)
278  {
279  auto mesh_prop_epsilon_xz = MeshLib::getOrCreateMeshProperty<double>(
280  const_cast<MeshLib::Mesh&>(mesh), "strain_xz",
282  mesh_prop_epsilon_xz->resize(mesh.getNumberOfElements());
283  _process_data.mesh_prop_strain_xz = mesh_prop_epsilon_xz;
284 
285  auto mesh_prop_epsilon_yz = MeshLib::getOrCreateMeshProperty<double>(
286  const_cast<MeshLib::Mesh&>(mesh), "strain_yz",
288  mesh_prop_epsilon_yz->resize(mesh.getNumberOfElements());
289  _process_data.mesh_prop_strain_yz = mesh_prop_epsilon_yz;
290  }
291 
292  auto mesh_prop_velocity = MeshLib::getOrCreateMeshProperty<double>(
293  const_cast<MeshLib::Mesh&>(mesh), "velocity",
294  MeshLib::MeshItemType::Cell, GlobalDim);
295  mesh_prop_velocity->resize(mesh.getNumberOfElements() * GlobalDim);
296  _process_data.mesh_prop_velocity = mesh_prop_velocity;
297 
298  if (!_vec_fracture_elements.empty())
299  {
300  auto mesh_prop_levelset = MeshLib::getOrCreateMeshProperty<double>(
301  const_cast<MeshLib::Mesh&>(mesh), "levelset1",
303  mesh_prop_levelset->resize(mesh.getNumberOfElements());
304  for (MeshLib::Element const* e : _mesh.getElements())
305  {
306  if (e->getDimension() < GlobalDim)
307  {
308  continue;
309  }
310 
311  std::vector<FractureProperty*> fracture_props(
312  {_process_data.fracture_property.get()});
313  std::vector<JunctionProperty*> junction_props;
314  std::unordered_map<int, int> fracID_to_local({{0, 0}});
315  std::vector<double> levelsets = uGlobalEnrichments(
316  fracture_props, junction_props, fracID_to_local,
317  Eigen::Vector3d(MeshLib::getCenterOfGravity(*e).data()));
318  (*mesh_prop_levelset)[e->getID()] = levelsets[0];
319  }
320 
321  auto mesh_prop_w_n = MeshLib::getOrCreateMeshProperty<double>(
322  const_cast<MeshLib::Mesh&>(mesh), "w_n",
324  mesh_prop_w_n->resize(mesh.getNumberOfElements());
325  auto mesh_prop_w_s = MeshLib::getOrCreateMeshProperty<double>(
326  const_cast<MeshLib::Mesh&>(mesh), "w_s",
328  mesh_prop_w_s->resize(mesh.getNumberOfElements());
329  _process_data.mesh_prop_w_n = mesh_prop_w_n;
330  _process_data.mesh_prop_w_s = mesh_prop_w_s;
331 
332  auto mesh_prop_b = MeshLib::getOrCreateMeshProperty<double>(
333  const_cast<MeshLib::Mesh&>(mesh), "aperture",
335  mesh_prop_b->resize(mesh.getNumberOfElements());
336  auto const mesh_prop_matid = materialIDs(mesh);
337  if (!mesh_prop_matid)
338  {
339  OGS_FATAL("Could not access MaterialIDs property from mesh.");
340  }
341  auto const& frac = _process_data.fracture_property;
342  for (MeshLib::Element const* e : _mesh.getElements())
343  {
344  if (e->getDimension() == GlobalDim)
345  {
346  continue;
347  }
348  if ((*mesh_prop_matid)[e->getID()] != frac->mat_id)
349  {
350  continue;
351  }
352  // Mean value for the element. This allows usage of node based
353  // properties for aperture.
354  (*mesh_prop_b)[e->getID()] =
355  frac->aperture0
356  .getNodalValuesOnElement(*e, /*time independent*/ 0)
357  .mean();
358  }
359  _process_data.mesh_prop_b = mesh_prop_b;
360 
361  auto mesh_prop_k_f = MeshLib::getOrCreateMeshProperty<double>(
362  const_cast<MeshLib::Mesh&>(mesh), "k_f",
364  mesh_prop_k_f->resize(mesh.getNumberOfElements());
365  _process_data.mesh_prop_k_f = mesh_prop_k_f;
366 
367  auto mesh_prop_fracture_stress_shear =
368  MeshLib::getOrCreateMeshProperty<double>(
369  const_cast<MeshLib::Mesh&>(mesh), "f_stress_s",
371  mesh_prop_fracture_stress_shear->resize(mesh.getNumberOfElements());
372  _process_data.mesh_prop_fracture_stress_shear =
373  mesh_prop_fracture_stress_shear;
374 
375  auto mesh_prop_fracture_stress_normal =
376  MeshLib::getOrCreateMeshProperty<double>(
377  const_cast<MeshLib::Mesh&>(mesh), "f_stress_n",
379  mesh_prop_fracture_stress_normal->resize(mesh.getNumberOfElements());
380  _process_data.mesh_prop_fracture_stress_normal =
381  mesh_prop_fracture_stress_normal;
382 
383  auto mesh_prop_fracture_shear_failure =
384  MeshLib::getOrCreateMeshProperty<double>(
385  const_cast<MeshLib::Mesh&>(mesh), "f_shear_failure",
387  mesh_prop_fracture_shear_failure->resize(mesh.getNumberOfElements());
388  _process_data.mesh_prop_fracture_shear_failure =
389  mesh_prop_fracture_shear_failure;
390 
391  auto mesh_prop_nodal_w = MeshLib::getOrCreateMeshProperty<double>(
392  const_cast<MeshLib::Mesh&>(mesh), "nodal_w",
393  MeshLib::MeshItemType::Node, GlobalDim);
394  mesh_prop_nodal_w->resize(mesh.getNumberOfNodes() * GlobalDim);
395  _process_data.mesh_prop_nodal_w = mesh_prop_nodal_w;
396 
397  auto mesh_prop_nodal_b = MeshLib::getOrCreateMeshProperty<double>(
398  const_cast<MeshLib::Mesh&>(mesh), "nodal_aperture",
400  mesh_prop_nodal_b->resize(mesh.getNumberOfNodes());
401  _process_data.mesh_prop_nodal_b = mesh_prop_nodal_b;
402 
403  if (GlobalDim == 3)
404  {
405  auto mesh_prop_w_s2 = MeshLib::getOrCreateMeshProperty<double>(
406  const_cast<MeshLib::Mesh&>(mesh), "w_s2",
408  mesh_prop_w_s2->resize(mesh.getNumberOfElements());
409  _process_data.mesh_prop_w_s2 = mesh_prop_w_s2;
410 
411  auto mesh_prop_fracture_stress_shear2 =
412  MeshLib::getOrCreateMeshProperty<double>(
413  const_cast<MeshLib::Mesh&>(mesh), "f_stress_s2",
415  mesh_prop_fracture_stress_shear2->resize(
416  mesh.getNumberOfElements());
417  _process_data.mesh_prop_fracture_stress_shear2 =
418  mesh_prop_fracture_stress_shear2;
419  }
420 
421  auto mesh_prop_nodal_p = MeshLib::getOrCreateMeshProperty<double>(
422  const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
424  mesh_prop_nodal_p->resize(mesh.getNumberOfNodes());
425  _process_data.mesh_prop_nodal_p = mesh_prop_nodal_p;
426 
427  _process_data.mesh_prop_nodal_forces =
428  MeshLib::getOrCreateMeshProperty<double>(
429  const_cast<MeshLib::Mesh&>(mesh), "NodalForces",
430  MeshLib::MeshItemType::Node, GlobalDim);
431  assert(_process_data.mesh_prop_nodal_forces->size() ==
432  GlobalDim * mesh.getNumberOfNodes());
433 
434  _process_data.mesh_prop_nodal_forces_jump =
435  MeshLib::getOrCreateMeshProperty<double>(
436  const_cast<MeshLib::Mesh&>(mesh), "NodalForcesJump",
437  MeshLib::MeshItemType::Node, GlobalDim);
438  assert(_process_data.mesh_prop_nodal_forces_jump->size() ==
439  GlobalDim * mesh.getNumberOfNodes());
440 
441  _process_data.mesh_prop_hydraulic_flow =
442  MeshLib::getOrCreateMeshProperty<double>(
443  const_cast<MeshLib::Mesh&>(mesh), "MassFlowRate",
445  assert(_process_data.mesh_prop_hydraulic_flow->size() ==
446  mesh.getNumberOfNodes());
447  }
448 }
449 
450 template <int GlobalDim>
452  std::vector<GlobalVector*> const& x, const double t, double const dt,
453  int const process_id)
454 {
455  if (process_id == 0)
456  {
457  DBUG("PostTimestep HydroMechanicsProcess.");
458  std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
459  auto const n_processes = x.size();
460  dof_tables.reserve(n_processes);
461  for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
462  {
463  dof_tables.push_back(&getDOFTable(process_id));
464  }
465 
466  ProcessLib::ProcessVariable const& pv =
467  getProcessVariables(process_id)[0];
470  _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, t, dt);
471  }
472 
473  DBUG("Compute the secondary variables for HydroMechanicsProcess.");
474 
475  const auto& dof_table = getDOFTable(process_id);
476 
477  // Copy displacement jumps in a solution vector to mesh property
478  // Remark: the copy is required because mesh properties for primary
479  // variables are set during output and are not ready yet when this function
480  // is called.
481  int g_variable_id = 0;
482  {
483  const int monolithic_process_id = 0;
484  auto const& pvs = getProcessVariables(monolithic_process_id);
485  auto const it =
486  std::find_if(pvs.begin(), pvs.end(),
487  [](ProcessVariable const& pv)
488  { return pv.getName() == "displacement_jump1"; });
489  if (it == pvs.end())
490  {
491  OGS_FATAL(
492  "Didn't find expected 'displacement_jump1' process variable.");
493  }
494  g_variable_id = static_cast<int>(std::distance(pvs.begin(), it));
495  }
496 
498 
499  const int monolithic_process_id = 0;
500  ProcessVariable& pv_g =
501  this->getProcessVariables(monolithic_process_id)[g_variable_id];
502  auto const num_comp = pv_g.getNumberOfGlobalComponents();
503  auto& mesh_prop_g = *MeshLib::getOrCreateMeshProperty<double>(
504  _mesh, pv_g.getName(), MeshLib::MeshItemType::Node, num_comp);
505  for (int component_id = 0; component_id < num_comp; ++component_id)
506  {
507  auto const& mesh_subset =
508  dof_table.getMeshSubset(g_variable_id, component_id);
509  auto const mesh_id = mesh_subset.getMeshID();
510  for (auto const* node : mesh_subset.getNodes())
511  {
513  node->getID());
514 
515  auto const global_index =
516  dof_table.getGlobalIndex(l, g_variable_id, component_id);
517  mesh_prop_g[node->getID() * num_comp + component_id] =
518  (*x[process_id])[global_index];
519  }
520  }
521 
522  // compute nodal w and aperture
523  auto const& R = _process_data.fracture_property->R;
524  auto const& b0 = _process_data.fracture_property->aperture0;
525  MeshLib::PropertyVector<double>& vec_w = *_process_data.mesh_prop_nodal_w;
526  MeshLib::PropertyVector<double>& vec_b = *_process_data.mesh_prop_nodal_b;
527 
528  auto compute_nodal_aperture =
529  [&](std::size_t const node_id, double const w_n)
530  {
531  // skip aperture computation for element-wise defined b0 because there
532  // are jumps on the nodes between the element's values.
533  if (dynamic_cast<ParameterLib::MeshElementParameter<double> const*>(
534  &b0))
535  {
536  return std::numeric_limits<double>::quiet_NaN();
537  }
538 
540  x.setNodeID(node_id);
541  return w_n + b0(/*time independent*/ 0, x)[0];
542  };
543 
544  Eigen::VectorXd g(GlobalDim);
545  Eigen::VectorXd w(GlobalDim);
546  for (MeshLib::Node const* node : _vec_fracture_nodes)
547  {
548  auto const node_id = node->getID();
549  g.setZero();
550  for (int k = 0; k < GlobalDim; k++)
551  {
552  g[k] = mesh_prop_g[node_id * GlobalDim + k];
553  }
554 
555  w.noalias() = R * g;
556  for (int k = 0; k < GlobalDim; k++)
557  {
558  vec_w[node_id * GlobalDim + k] = w[k];
559  }
560 
561  vec_b[node_id] = compute_nodal_aperture(node_id, w[GlobalDim - 1]);
562  }
563 }
564 
565 template <int GlobalDim>
567 {
568  return false;
569 }
570 
571 template <int GlobalDim>
573  const double t, double const dt, std::vector<GlobalVector*> const& x,
574  std::vector<GlobalVector*> const& xdot, int const process_id,
576 {
577  DBUG("Assemble HydroMechanicsProcess.");
578 
579  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
580  dof_table = {std::ref(*_local_to_global_index_map)};
581  // Call global assembler for each local assembly item.
583  _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
584  dof_table, t, dt, x, xdot, process_id, M, K, b);
585 }
586 
587 template <int GlobalDim>
589  const double t, double const dt, std::vector<GlobalVector*> const& x,
590  std::vector<GlobalVector*> const& xdot, int const process_id,
592 {
593  DBUG("AssembleWithJacobian HydroMechanicsProcess.");
594 
595  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
596 
597  // Call global assembler for each local assembly item.
598  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
599  dof_table = {std::ref(*_local_to_global_index_map)};
602  _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
603  process_id, M, K, b, Jac);
604 
605  auto copyRhs = [&](int const variable_id, auto& output_vector)
606  {
607  transformVariableFromGlobalVector(b, variable_id,
608  *_local_to_global_index_map,
609  output_vector, std::negate<double>());
610  };
611  copyRhs(0, *_process_data.mesh_prop_hydraulic_flow);
612  copyRhs(1, *_process_data.mesh_prop_nodal_forces);
613  copyRhs(2, *_process_data.mesh_prop_nodal_forces_jump);
614 }
615 
616 template <int GlobalDim>
618  std::vector<GlobalVector*> const& x, double const t, double const dt,
619  int const process_id)
620 {
621  DBUG("PreTimestep HydroMechanicsProcess.");
622 
623  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
624 
627  pv.getActiveElementIDs(), *_local_to_global_index_map, *x[process_id],
628  t, dt);
629 }
630 
631 // ------------------------------------------------------------------------------------
632 // template instantiation
633 // ------------------------------------------------------------------------------------
634 template class HydroMechanicsProcess<2>;
635 template class HydroMechanicsProcess<3>;
636 
637 } // namespace HydroMechanics
638 } // namespace LIE
639 } // namespace ProcessLib
Definition of the ElementStatus class.
#define OGS_FATAL(...)
Definition: Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:34
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:29
Definition of the MeshInformation class.
Definition of the class Properties that implements a container of properties.
Definition of the Mesh class.
Global vector based on Eigen vector.
Definition: EigenVector.h:28
const double * data() const
Definition: Point3d.h:63
static std::optional< std::pair< T, T > > const getValueBounds(PropertyVector< T > const &property)
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
const std::string getName() const
Get name of the mesh.
Definition: Mesh.h:97
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:94
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:91
void setNodeID(std::size_t node_id)
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, int const process_id) override
std::vector< MeshLib::Element * > _vec_fracture_matrix_elements
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
HydroMechanicsProcess(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, HydroMechanicsProcessData< GlobalDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, int const process_id) override
HydroMechanicsProcessData< GlobalDim > _process_data
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
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
virtual void postTimestep(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
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)
std::string const & getName() const
ParameterLib::Parameter< double > const & getInitialCondition() const
int getNumberOfGlobalComponents() const
Returns the number of components of the process variable.
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:148
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)
void setLocalAccessibleVector(PETScVector const &x)
Definition: LinAlg.cpp:27
static const double t
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition: Mesh.cpp:264
MathLib::Point3d getCenterOfGravity(Element const &element)
Calculates the center of gravity for the mesh element.
Definition: Element.cpp:124
std::vector< Node * > getBaseNodes(std::vector< Element * > const &elements)
Definition: Utils.h:26
@ BY_COMPONENT
Ordering data by component type.
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 createLocalAssemblers(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)
void setFractureProperty(int const dim, MeshLib::Element const &e, FractureProperty &frac_prop)
std::vector< double > uGlobalEnrichments(std::vector< FractureProperty * > const &frac_props, std::vector< JunctionProperty * > const &junction_props, std::unordered_map< int, int > const &fracID_to_local, Eigen::Vector3d const &x)
void getFractureMatrixDataInMesh(MeshLib::Mesh const &mesh, std::vector< MeshLib::Element * > &vec_matrix_elements, std::vector< int > &vec_fracture_mat_IDs, std::vector< std::vector< MeshLib::Element * >> &vec_fracture_elements, std::vector< std::vector< MeshLib::Element * >> &vec_fracture_matrix_elements, std::vector< std::vector< MeshLib::Node * >> &vec_fracture_nodes, std::vector< std::pair< std::size_t, std::vector< int >>> &vec_branch_nodeID_matIDs, std::vector< std::pair< std::size_t, std::vector< int >>> &vec_junction_nodeID_matIDs)
Definition: MeshUtils.cpp:213
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 executeMemberDereferenced(Object &object, Method method, Container const &container, Args &&... args)
A parameter represented by a mesh property vector.