13#include <range/v3/view/join.hpp>
37namespace HydroMechanics
39template <
int DisplacementDim>
43 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
44 std::vector<std::unique_ptr<ParameterLib::ParameterBase>>
const& parameters,
45 unsigned const integration_order,
46 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
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), use_monolithic_scheme),
54 _process_data(std::move(process_data))
56 INFO(
"[LIE/HM] looking for fracture elements in the given mesh");
57 std::vector<std::pair<std::size_t, std::vector<int>>>
58 vec_branch_nodeID_matIDs;
59 std::vector<std::pair<std::size_t, std::vector<int>>>
60 vec_junction_nodeID_matIDs;
65 vec_junction_nodeID_matIDs);
71 "The number of the given fracture properties ({:d}) are not "
72 "consistent with the number of fracture groups in a mesh ({:d}).",
80 _process_data.map_materialID_to_fractureID.resize(*max_frac_mat_id + 1);
94 _process_data.vec_ele_connected_fractureIDs[e->getID()].push_back(
100 for (
auto& fracture_prop :
_process_data.fracture_properties)
110 for (
auto const& [vec_branch_nodeID, matID] : vec_branch_nodeID_matIDs)
112 auto const master_matId = matID[0];
113 auto const slave_matId = matID[1];
122 *mesh.
getNode(vec_branch_nodeID), master_frac, slave_frac));
125 *mesh.
getNode(vec_branch_nodeID), master_frac, slave_frac));
129 transform(cbegin(vec_junction_nodeID_matIDs),
130 cend(vec_junction_nodeID_matIDs),
132 [&](
auto& vec_junction_nodeID_matID)
138 for (std::size_t i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
140 auto const& material_ids = vec_junction_nodeID_matIDs[i].second;
141 assert(material_ids.size() == 2);
142 std::array<int, 2> fracture_ids{
143 {
_process_data.map_materialID_to_fractureID[material_ids[0]],
144 _process_data.map_materialID_to_fractureID[material_ids[1]]}};
147 i, *mesh.
getNode(vec_junction_nodeID_matIDs[i].first),
154 for (
unsigned i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
156 auto node = mesh.
getNode(vec_junction_nodeID_matIDs[i].first);
160 _process_data.vec_ele_connected_junctionIDs[id].push_back(i);
166 vec_junction_nodeID_matIDs.size());
167 for (
unsigned i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
169 auto node = mesh.
getNode(vec_junction_nodeID_matIDs[i].first);
187 std::make_unique<MeshLib::ElementStatus>(&mesh);
196 "Could not get minimum/maximum ranges values for the "
197 "MaterialIDs property in the mesh '{:s}'.",
201 std::vector<int> vec_p_inactive_matIDs;
202 for (
int matID = range->first; matID <= range->second; matID++)
208 vec_p_inactive_matIDs.push_back(matID);
212 std::make_unique<MeshLib::ElementStatus>(&mesh,
213 vec_p_inactive_matIDs);
215 const int monolithic_process_id = 0;
225template <
int DisplacementDim>
232 _mesh_subset_all_nodes =
233 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
236 _process_data.p_element_status->getActiveElements());
237 _mesh_subset_nodes_p =
238 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh_nodes_p);
240 _mesh_subset_matrix_nodes =
241 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
243 for (
unsigned i = 0; i < _vec_fracture_nodes.size(); i++)
245 _mesh_subset_fracture_nodes.push_back(
246 std::make_unique<MeshLib::MeshSubset const>(
247 _mesh, _vec_fracture_nodes[i]));
250 _mesh_subset_junction_nodes =
251 std::make_unique<MeshLib::MeshSubset>(_mesh, _vec_junction_nodes);
255 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
256 all_mesh_subsets.emplace_back(*_mesh_subset_nodes_p);
257 std::generate_n(std::back_inserter(all_mesh_subsets), DisplacementDim,
258 [&]() {
return *_mesh_subset_matrix_nodes; });
259 for (
auto const& ms : _mesh_subset_fracture_nodes)
261 std::generate_n(std::back_inserter(all_mesh_subsets),
263 [&]() {
return *ms; });
265 std::generate_n(std::back_inserter(all_mesh_subsets),
267 [&]() {
return *_mesh_subset_junction_nodes; });
271 std::vector<int> vec_n_components;
272 vec_n_components.push_back(1);
273 vec_n_components.insert(
274 vec_n_components.end(),
275 1 + _vec_fracture_mat_IDs.size() + _vec_junction_nodes.size(),
278 auto const all_fracture_matrix_elements = _vec_fracture_matrix_elements |
279 ranges::views::join |
280 ranges::to<std::vector>();
281 std::vector<std::vector<MeshLib::Element*>
const*> vec_var_elements;
282 if (!_process_data.deactivate_matrix_in_flow)
284 vec_var_elements.push_back(&_mesh.getElements());
290 vec_var_elements.push_back(&all_fracture_matrix_elements);
292 vec_var_elements.push_back(&_vec_matrix_elements);
293 for (
unsigned i = 0; i < _vec_fracture_matrix_elements.size(); i++)
295 vec_var_elements.push_back(&_vec_fracture_matrix_elements[i]);
297 for (
unsigned i = 0; i < _vec_junction_fracture_matrix_elements.size(); i++)
299 vec_var_elements.push_back(&_vec_junction_fracture_matrix_elements[i]);
302 INFO(
"[LIE/HM] creating a DoF table");
303 _local_to_global_index_map =
304 std::make_unique<NumLib::LocalToGlobalIndexMap>(
305 std::move(all_mesh_subsets),
310 DBUG(
"[LIE/HM] created {:d} DoF", _local_to_global_index_map->size());
313template <
int DisplacementDim>
317 Eigen::Vector3d
const pt(getCenterOfGravity(e).asEigenVector3d());
318 std::vector<FractureProperty*> e_fracture_props;
319 std::unordered_map<int, int> e_fracID_to_local;
321 for (
auto fid : _process_data.vec_ele_connected_fractureIDs[e.
getID()])
323 e_fracture_props.push_back(&_process_data.fracture_properties[fid]);
324 e_fracID_to_local.insert({fid, tmpi++});
326 std::vector<JunctionProperty*> e_junction_props;
327 std::unordered_map<int, int> e_juncID_to_local;
329 for (
auto fid : _process_data.vec_ele_connected_junctionIDs[e.
getID()])
331 e_junction_props.push_back(&_process_data.junction_properties[fid]);
332 e_juncID_to_local.insert({fid, tmpi++});
335 e_fracture_props, e_junction_props, e_fracID_to_local, pt));
337 auto update_levelset_property = [&](
unsigned const i,
int const id,
338 unsigned const levelset_idx_offset,
339 unsigned const name_offset)
343 "levelset" + std::to_string(
id + 1 + name_offset),
346 (*levelset_property)[e.
getID()] = levelsets[i + levelset_idx_offset];
349 for (
unsigned i = 0; i < e_fracture_props.size(); i++)
351 update_levelset_property(i, e_fracture_props[i]->fracture_id, 0, 0);
353 for (
unsigned i = 0; i < e_junction_props.size(); i++)
355 update_levelset_property(i, e_junction_props[i]->junction_id,
356 e_fracture_props.size(),
357 _process_data.fracture_properties.size());
361template <
int DisplacementDim>
365 unsigned const integration_order)
368 INFO(
"[LIE/HM] creating local assemblers");
377 auto add_secondary_variable = [&](std::string
const& name,
378 int const num_components,
379 auto get_ip_values_function)
381 _secondary_variables.addSecondaryVariable(
385 std::move(get_ip_values_function)));
388 add_secondary_variable(
"sigma",
390 DisplacementDim>::RowsAtCompileTime,
391 &LocalAssemblerInterface::getIntPtSigma);
393 add_secondary_variable(
"epsilon",
395 DisplacementDim>::RowsAtCompileTime,
396 &LocalAssemblerInterface::getIntPtEpsilon);
398 add_secondary_variable(
"velocity", DisplacementDim,
399 &LocalAssemblerInterface::getIntPtDarcyVelocity);
401 add_secondary_variable(
"fracture_velocity", DisplacementDim,
402 &LocalAssemblerInterface::getIntPtFractureVelocity);
404 add_secondary_variable(
"fracture_stress", DisplacementDim,
405 &LocalAssemblerInterface::getIntPtFractureStress);
407 add_secondary_variable(
"fracture_aperture", 1,
408 &LocalAssemblerInterface::getIntPtFractureAperture);
410 add_secondary_variable(
411 "fracture_permeability", 1,
412 &LocalAssemblerInterface::getIntPtFracturePermeability);
418 DisplacementDim>::RowsAtCompileTime);
426 if (e->getDimension() < DisplacementDim)
431 updateElementLevelSets(*e,
const_cast<MeshLib::Mesh&
>(mesh));
434 _process_data.element_local_jumps =
439 _process_data.element_fracture_stresses =
444 _process_data.element_fracture_velocities =
454 auto const& mesh_prop_matid = *_process_data.mesh_prop_materialIDs;
455 for (
auto const& fracture_prop : _process_data.fracture_properties)
459 if (e->getDimension() == DisplacementDim)
463 if (mesh_prop_matid[e->getID()] != fracture_prop.mat_id)
469 (*mesh_prop_b)[e->getID()] =
470 fracture_prop.aperture0
471 .getNodalValuesOnElement(*e, 0)
477 const_cast<MeshLib::Mesh&
>(mesh),
"fracture_permeability_avg",
480 _process_data.mesh_prop_k_f = mesh_prop_k_f;
482 auto mesh_prop_fracture_shear_failure =
487 _process_data.mesh_prop_fracture_shear_failure =
488 mesh_prop_fracture_shear_failure;
494 _process_data.mesh_prop_nodal_p = mesh_prop_nodal_p;
496 _process_data.mesh_prop_nodal_forces =
500 assert(_process_data.mesh_prop_nodal_forces->size() ==
503 _process_data.mesh_prop_nodal_forces_jump =
507 assert(_process_data.mesh_prop_nodal_forces_jump->size() ==
510 _process_data.mesh_prop_hydraulic_flow =
514 assert(_process_data.mesh_prop_hydraulic_flow->size() ==
516 _process_data.mesh_prop_b = mesh_prop_b;
519template <
int DisplacementDim>
521 std::vector<GlobalVector*>
const& x,
522 std::vector<GlobalVector*>
const& x_prev,
const double t,
double const dt,
523 int const process_id)
527 DBUG(
"PostTimestep HydroMechanicsProcess.");
531 _local_assemblers, getActiveElementIDs(), getDOFTables(x.size()), x,
532 x_prev, t, dt, process_id);
535 DBUG(
"Compute the secondary variables for HydroMechanicsProcess.");
537 const auto& dof_table = getDOFTable(process_id);
543 int g_variable_id = 0;
545 const int monolithic_process_id = 0;
546 auto const& pvs = getProcessVariables(monolithic_process_id);
548 std::find_if(pvs.begin(), pvs.end(),
550 { return pv.getName() ==
"displacement_jump1"; });
554 "Didn't find expected 'displacement_jump1' process variable.");
556 g_variable_id =
static_cast<int>(std::distance(pvs.begin(), it));
561 const int monolithic_process_id = 0;
563 this->getProcessVariables(monolithic_process_id)[g_variable_id];
567 for (
int component_id = 0; component_id < num_comp; ++component_id)
569 auto const& mesh_subset =
570 dof_table.getMeshSubset(g_variable_id, component_id);
574 auto const global_index =
575 dof_table.getGlobalIndex(l, g_variable_id, component_id);
576 auto const node_id = l.item_id;
577 mesh_prop_g[node_id * num_comp + component_id] =
578 (*x[process_id])[global_index];
583template <
int DisplacementDim>
589template <
int DisplacementDim>
591 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
592 std::vector<GlobalVector*>
const& x_prev,
int const process_id,
595 DBUG(
"Assemble HydroMechanicsProcess.");
597 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
598 _local_to_global_index_map.get()};
602 dof_table, t, dt, x, x_prev, process_id, &M, &K, &b);
605template <
int DisplacementDim>
608 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
609 std::vector<GlobalVector*>
const& x_prev,
int const process_id,
612 DBUG(
"AssembleWithJacobian HydroMechanicsProcess.");
615 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
616 _local_to_global_index_map.get()};
619 _local_assemblers, getActiveElementIDs(), dof_table, t, dt, x, x_prev,
620 process_id, &b, &Jac);
622 auto copyRhs = [&](
int const variable_id,
auto& output_vector)
624 transformVariableFromGlobalVector(b, variable_id,
625 *_local_to_global_index_map,
626 output_vector, std::negate<double>());
628 copyRhs(0, *_process_data.mesh_prop_hydraulic_flow);
629 copyRhs(1, *_process_data.mesh_prop_nodal_forces);
630 copyRhs(2, *_process_data.mesh_prop_nodal_forces_jump);
633template <
int DisplacementDim>
635 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
636 int const process_id)
638 DBUG(
"PreTimestep HydroMechanicsProcess.");
642 getActiveElementIDs(), *_local_to_global_index_map, *x[process_id], t,
Definition of the ElementStatus class.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the class Properties that implements a container of properties.
Definition of the Mesh class.
Global vector based on Eigen vector.
std::size_t getID() const
Returns the ID of the element.
bool isAxiallySymmetric() const
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
const Node * getNode(std::size_t idx) const
Get the node with the given index.
Properties & getProperties()
const std::string getName() const
Get name of the mesh.
std::size_t getNumberOfNodes() const
Get the number of nodes.
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
std::size_t getNumberOfElements() const
Get the number of elements.
PropertyVector< T > const * getPropertyVector(std::string_view name) const
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< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id) override
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
std::vector< std::vector< MeshLib::Node * > > _vec_fracture_nodes
std::vector< MeshLib::Element * > _vec_matrix_elements
std::vector< int > _vec_fracture_mat_IDs
std::vector< std::vector< MeshLib::Element * > > _vec_junction_fracture_matrix_elements
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, int const process_id) override
std::vector< MeshLib::Node * > _vec_junction_nodes
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
void constructDofTable() override
void updateElementLevelSets(MeshLib::Element const &e, MeshLib::Mesh &mesh)
std::vector< std::vector< MeshLib::Element * > > _vec_fracture_elements
HydroMechanicsProcessData< DisplacementDim > _process_data
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector &b, GlobalMatrix &Jac) override
bool isLinear() const override
std::vector< std::vector< MeshLib::Element * > > _vec_fracture_matrix_elements
virtual void postTimestep(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
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)
ParameterLib::Parameter< double > const & getInitialCondition() const
std::string const & getName() const
int getNumberOfGlobalComponents() const
Returns the number of components of the process variable.
std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > const & getProcessVariables() const
Handles configuration of several secondary variables from the project file.
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix *M, GlobalMatrix *K, GlobalVector *b)
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector *b, GlobalMatrix *Jac)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
void setLocalAccessibleVector(PETScVector const &x)
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
auto meshLocations(Mesh const &mesh, MeshItemType const item_type)
PropertyVector< T > * getOrCreateMeshProperty(Mesh &mesh, std::string const &property_name, MeshItemType const item_type, int const number_of_components)
std::vector< Node * > getBaseNodes(std::vector< Element * > const &elements)
@ BY_COMPONENT
Ordering data by component type.
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)
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)
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)
BranchProperty createBranchProperty(MeshLib::Node const &branchNode, FractureProperty const &master_frac, FractureProperty const &slave_frac)
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 executeMemberDereferenced(Object &object, Method method, Container const &container, Args &&... args)