Loading [MathJax]/extensions/tex2jax.js
OGS
RichardsMechanicsProcess.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14
25
26namespace ProcessLib
27{
28namespace RichardsMechanics
29{
30template <int DisplacementDim>
32 std::string name,
33 MeshLib::Mesh& mesh,
34 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
35 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
36 unsigned const integration_order,
37 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
38 process_variables,
40 SecondaryVariableCollection&& secondary_variables,
41 bool const use_monolithic_scheme)
42 : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
43 integration_order, std::move(process_variables),
44 std::move(secondary_variables), use_monolithic_scheme),
47 process_data_(std::move(process_data))
48{
50 mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
51
53 mesh, "MassFlowRate", MeshLib::MeshItemType::Node, 1);
54
56 DisplacementDim>(LocalAssemblerIF::getReflectionDataForOutput(),
57 _integration_point_writer, integration_order,
58 local_assemblers_);
59}
60
61template <int DisplacementDim>
63{
64 return false;
65}
66
67template <int DisplacementDim>
70 const int process_id) const
71{
72 // For the monolithic scheme or the M process (deformation) in the staggered
73 // scheme.
74 if (_use_monolithic_scheme || process_id == 1)
75 {
76 auto const& l = *_local_to_global_index_map;
77 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
78 &l.getGhostIndices(), &this->_sparsity_pattern};
79 }
80
81 // For staggered scheme and H process (pressure).
82 auto const& l = *local_to_global_index_map_with_base_nodes_;
83 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
84 &l.getGhostIndices(), &sparsity_pattern_with_linear_element_};
85}
86
87template <int DisplacementDim>
89{
90 // Create single component dof in every of the mesh's nodes.
91 _mesh_subset_all_nodes =
92 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
93 // Create single component dof in the mesh's base nodes.
94 base_nodes_ = MeshLib::getBaseNodes(_mesh.getElements());
95 mesh_subset_base_nodes_ =
96 std::make_unique<MeshLib::MeshSubset>(_mesh, base_nodes_);
97
98 // TODO move the two data members somewhere else.
99 // for extrapolation of secondary variables of stress or strain
100 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
101 *_mesh_subset_all_nodes};
102 local_to_global_index_map_single_component_ =
103 std::make_unique<NumLib::LocalToGlobalIndexMap>(
104 std::move(all_mesh_subsets_single_component),
105 // by location order is needed for output
107
108 if (_use_monolithic_scheme)
109 {
110 // For pressure, which is the first
111 std::vector<MeshLib::MeshSubset> all_mesh_subsets{
112 *mesh_subset_base_nodes_};
113
114 // For displacement.
115 const int monolithic_process_id = 0;
116 std::generate_n(std::back_inserter(all_mesh_subsets),
117 getProcessVariables(monolithic_process_id)[1]
118 .get()
119 .getNumberOfGlobalComponents(),
120 [&]() { return *_mesh_subset_all_nodes; });
121
122 std::vector<int> const vec_n_components{1, DisplacementDim};
123 _local_to_global_index_map =
124 std::make_unique<NumLib::LocalToGlobalIndexMap>(
125 std::move(all_mesh_subsets), vec_n_components,
127 assert(_local_to_global_index_map);
128 }
129 else
130 {
131 // For displacement equation.
132 const int process_id = 1;
133 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
134 std::generate_n(std::back_inserter(all_mesh_subsets),
135 getProcessVariables(process_id)[0]
136 .get()
137 .getNumberOfGlobalComponents(),
138 [&]() { return *_mesh_subset_all_nodes; });
139
140 std::vector<int> const vec_n_components{DisplacementDim};
141 _local_to_global_index_map =
142 std::make_unique<NumLib::LocalToGlobalIndexMap>(
143 std::move(all_mesh_subsets), vec_n_components,
145
146 // For pressure equation.
147 // Collect the mesh subsets with base nodes in a vector.
148 std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
149 *mesh_subset_base_nodes_};
150 local_to_global_index_map_with_base_nodes_ =
151 std::make_unique<NumLib::LocalToGlobalIndexMap>(
152 std::move(all_mesh_subsets_base_nodes),
153 // by location order is needed for output
155
156 sparsity_pattern_with_linear_element_ = NumLib::computeSparsityPattern(
157 *local_to_global_index_map_with_base_nodes_, _mesh);
158
159 assert(_local_to_global_index_map);
160 assert(local_to_global_index_map_with_base_nodes_);
161 }
162}
163
164template <int DisplacementDim>
166 NumLib::LocalToGlobalIndexMap const& dof_table,
167 MeshLib::Mesh const& mesh,
168 unsigned const integration_order)
169{
172 mesh.getElements(), dof_table, local_assemblers_,
173 NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
174 process_data_);
175
177 LocalAssemblerIF::getReflectionDataForOutput(), _secondary_variables,
178 getExtrapolator(), local_assemblers_);
179
180 auto add_secondary_variable = [&](std::string const& name,
181 int const num_components,
182 auto get_ip_values_function)
183 {
184 _secondary_variables.addSecondaryVariable(
185 name,
186 makeExtrapolator(num_components, getExtrapolator(),
187 local_assemblers_,
188 std::move(get_ip_values_function)));
189 };
190
191 //
192 // enable output of internal variables defined by material models
193 //
195 LocalAssemblerIF>(process_data_.solid_materials,
196 add_secondary_variable);
197
200 process_data_.solid_materials, local_assemblers_,
201 _integration_point_writer, integration_order);
202
203 process_data_.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
204 const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
206
207 process_data_.element_porosity = MeshLib::getOrCreateMeshProperty<double>(
208 const_cast<MeshLib::Mesh&>(mesh), "porosity_avg",
210
211 process_data_.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
212 const_cast<MeshLib::Mesh&>(mesh), "stress_avg",
215 DisplacementDim>::RowsAtCompileTime);
216
217 process_data_.pressure_interpolated =
219 const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
221
222 setIPDataInitialConditions(_integration_point_writer, mesh.getProperties(),
223 local_assemblers_);
224
225 // Initialize local assemblers after all variables have been set.
226 GlobalExecutor::executeMemberOnDereferenced(&LocalAssemblerIF::initialize,
227 local_assemblers_,
228 *_local_to_global_index_map);
229}
230
231template <int DisplacementDim>
233 std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
234{
235 if (_use_monolithic_scheme)
236 {
237 const int monolithic_process_id = 0;
238 initializeProcessBoundaryConditionsAndSourceTerms(
239 *_local_to_global_index_map, monolithic_process_id, media);
240 return;
241 }
242
243 // Staggered scheme:
244 // for the equations of pressure
245 const int hydraulic_process_id = 0;
246 initializeProcessBoundaryConditionsAndSourceTerms(
247 *local_to_global_index_map_with_base_nodes_, hydraulic_process_id,
248 media);
249
250 // for the equations of deformation.
251 const int mechanical_process_id = 1;
252 initializeProcessBoundaryConditionsAndSourceTerms(
253 *_local_to_global_index_map, mechanical_process_id, media);
254}
255
256template <int DisplacementDim>
258 setInitialConditionsConcreteProcess(std::vector<GlobalVector*>& x,
259 double const t,
260 int const process_id)
261{
262 if (process_id != 0)
263 {
264 return;
265 }
266
267 DBUG("SetInitialConditions RichardsMechanicsProcess.");
268
270 &LocalAssemblerIF::setInitialConditions, local_assemblers_,
271 getActiveElementIDs(), getDOFTables(x.size()), x, t, process_id);
272}
273
274template <int DisplacementDim>
276 const double t, double const dt, std::vector<GlobalVector*> const& x,
277 std::vector<GlobalVector*> const& x_prev, int const process_id,
279{
280 DBUG("Assemble the equations for RichardsMechanics");
281
283 t, dt, x, x_prev, process_id, M, K, b);
284}
285
286template <int DisplacementDim>
289 const double t, double const dt, std::vector<GlobalVector*> const& x,
290 std::vector<GlobalVector*> const& x_prev, int const process_id,
291 GlobalVector& b, GlobalMatrix& Jac)
292{
293 // For the monolithic scheme
294 if (_use_monolithic_scheme)
295 {
296 DBUG(
297 "Assemble the Jacobian of RichardsMechanics for the monolithic"
298 " scheme.");
299 }
300 else
301 {
302 // For the staggered scheme
303 if (process_id == 0)
304 {
305 DBUG(
306 "Assemble the Jacobian equations of liquid fluid process in "
307 "RichardsMechanics for the staggered scheme.");
308 }
309 else
310 {
311 DBUG(
312 "Assemble the Jacobian equations of mechanical process in "
313 "RichardsMechanics for the staggered scheme.");
314 }
315 }
316
317 auto const dof_tables = getDOFTables(x.size());
319 assembleWithJacobian(t, dt, x, x_prev, process_id, b, Jac);
320
321 auto copyRhs = [&](int const variable_id, auto& output_vector)
322 {
323 if (_use_monolithic_scheme)
324 {
325 transformVariableFromGlobalVector(b, variable_id, *dof_tables[0],
326 output_vector,
327 std::negate<double>());
328 }
329 else
330 {
331 transformVariableFromGlobalVector(b, 0, *dof_tables[process_id],
332 output_vector,
333 std::negate<double>());
334 }
335 };
336 if (_use_monolithic_scheme || process_id == 0)
337 {
338 copyRhs(0, *hydraulic_flow_);
339 }
340 if (_use_monolithic_scheme || process_id == 1)
341 {
342 copyRhs(1, *nodal_forces_);
343 }
344}
345
346template <int DisplacementDim>
348 std::vector<GlobalVector*> const& x, double const t, double const dt,
349 const int process_id)
350{
351 DBUG("PreTimestep RichardsMechanicsProcess.");
352
354 &LocalAssemblerIF::preTimestep, local_assemblers_,
355 getActiveElementIDs(), *_local_to_global_index_map, *x[process_id], t,
356 dt);
357
359 RichardsMechanicsProcess<DisplacementDim>>::updateActiveElements();
360}
361
362template <int DisplacementDim>
364 std::vector<GlobalVector*> const& x,
365 std::vector<GlobalVector*> const& x_prev, double const t, double const dt,
366 const int process_id)
367{
368 if (hasMechanicalProcess(process_id))
369 {
370 DBUG("PostTimestep RichardsMechanicsProcess.");
371
373 &LocalAssemblerIF::postTimestep, local_assemblers_,
374 getActiveElementIDs(), getDOFTables(x.size()), x, x_prev, t, dt,
375 process_id);
376 }
377}
378
379template <int DisplacementDim>
380std::vector<std::vector<std::string>>
382 std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes)
383{
384 INFO("RichardsMechanics process initializeSubmeshOutput().");
385 std::vector<std::vector<std::string>> residuum_names{
386 {"MassFlowRate", "NodalForces"}};
387
389 initializeAssemblyOnSubmeshes(meshes, residuum_names);
390
391 return residuum_names;
392}
393
394template <int DisplacementDim>
396 computeSecondaryVariableConcrete(const double t, const double dt,
397 std::vector<GlobalVector*> const& x,
398 GlobalVector const& x_prev,
399 int const process_id)
400{
401 if (process_id != 0)
402 {
403 return;
404 }
405
406 DBUG("Compute the secondary variables for RichardsMechanicsProcess.");
407
409 &LocalAssemblerIF::computeSecondaryVariable, local_assemblers_,
410 getActiveElementIDs(), getDOFTables(x.size()), t, dt, x, x_prev,
411 process_id);
412}
413
414template <int DisplacementDim>
415std::tuple<NumLib::LocalToGlobalIndexMap*, bool> RichardsMechanicsProcess<
416 DisplacementDim>::getDOFTableForExtrapolatorData() const
417{
418 const bool manage_storage = false;
419 return std::make_tuple(local_to_global_index_map_single_component_.get(),
420 manage_storage);
421}
422
423template <int DisplacementDim>
426 const int process_id) const
427{
428 if (hasMechanicalProcess(process_id))
429 {
430 return *_local_to_global_index_map;
431 }
432
433 // For the equation of pressure
434 return *local_to_global_index_map_with_base_nodes_;
435}
436
437template class RichardsMechanicsProcess<2>;
438template class RichardsMechanicsProcess<3>;
439
440} // namespace RichardsMechanics
441} // namespace ProcessLib
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
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:139
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
Properties & getProperties()
Definition Mesh.h:136
std::unique_ptr< ProcessLib::AbstractJacobianAssembler > _jacobian_assembler
Definition Process.h:376
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 preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
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().
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id) override
RichardsMechanicsProcess(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, RichardsMechanicsProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &x, double const t, int const process_id) override
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
void initializeBoundaryConditions(std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media) override
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, const int process_id) override
std::vector< std::vector< std::string > > initializeAssemblyOnSubmeshes(std::vector< std::reference_wrapper< MeshLib::Mesh > > const &meshes) override
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
Handles configuration of several secondary variables from the project file.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
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)
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 solidMaterialInternalVariablesToIntegrationPointWriter(std::map< int, std::shared_ptr< SolidMaterial > > const &solid_materials, std::vector< std::unique_ptr< LocalAssemblerInterface > > const &local_assemblers, std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > &integration_point_writer, int const integration_order)
void solidMaterialInternalToSecondaryVariables(std::map< int, std::shared_ptr< SolidMaterial > > const &solid_materials, AddSecondaryVariableCallback const &add_secondary_variable)
void addReflectedSecondaryVariables(ReflData const &reflection_data, SecondaryVariableCollection &secondary_variables, NumLib::Extrapolator &extrapolator, std::vector< std::unique_ptr< LocAsmIF > > const &local_assemblers)
void addReflectedIntegrationPointWriters(ReflData const &reflection_data, std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > &integration_point_writers, unsigned const integration_order, std::vector< std::unique_ptr< LocAsmIF > > const &local_assemblers)
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 executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)