13 #include "LocalAssembler/CreateLocalAssemblers.h"
34 namespace HydroMechanics
36 template <
int GlobalDim>
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>>>&&
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))
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;
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());
97 std::make_unique<MeshLib::ElementStatus>(&mesh);
106 "Could not get minimum/maximum ranges values for the "
107 "MaterialIDs property in the mesh '{:s}'.",
111 std::vector<int> vec_p_inactive_matIDs;
112 for (
int matID = range->first; matID <= range->second; matID++)
114 if (std::find(vec_fracture_mat_IDs.begin(),
115 vec_fracture_mat_IDs.end(),
116 matID) == vec_fracture_mat_IDs.end())
118 vec_p_inactive_matIDs.push_back(matID);
122 std::make_unique<MeshLib::ElementStatus>(&mesh,
123 vec_p_inactive_matIDs);
125 const int monolithic_process_id = 0;
132 template <
int GlobalDim>
139 _mesh_subset_all_nodes =
140 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
143 _process_data.p_element_status->getActiveElements());
144 _mesh_subset_nodes_p =
145 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh_nodes_p);
147 _mesh_subset_matrix_nodes =
148 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
149 if (!_vec_fracture_nodes.empty())
152 _mesh_subset_fracture_nodes =
153 std::make_unique<MeshLib::MeshSubset>(_mesh, _vec_fracture_nodes);
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;
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)
165 vec_var_elements.push_back(&_mesh.getElements());
171 vec_var_elements.push_back(&_vec_fracture_matrix_elements);
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())
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);
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),
195 DBUG(
"created {:d} DoF", _local_to_global_index_map->size());
198 template <
int GlobalDim>
202 unsigned const integration_order)
205 INFO(
"[LIE/HM] creating local assemblers");
206 const int monolithic_process_id = 0;
213 getProcessVariables(monolithic_process_id)[1]
215 .getShapeFunctionOrder(),
219 auto mesh_prop_sigma_xx = MeshLib::getOrCreateMeshProperty<double>(
223 _process_data.mesh_prop_stress_xx = mesh_prop_sigma_xx;
225 auto mesh_prop_sigma_yy = MeshLib::getOrCreateMeshProperty<double>(
229 _process_data.mesh_prop_stress_yy = mesh_prop_sigma_yy;
231 auto mesh_prop_sigma_zz = MeshLib::getOrCreateMeshProperty<double>(
235 _process_data.mesh_prop_stress_zz = mesh_prop_sigma_zz;
237 auto mesh_prop_sigma_xy = MeshLib::getOrCreateMeshProperty<double>(
241 _process_data.mesh_prop_stress_xy = mesh_prop_sigma_xy;
245 auto mesh_prop_sigma_xz = MeshLib::getOrCreateMeshProperty<double>(
249 _process_data.mesh_prop_stress_xz = mesh_prop_sigma_xz;
251 auto mesh_prop_sigma_yz = MeshLib::getOrCreateMeshProperty<double>(
255 _process_data.mesh_prop_stress_yz = mesh_prop_sigma_yz;
258 auto mesh_prop_epsilon_xx = MeshLib::getOrCreateMeshProperty<double>(
262 _process_data.mesh_prop_strain_xx = mesh_prop_epsilon_xx;
264 auto mesh_prop_epsilon_yy = MeshLib::getOrCreateMeshProperty<double>(
268 _process_data.mesh_prop_strain_yy = mesh_prop_epsilon_yy;
270 auto mesh_prop_epsilon_zz = MeshLib::getOrCreateMeshProperty<double>(
274 _process_data.mesh_prop_strain_zz = mesh_prop_epsilon_zz;
276 auto mesh_prop_epsilon_xy = MeshLib::getOrCreateMeshProperty<double>(
280 _process_data.mesh_prop_strain_xy = mesh_prop_epsilon_xy;
284 auto mesh_prop_epsilon_xz = MeshLib::getOrCreateMeshProperty<double>(
288 _process_data.mesh_prop_strain_xz = mesh_prop_epsilon_xz;
290 auto mesh_prop_epsilon_yz = MeshLib::getOrCreateMeshProperty<double>(
294 _process_data.mesh_prop_strain_yz = mesh_prop_epsilon_yz;
297 auto mesh_prop_velocity = MeshLib::getOrCreateMeshProperty<double>(
301 _process_data.mesh_prop_velocity = mesh_prop_velocity;
303 if (!_vec_fracture_elements.empty())
305 auto mesh_prop_levelset = MeshLib::getOrCreateMeshProperty<double>(
311 if (e->getDimension() < GlobalDim)
316 std::vector<FractureProperty*> fracture_props(
317 {_process_data.fracture_property.get()});
318 std::vector<JunctionProperty*> junction_props;
319 std::unordered_map<int, int> fracID_to_local({{0, 0}});
321 fracture_props, junction_props, fracID_to_local,
323 (*mesh_prop_levelset)[e->getID()] = levelsets[0];
326 auto mesh_prop_w_n = MeshLib::getOrCreateMeshProperty<double>(
330 auto mesh_prop_w_s = MeshLib::getOrCreateMeshProperty<double>(
334 _process_data.mesh_prop_w_n = mesh_prop_w_n;
335 _process_data.mesh_prop_w_s = mesh_prop_w_s;
337 auto mesh_prop_b = MeshLib::getOrCreateMeshProperty<double>(
342 if (!mesh_prop_matid)
344 OGS_FATAL(
"Could not access MaterialIDs property from mesh.");
346 auto const& frac = _process_data.fracture_property;
349 if (e->getDimension() == GlobalDim)
353 if ((*mesh_prop_matid)[e->getID()] != frac->mat_id)
359 (*mesh_prop_b)[e->getID()] =
361 .getNodalValuesOnElement(*e, 0)
364 _process_data.mesh_prop_b = mesh_prop_b;
366 auto mesh_prop_k_f = MeshLib::getOrCreateMeshProperty<double>(
370 _process_data.mesh_prop_k_f = mesh_prop_k_f;
372 auto mesh_prop_fracture_stress_shear =
373 MeshLib::getOrCreateMeshProperty<double>(
377 _process_data.mesh_prop_fracture_stress_shear =
378 mesh_prop_fracture_stress_shear;
380 auto mesh_prop_fracture_stress_normal =
381 MeshLib::getOrCreateMeshProperty<double>(
385 _process_data.mesh_prop_fracture_stress_normal =
386 mesh_prop_fracture_stress_normal;
388 auto mesh_prop_fracture_shear_failure =
389 MeshLib::getOrCreateMeshProperty<double>(
393 _process_data.mesh_prop_fracture_shear_failure =
394 mesh_prop_fracture_shear_failure;
396 auto mesh_prop_nodal_w = MeshLib::getOrCreateMeshProperty<double>(
400 _process_data.mesh_prop_nodal_w = mesh_prop_nodal_w;
402 auto mesh_prop_nodal_b = MeshLib::getOrCreateMeshProperty<double>(
406 _process_data.mesh_prop_nodal_b = mesh_prop_nodal_b;
410 auto mesh_prop_w_s2 = MeshLib::getOrCreateMeshProperty<double>(
414 _process_data.mesh_prop_w_s2 = mesh_prop_w_s2;
416 auto mesh_prop_fracture_stress_shear2 =
417 MeshLib::getOrCreateMeshProperty<double>(
420 mesh_prop_fracture_stress_shear2->resize(
422 _process_data.mesh_prop_fracture_stress_shear2 =
423 mesh_prop_fracture_stress_shear2;
426 auto mesh_prop_nodal_p = MeshLib::getOrCreateMeshProperty<double>(
430 _process_data.mesh_prop_nodal_p = mesh_prop_nodal_p;
432 _process_data.mesh_prop_nodal_forces =
433 MeshLib::getOrCreateMeshProperty<double>(
436 assert(_process_data.mesh_prop_nodal_forces->size() ==
439 _process_data.mesh_prop_nodal_forces_jump =
440 MeshLib::getOrCreateMeshProperty<double>(
443 assert(_process_data.mesh_prop_nodal_forces_jump->size() ==
446 _process_data.mesh_prop_hydraulic_flow =
447 MeshLib::getOrCreateMeshProperty<double>(
450 assert(_process_data.mesh_prop_hydraulic_flow->size() ==
455 template <
int GlobalDim>
457 std::vector<GlobalVector*>
const& x,
const double t,
double const dt,
458 int const process_id)
462 DBUG(
"PostTimestep HydroMechanicsProcess.");
463 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
464 auto const n_processes = x.size();
465 dof_tables.reserve(n_processes);
466 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
468 dof_tables.push_back(&getDOFTable(process_id));
472 getProcessVariables(process_id)[0];
478 DBUG(
"Compute the secondary variables for HydroMechanicsProcess.");
480 const auto& dof_table = getDOFTable(process_id);
486 int g_variable_id = 0;
488 const int monolithic_process_id = 0;
489 auto const& pvs = getProcessVariables(monolithic_process_id);
491 std::find_if(pvs.begin(), pvs.end(),
493 { return pv.getName() ==
"displacement_jump1"; });
497 "Didn't find expected 'displacement_jump1' process variable.");
499 g_variable_id =
static_cast<int>(std::distance(pvs.begin(), it));
504 const int monolithic_process_id = 0;
506 this->getProcessVariables(monolithic_process_id)[g_variable_id];
508 auto& mesh_prop_g = *MeshLib::getOrCreateMeshProperty<double>(
510 for (
int component_id = 0; component_id < num_comp; ++component_id)
512 auto const& mesh_subset =
513 dof_table.getMeshSubset(g_variable_id, component_id);
514 auto const mesh_id = mesh_subset.getMeshID();
515 for (
auto const* node : mesh_subset.getNodes())
520 auto const global_index =
521 dof_table.getGlobalIndex(l, g_variable_id, component_id);
522 mesh_prop_g[node->getID() * num_comp + component_id] =
523 (*x[process_id])[global_index];
528 auto const& R = _process_data.fracture_property->R;
529 auto const& b0 = _process_data.fracture_property->aperture0;
533 auto compute_nodal_aperture =
534 [&](std::size_t
const node_id,
double const w_n)
541 return std::numeric_limits<double>::quiet_NaN();
546 return w_n + b0( 0, x)[0];
549 Eigen::VectorXd g(GlobalDim);
550 Eigen::VectorXd w(GlobalDim);
553 auto const node_id = node->getID();
555 for (
int k = 0; k < GlobalDim; k++)
557 g[k] = mesh_prop_g[node_id * GlobalDim + k];
561 for (
int k = 0; k < GlobalDim; k++)
563 vec_w[node_id * GlobalDim + k] = w[k];
566 vec_b[node_id] = compute_nodal_aperture(node_id, w[GlobalDim - 1]);
570 template <
int GlobalDim>
576 template <
int GlobalDim>
578 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
579 std::vector<GlobalVector*>
const& xdot,
int const process_id,
582 DBUG(
"Assemble HydroMechanicsProcess.");
584 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
585 dof_table = {std::ref(*_local_to_global_index_map)};
589 dof_table, t, dt, x, xdot, process_id, M, K, b);
592 template <
int GlobalDim>
594 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
595 std::vector<GlobalVector*>
const& xdot,
const double dxdot_dx,
599 DBUG(
"AssembleWithJacobian HydroMechanicsProcess.");
604 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
605 dof_table = {std::ref(*_local_to_global_index_map)};
609 dxdot_dx, dx_dx, process_id, M, K, b, Jac);
611 auto copyRhs = [&](
int const variable_id,
auto& output_vector)
614 *_local_to_global_index_map,
615 output_vector, std::negate<double>());
617 copyRhs(0, *_process_data.mesh_prop_hydraulic_flow);
618 copyRhs(1, *_process_data.mesh_prop_nodal_forces);
619 copyRhs(2, *_process_data.mesh_prop_nodal_forces_jump);
622 template <
int GlobalDim>
624 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
625 int const process_id)
627 DBUG(
"PreTimestep HydroMechanicsProcess.");
Definition of the ElementStatus class.
void INFO(char const *fmt, Args const &... args)
void DBUG(char const *fmt, Args const &... args)
Definition of the class Properties that implements a container of properties.
Definition of the Mesh class.
Global vector based on Eigen vector.
const T * getCoords() const
bool isAxiallySymmetric() const
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
const std::string getName() const
Get name of the mesh.
std::size_t getNumberOfNodes() const
Get the number of nodes.
std::size_t getNumberOfElements() const
Get the number of elements.
void setNodeID(std::size_t node_id)
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, const double dxdot_dx, const double dx_dx, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, int const process_id) override
std::vector< MeshLib::Node * > _vec_fracture_nodes
std::vector< MeshLib::Element * > _vec_matrix_elements
bool isLinear() const override
std::vector< MeshLib::Element * > _vec_fracture_matrix_elements
std::vector< MeshLib::Element * > _vec_fracture_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 ¶meters, 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 constructDofTable() 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
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, const double dxdot_dx, const double dx_dx, 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)
PropertyVector< int > const * materialIDs(Mesh const &mesh)
MeshLib::Node getCenterOfGravity(Element const &element)
Calculates the center of gravity for the mesh element.
std::vector< Node * > getBaseNodes(std::vector< Element * > const &elements)
@ 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)
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)
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)
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.