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