OGS
HydroMechanicsProcess.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14
15#include "HydroMechanicsFEM.h"
21#include "ProcessLib/Process.h"
24
25namespace ProcessLib
26{
27namespace HydroMechanics
28{
29template <int DisplacementDim>
31 std::string name,
32 MeshLib::Mesh& mesh,
33 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
34 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
35 unsigned const integration_order,
36 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
37 process_variables,
39 SecondaryVariableCollection&& secondary_variables,
40 bool const use_monolithic_scheme)
41 : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
42 integration_order, std::move(process_variables),
43 std::move(secondary_variables), use_monolithic_scheme),
44 _process_data(std::move(process_data))
45{
46 _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
47 mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
48
49 _hydraulic_flow = MeshLib::getOrCreateMeshProperty<double>(
50 mesh, "MassFlowRate", MeshLib::MeshItemType::Node, 1);
51
52 _integration_point_writer.emplace_back(
53 std::make_unique<MeshLib::IntegrationPointWriter>(
54 "sigma_ip",
55 static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
56 integration_order, _local_assemblers, &LocalAssemblerIF::getSigma));
57
58 _integration_point_writer.emplace_back(
59 std::make_unique<MeshLib::IntegrationPointWriter>(
60 "epsilon_ip",
61 static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
62 integration_order, _local_assemblers,
64
66 {
67 _integration_point_writer.emplace_back(
68 std::make_unique<MeshLib::IntegrationPointWriter>(
69 "strain_rate_variable_ip", 1, integration_order,
71 }
72}
73
74template <int DisplacementDim>
76{
77 return false;
78}
79
80template <int DisplacementDim>
83 const int process_id) const
84{
85 // For the monolithic scheme or the M process (deformation) in the staggered
86 // scheme.
87 if (process_id == _process_data.mechanics_related_process_id)
88 {
89 auto const& l = *_local_to_global_index_map;
90 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
91 &l.getGhostIndices(), &this->_sparsity_pattern};
92 }
93
94 // For staggered scheme and H process (pressure).
95 auto const& l = *_local_to_global_index_map_with_base_nodes;
96 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
97 &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
98}
99
100template <int DisplacementDim>
102{
103 // Create single component dof in every of the mesh's nodes.
104 _mesh_subset_all_nodes = std::make_unique<MeshLib::MeshSubset>(
105 _mesh, _mesh.getNodes(), _process_data.use_taylor_hood_elements);
106
107 // Create single component dof in the mesh's base nodes.
108 _base_nodes = MeshLib::getBaseNodes(_mesh.getElements());
109 _mesh_subset_base_nodes = std::make_unique<MeshLib::MeshSubset>(
110 _mesh, _base_nodes, _process_data.use_taylor_hood_elements);
111
112 // TODO move the two data members somewhere else.
113 // for extrapolation of secondary variables of stress or strain
114 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
115 *_mesh_subset_all_nodes};
116 _local_to_global_index_map_single_component =
117 std::make_unique<NumLib::LocalToGlobalIndexMap>(
118 std::move(all_mesh_subsets_single_component),
119 // by location order is needed for output
121
122 if (_process_data.isMonolithicSchemeUsed())
123 {
124 // For pressure, which is the first
125 std::vector<MeshLib::MeshSubset> all_mesh_subsets{
126 *_mesh_subset_base_nodes};
127
128 // For displacement.
129 const int monolithic_process_id = 0;
130 std::generate_n(std::back_inserter(all_mesh_subsets),
131 getProcessVariables(monolithic_process_id)[1]
132 .get()
133 .getNumberOfGlobalComponents(),
134 [&]() { return *_mesh_subset_all_nodes; });
135
136 std::vector<int> const vec_n_components{1, DisplacementDim};
137 _local_to_global_index_map =
138 std::make_unique<NumLib::LocalToGlobalIndexMap>(
139 std::move(all_mesh_subsets), vec_n_components,
141 assert(_local_to_global_index_map);
142 }
143 else
144 {
145 // For displacement equation.
146 const int process_id = 1;
147 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
148 std::generate_n(std::back_inserter(all_mesh_subsets),
149 getProcessVariables(process_id)[0]
150 .get()
151 .getNumberOfGlobalComponents(),
152 [&]() { return *_mesh_subset_all_nodes; });
153
154 std::vector<int> const vec_n_components{DisplacementDim};
155 _local_to_global_index_map =
156 std::make_unique<NumLib::LocalToGlobalIndexMap>(
157 std::move(all_mesh_subsets), vec_n_components,
159
160 // For pressure equation.
161 // Collect the mesh subsets with base nodes in a vector.
162 std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
163 *_mesh_subset_base_nodes};
164 _local_to_global_index_map_with_base_nodes =
165 std::make_unique<NumLib::LocalToGlobalIndexMap>(
166 std::move(all_mesh_subsets_base_nodes),
167 // by location order is needed for output
169
170 _sparsity_pattern_with_linear_element = NumLib::computeSparsityPattern(
171 *_local_to_global_index_map_with_base_nodes, _mesh);
172
173 assert(_local_to_global_index_map);
174 assert(_local_to_global_index_map_with_base_nodes);
175 }
176}
177
178template <int DisplacementDim>
180 NumLib::LocalToGlobalIndexMap const& dof_table,
181 MeshLib::Mesh const& mesh,
182 unsigned const integration_order)
183{
186 mesh.getElements(), dof_table, _local_assemblers,
187 NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
188 _process_data);
189
190 auto add_secondary_variable = [&](std::string const& name,
191 int const num_components,
192 auto get_ip_values_function)
193 {
194 _secondary_variables.addSecondaryVariable(
195 name,
196 makeExtrapolator(num_components, getExtrapolator(),
197 _local_assemblers,
198 std::move(get_ip_values_function)));
199 };
200
201 add_secondary_variable("sigma",
203 DisplacementDim>::RowsAtCompileTime,
204 &LocalAssemblerIF::getIntPtSigma);
205
206 add_secondary_variable("epsilon",
208 DisplacementDim>::RowsAtCompileTime,
209 &LocalAssemblerIF::getIntPtEpsilon);
210
211 add_secondary_variable("velocity", DisplacementDim,
212 &LocalAssemblerIF::getIntPtDarcyVelocity);
213
214 //
215 // enable output of internal variables defined by material models
216 //
218 LocalAssemblerIF>(_process_data.solid_materials,
219 add_secondary_variable);
220
221 _process_data.pressure_interpolated =
222 MeshLib::getOrCreateMeshProperty<double>(
223 const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
225
226 _process_data.principal_stress_vector[0] =
227 MeshLib::getOrCreateMeshProperty<double>(
228 const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_1",
230
231 _process_data.principal_stress_vector[1] =
232 MeshLib::getOrCreateMeshProperty<double>(
233 const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_2",
235
236 _process_data.principal_stress_vector[2] =
237 MeshLib::getOrCreateMeshProperty<double>(
238 const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_3",
240
241 _process_data.principal_stress_values =
242 MeshLib::getOrCreateMeshProperty<double>(
243 const_cast<MeshLib::Mesh&>(mesh), "principal_stress_values",
245
246 _process_data.permeability = MeshLib::getOrCreateMeshProperty<double>(
247 const_cast<MeshLib::Mesh&>(mesh), "permeability",
250
251 setIPDataInitialConditions(_integration_point_writer, mesh.getProperties(),
252 _local_assemblers);
253
254 // Initialize local assemblers after all variables have been set.
255 GlobalExecutor::executeMemberOnDereferenced(&LocalAssemblerIF::initialize,
256 _local_assemblers,
257 *_local_to_global_index_map);
258}
259
260template <int DisplacementDim>
262 std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
263{
264 if (_process_data.isMonolithicSchemeUsed())
265 {
266 const int process_id_of_hydromechanics = 0;
267 initializeProcessBoundaryConditionsAndSourceTerms(
268 *_local_to_global_index_map, process_id_of_hydromechanics, media);
269 return;
270 }
271
272 // Staggered scheme:
273 // for the equations of pressure
274 const int hydraulic_process_id = 0;
275 initializeProcessBoundaryConditionsAndSourceTerms(
276 *_local_to_global_index_map_with_base_nodes, hydraulic_process_id,
277 media);
278
279 // for the equations of deformation.
280 const int mechanical_process_id = 1;
281 initializeProcessBoundaryConditionsAndSourceTerms(
282 *_local_to_global_index_map, mechanical_process_id, media);
283}
284
285template <int DisplacementDim>
287 const double t, double const dt, std::vector<GlobalVector*> const& x,
288 std::vector<GlobalVector*> const& x_prev, int const process_id,
290{
291 DBUG("Assemble the equations for HydroMechanics");
292
293 // Note: This assembly function is for the Picard nonlinear solver. Since
294 // only the Newton-Raphson method is employed to simulate coupled HM
295 // processes in this class, this function is actually not used so far.
296
297 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
298 _local_to_global_index_map.get()};
299
301 _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
302 getActiveElementIDs(), dof_table, t, dt, x, x_prev, process_id, M, K,
303 b);
304}
305
306template <int DisplacementDim>
309 const double t, double const dt, std::vector<GlobalVector*> const& x,
310 std::vector<GlobalVector*> const& x_prev, int const process_id,
312{
313 // For the monolithic scheme
314 bool const use_monolithic_scheme = _process_data.isMonolithicSchemeUsed();
315 if (use_monolithic_scheme)
316 {
317 DBUG(
318 "Assemble the Jacobian of HydroMechanics for the monolithic "
319 "scheme.");
320 }
321 else
322 {
323 // For the staggered scheme
324 if (process_id == _process_data.hydraulic_process_id)
325 {
326 DBUG(
327 "Assemble the Jacobian equations of liquid fluid process in "
328 "HydroMechanics for the staggered scheme.");
329 }
330 else
331 {
332 DBUG(
333 "Assemble the Jacobian equations of mechanical process in "
334 "HydroMechanics for the staggered scheme.");
335 }
336 }
337
338 auto const dof_tables = getDOFTables(x.size());
341 _local_assemblers, getActiveElementIDs(), dof_tables, t, dt, x, x_prev,
342 process_id, M, K, b, Jac);
343
344 auto copyRhs = [&](int const variable_id, auto& output_vector)
345 {
346 if (use_monolithic_scheme)
347 {
348 transformVariableFromGlobalVector(b, variable_id, *dof_tables[0],
349 output_vector,
350 std::negate<double>());
351 }
352 else
353 {
354 transformVariableFromGlobalVector(b, 0, *dof_tables[process_id],
355 output_vector,
356 std::negate<double>());
357 }
358 };
359 if (process_id == _process_data.hydraulic_process_id)
360 {
361 copyRhs(0, *_hydraulic_flow);
362 }
363 if (process_id == _process_data.mechanics_related_process_id)
364 {
365 copyRhs(1, *_nodal_forces);
366 }
367}
368
369template <int DisplacementDim>
371 std::vector<GlobalVector*> const& x, double const t, double const dt,
372 const int process_id)
373{
374 DBUG("PreTimestep HydroMechanicsProcess.");
375
376 if (hasMechanicalProcess(process_id))
377 {
379 &LocalAssemblerIF::preTimestep, _local_assemblers,
380 getActiveElementIDs(), *_local_to_global_index_map, *x[process_id],
381 t, dt);
382 }
383}
384
385template <int DisplacementDim>
387 std::vector<GlobalVector*> const& x,
388 std::vector<GlobalVector*> const& x_prev, double const t, double const dt,
389 const int process_id)
390{
391 if (process_id != _process_data.hydraulic_process_id)
392 {
393 return;
394 }
395
396 DBUG("PostTimestep HydroMechanicsProcess.");
397
399 &LocalAssemblerIF::postTimestep, _local_assemblers,
400 getActiveElementIDs(), getDOFTables(x.size()), x, x_prev, t, dt,
401 process_id);
402}
403
404template <int DisplacementDim>
406 std::vector<GlobalVector*> const& x,
407 std::vector<GlobalVector*> const& x_prev, const double t, double const dt,
408 const int process_id)
409{
410 DBUG("PostNonLinearSolver HydroMechanicsProcess.");
411
412 // Calculate strain, stress or other internal variables of mechanics.
414 &LocalAssemblerIF::postNonLinearSolver, _local_assemblers,
415 getActiveElementIDs(), getDOFTables(x.size()), x, x_prev, t, dt,
416 process_id);
417}
418
419template <int DisplacementDim>
421 setInitialConditionsConcreteProcess(std::vector<GlobalVector*>& x,
422 double const t,
423 int const process_id)
424{
425 // So far, this function only sets the initial stress using the input data.
426 if (process_id != _process_data.mechanics_related_process_id)
427 {
428 return;
429 }
430
431 DBUG("Set initial conditions of HydroMechanicsProcess.");
432
434 &LocalAssemblerIF::setInitialConditions, _local_assemblers,
435 getActiveElementIDs(), getDOFTables(x.size()), x, t, process_id);
436}
437
438template <int DisplacementDim>
440 double const t, double const dt, std::vector<GlobalVector*> const& x,
441 GlobalVector const& x_prev, const int process_id)
442{
443 if (process_id != _process_data.hydraulic_process_id)
444 {
445 return;
446 }
447
448 DBUG("Compute the secondary variables for HydroMechanicsProcess.");
449
451 &LocalAssemblerIF::computeSecondaryVariable, _local_assemblers,
452 getActiveElementIDs(), getDOFTables(x.size()), t, dt, x, x_prev,
453 process_id);
454}
455
456template <int DisplacementDim>
457std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
459{
460 const bool manage_storage = false;
461 return std::make_tuple(_local_to_global_index_map_single_component.get(),
462 manage_storage);
463}
464
465template <int DisplacementDim>
468{
469 if (hasMechanicalProcess(process_id))
470 {
471 return *_local_to_global_index_map;
472 }
473
474 // For the equation of pressure
475 return *_local_to_global_index_map_with_base_nodes;
476}
477
478template class HydroMechanicsProcess<2>;
479template class HydroMechanicsProcess<3>;
480
481} // namespace HydroMechanics
482} // namespace ProcessLib
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Global vector based on Eigen vector.
Definition EigenVector.h:25
bool isAxiallySymmetric() const
Definition Mesh.h:137
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:109
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition Mesh.h:88
Properties & getProperties()
Definition Mesh.h:134
void assembleConcreteProcess(const double t, double const, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) 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 &parameters, 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 computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, const int process_id) override
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
void initializeBoundaryConditions(std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media) override
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
void assembleWithJacobianConcreteProcess(const double t, double const, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
void postNonLinearSolverConcreteProcess(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, const double t, double const dt, int const process_id) override
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 initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData() const override
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, const double t, const double dt, int const process_id) override
std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > _integration_point_writer
Definition Process.h:383
const bool _use_monolithic_scheme
Definition Process.h:372
Handles configuration of several secondary variables from the project file.
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, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac)
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)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
std::vector< Node * > getBaseNodes(std::vector< Element * > const &elements)
Definition Utils.h:26
@ 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 solidMaterialInternalToSecondaryVariables(std::map< int, std::unique_ptr< SolidMaterial > > const &solid_materials, AddSecondaryVariableCallback const &add_secondary_variable)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
void setIPDataInitialConditions(std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > const &_integration_point_writer, MeshLib::Properties const &mesh_properties, LocalAssemblersVector &local_assemblers)
void createLocalAssemblersHM(std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface > > &local_assemblers, ProviderOrOrder const &provider_or_order, ExtraCtorArgs &&... extra_ctor_args)
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 > getSigma() const =0
virtual std::vector< double > getStrainRateVariable() const =0
virtual std::vector< double > getEpsilon() const =0