24 namespace HydroMechanics
26 template <
int DisplacementDim>
30 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
31 std::vector<std::unique_ptr<ParameterLib::ParameterBase>>
const& parameters,
32 unsigned const integration_order,
33 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
37 bool const use_monolithic_scheme)
38 :
Process(std::move(
name), mesh, std::move(jacobian_assembler), parameters,
39 integration_order, std::move(process_variables),
40 std::move(secondary_variables), use_monolithic_scheme),
41 _process_data(std::move(process_data))
50 std::make_unique<IntegrationPointWriter>(
56 std::make_unique<IntegrationPointWriter>(
63 template <
int DisplacementDim>
69 template <
int DisplacementDim>
72 const int process_id)
const
76 if (_use_monolithic_scheme || process_id == 1)
78 auto const& l = *_local_to_global_index_map;
79 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
80 &l.getGhostIndices(), &this->_sparsity_pattern};
84 auto const& l = *_local_to_global_index_map_with_base_nodes;
85 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
86 &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
89 template <
int DisplacementDim>
93 _mesh_subset_all_nodes =
94 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
97 _mesh_subset_base_nodes =
98 std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
102 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
103 *_mesh_subset_all_nodes};
104 _local_to_global_index_map_single_component =
105 std::make_unique<NumLib::LocalToGlobalIndexMap>(
106 std::move(all_mesh_subsets_single_component),
110 if (_use_monolithic_scheme)
113 std::vector<MeshLib::MeshSubset> all_mesh_subsets{
114 *_mesh_subset_base_nodes};
117 const int monolithic_process_id = 0;
118 std::generate_n(std::back_inserter(all_mesh_subsets),
119 getProcessVariables(monolithic_process_id)[1]
121 .getNumberOfGlobalComponents(),
122 [&]() {
return *_mesh_subset_all_nodes; });
124 std::vector<int>
const vec_n_components{1, DisplacementDim};
125 _local_to_global_index_map =
126 std::make_unique<NumLib::LocalToGlobalIndexMap>(
127 std::move(all_mesh_subsets), vec_n_components,
129 assert(_local_to_global_index_map);
134 const int process_id = 1;
135 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
136 std::generate_n(std::back_inserter(all_mesh_subsets),
137 getProcessVariables(process_id)[0]
139 .getNumberOfGlobalComponents(),
140 [&]() {
return *_mesh_subset_all_nodes; });
142 std::vector<int>
const vec_n_components{DisplacementDim};
143 _local_to_global_index_map =
144 std::make_unique<NumLib::LocalToGlobalIndexMap>(
145 std::move(all_mesh_subsets), vec_n_components,
150 std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
151 *_mesh_subset_base_nodes};
152 _local_to_global_index_map_with_base_nodes =
153 std::make_unique<NumLib::LocalToGlobalIndexMap>(
154 std::move(all_mesh_subsets_base_nodes),
159 *_local_to_global_index_map_with_base_nodes, _mesh);
161 assert(_local_to_global_index_map);
162 assert(_local_to_global_index_map_with_base_nodes);
166 template <
int DisplacementDim>
170 unsigned const integration_order)
172 const int mechanical_process_id = _use_monolithic_scheme ? 0 : 1;
173 const int deformation_variable_id = _use_monolithic_scheme ? 1 : 0;
178 getProcessVariables(mechanical_process_id)[deformation_variable_id]
180 .getShapeFunctionOrder(),
184 auto add_secondary_variable = [&](std::string
const&
name,
185 int const num_components,
186 auto get_ip_values_function)
188 _secondary_variables.addSecondaryVariable(
192 std::move(get_ip_values_function)));
195 add_secondary_variable(
"sigma",
197 DisplacementDim>::RowsAtCompileTime,
198 &LocalAssemblerIF::getIntPtSigma);
200 add_secondary_variable(
"epsilon",
202 DisplacementDim>::RowsAtCompileTime,
203 &LocalAssemblerIF::getIntPtEpsilon);
205 add_secondary_variable(
"velocity", DisplacementDim,
206 &LocalAssemblerIF::getIntPtDarcyVelocity);
213 add_secondary_variable);
215 _process_data.pressure_interpolated =
216 MeshLib::getOrCreateMeshProperty<double>(
220 _process_data.principal_stress_vector[0] =
221 MeshLib::getOrCreateMeshProperty<double>(
222 const_cast<MeshLib::Mesh&
>(mesh),
"principal_stress_vector_1",
225 _process_data.principal_stress_vector[1] =
226 MeshLib::getOrCreateMeshProperty<double>(
227 const_cast<MeshLib::Mesh&
>(mesh),
"principal_stress_vector_2",
230 _process_data.principal_stress_vector[2] =
231 MeshLib::getOrCreateMeshProperty<double>(
232 const_cast<MeshLib::Mesh&
>(mesh),
"principal_stress_vector_3",
235 _process_data.principal_stress_values =
236 MeshLib::getOrCreateMeshProperty<double>(
240 _process_data.permeability = MeshLib::getOrCreateMeshProperty<double>(
246 for (
auto const& ip_writer : _integration_point_writer)
249 auto const&
name = ip_writer->name();
254 auto const& mesh_property =
258 if (mesh_property.getMeshItemType() !=
267 if (ip_meta_data.n_components !=
268 mesh_property.getNumberOfGlobalComponents())
271 "Different number of components in meta data ({:d}) than in "
272 "the integration point field data for '{:s}': {:d}.",
273 ip_meta_data.n_components,
name,
274 mesh_property.getNumberOfGlobalComponents());
279 std::size_t position = 0;
280 for (
auto& local_asm : _local_assemblers)
282 std::size_t
const integration_points_read =
283 local_asm->setIPDataInitialConditions(
284 name, &mesh_property[position],
285 ip_meta_data.integration_order);
286 if (integration_points_read == 0)
289 "No integration points read in the integration point "
290 "initial conditions set function.");
292 position += integration_points_read * ip_meta_data.n_components;
299 *_local_to_global_index_map);
302 template <
int DisplacementDim>
305 if (_use_monolithic_scheme)
307 const int process_id_of_hydromechanics = 0;
308 initializeProcessBoundaryConditionsAndSourceTerms(
309 *_local_to_global_index_map, process_id_of_hydromechanics);
315 const int hydraulic_process_id = 0;
316 initializeProcessBoundaryConditionsAndSourceTerms(
317 *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
320 const int mechanical_process_id = 1;
321 initializeProcessBoundaryConditionsAndSourceTerms(
322 *_local_to_global_index_map, mechanical_process_id);
325 template <
int DisplacementDim>
327 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
328 std::vector<GlobalVector*>
const& xdot,
int const process_id,
331 DBUG(
"Assemble the equations for HydroMechanics");
337 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
338 dof_table = {std::ref(*_local_to_global_index_map)};
348 template <
int DisplacementDim>
351 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
352 std::vector<GlobalVector*>
const& xdot,
const double dxdot_dx,
353 const double dx_dx,
int const process_id,
GlobalMatrix& M,
356 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
359 if (_use_monolithic_scheme)
362 "Assemble the Jacobian of HydroMechanics for the monolithic "
364 dof_tables.emplace_back(*_local_to_global_index_map);
372 "Assemble the Jacobian equations of liquid fluid process in "
373 "HydroMechanics for the staggered scheme.");
378 "Assemble the Jacobian equations of mechanical process in "
379 "HydroMechanics for the staggered scheme.");
381 dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
382 dof_tables.emplace_back(*_local_to_global_index_map);
390 dxdot_dx, dx_dx, process_id, M, K, b, Jac);
392 auto copyRhs = [&](
int const variable_id,
auto& output_vector)
394 if (_use_monolithic_scheme)
398 std::negate<double>());
404 std::negate<double>());
407 if (_use_monolithic_scheme || process_id == 0)
409 copyRhs(0, *_hydraulic_flow);
411 if (_use_monolithic_scheme || process_id == 1)
413 copyRhs(1, *_nodal_forces);
417 template <
int DisplacementDim>
419 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
420 const int process_id)
422 DBUG(
"PreTimestep HydroMechanicsProcess.");
424 if (hasMechanicalProcess(process_id))
427 getProcessVariables(process_id)[0];
429 &LocalAssemblerIF::preTimestep, _local_assemblers,
431 *x[process_id], t, dt);
435 template <
int DisplacementDim>
437 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
438 const int process_id)
445 DBUG(
"PostTimestep HydroMechanicsProcess.");
446 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
447 auto const n_processes = x.size();
448 dof_tables.reserve(n_processes);
449 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
451 dof_tables.push_back(&getDOFTable(process_id));
456 &LocalAssemblerIF::postTimestep, _local_assemblers,
460 template <
int DisplacementDim>
463 double const dt,
const int process_id)
465 DBUG(
"PostNonLinearSolver HydroMechanicsProcess.");
469 &LocalAssemblerIF::postNonLinearSolver, _local_assemblers,
471 _use_monolithic_scheme, process_id);
474 template <
int DisplacementDim>
478 int const process_id)
480 if (process_id != _process_data.hydraulic_process_id)
485 DBUG(
"Set initial conditions of HydroMechanicsProcess.");
491 _use_monolithic_scheme, process_id);
494 template <
int DisplacementDim>
496 double const t,
double const dt, std::vector<GlobalVector*>
const& x,
504 DBUG(
"Compute the secondary variables for HydroMechanicsProcess.");
505 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
506 auto const n_processes = x.size();
507 dof_tables.reserve(n_processes);
508 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
510 dof_tables.push_back(&getDOFTable(process_id));
515 &LocalAssemblerIF::computeSecondaryVariable, _local_assemblers,
519 template <
int DisplacementDim>
520 std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
523 const bool manage_storage =
false;
524 return std::make_tuple(_local_to_global_index_map_single_component.get(),
528 template <
int DisplacementDim>
532 if (hasMechanicalProcess(process_id))
534 return *_local_to_global_index_map;
538 return *_local_to_global_index_map_with_base_nodes;
void DBUG(char const *fmt, Args const &... args)
Global vector based on Eigen vector.
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.
Properties & getProperties()
bool existsPropertyVector(std::string const &name) const
void postNonLinearSolverConcreteProcess(GlobalVector const &x, GlobalVector const &xdot, const double t, double const dt, int const process_id) override
std::vector< std::unique_ptr< LocalAssemblerIF > > _local_assemblers
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 preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double dt, int const process_id) override
void constructDofTable() override
void initializeBoundaryConditions() override
bool isLinear() const override
MeshLib::PropertyVector< double > * _nodal_forces
void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &x, double const t, int const process_id) override
MeshLib::PropertyVector< double > * _hydraulic_flow
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
void assembleWithJacobianConcreteProcess(const double t, double const, 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 initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
void assembleConcreteProcess(const double t, double const, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, const int process_id) override
std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData() const override
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::unique_ptr< IntegrationPointWriter > > _integration_point_writer
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)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
std::vector< Node * > getBaseNodes(std::vector< Element * > const &elements)
@ BY_LOCATION
Ordering data by spatial location.
GlobalSparsityPattern computeSparsityPattern(LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
Computes a sparsity pattern for the given inputs.
void transformVariableFromGlobalVector(GlobalVector const &input_vector, int const variable_id, NumLib::LocalToGlobalIndexMap const &local_to_global_index_map, MeshLib::PropertyVector< double > &output_vector, Functor mapFunction)
void createLocalAssemblers(const unsigned, 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)
IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Mesh const &mesh, std::string const &name)
std::pair< std::vector< GlobalVector * >, std::vector< GlobalVector * > > setInitialConditions(double const t0, std::vector< std::unique_ptr< ProcessData >> const &per_process_data)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)
virtual std::vector< double > getEpsilon() const =0
virtual std::vector< double > getSigma() const =0