24 namespace RichardsMechanics
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))
53 std::make_unique<IntegrationPointWriter>(
59 std::make_unique<IntegrationPointWriter>(
60 "saturation_ip", 1 , integration_order,
64 std::make_unique<IntegrationPointWriter>(
65 "porosity_ip", 1 , integration_order,
69 std::make_unique<IntegrationPointWriter>(
70 "transport_porosity_ip", 1 , integration_order,
74 std::make_unique<IntegrationPointWriter>(
81 std::make_unique<IntegrationPointWriter>(
88 template <
int DisplacementDim>
94 template <
int DisplacementDim>
97 const int process_id)
const
101 if (_use_monolithic_scheme || process_id == 1)
103 auto const& l = *_local_to_global_index_map;
104 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
105 &l.getGhostIndices(), &this->_sparsity_pattern};
109 auto const& l = *_local_to_global_index_map_with_base_nodes;
110 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
111 &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
114 template <
int DisplacementDim>
118 _mesh_subset_all_nodes =
119 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
122 _mesh_subset_base_nodes =
123 std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
127 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
128 *_mesh_subset_all_nodes};
129 _local_to_global_index_map_single_component =
130 std::make_unique<NumLib::LocalToGlobalIndexMap>(
131 std::move(all_mesh_subsets_single_component),
135 if (_use_monolithic_scheme)
138 std::vector<MeshLib::MeshSubset> all_mesh_subsets{
139 *_mesh_subset_base_nodes};
142 const int monolithic_process_id = 0;
143 std::generate_n(std::back_inserter(all_mesh_subsets),
144 getProcessVariables(monolithic_process_id)[1]
146 .getNumberOfGlobalComponents(),
147 [&]() {
return *_mesh_subset_all_nodes; });
149 std::vector<int>
const vec_n_components{1, DisplacementDim};
150 _local_to_global_index_map =
151 std::make_unique<NumLib::LocalToGlobalIndexMap>(
152 std::move(all_mesh_subsets), vec_n_components,
154 assert(_local_to_global_index_map);
159 const int process_id = 1;
160 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
161 std::generate_n(std::back_inserter(all_mesh_subsets),
162 getProcessVariables(process_id)[0]
164 .getNumberOfGlobalComponents(),
165 [&]() {
return *_mesh_subset_all_nodes; });
167 std::vector<int>
const vec_n_components{DisplacementDim};
168 _local_to_global_index_map =
169 std::make_unique<NumLib::LocalToGlobalIndexMap>(
170 std::move(all_mesh_subsets), vec_n_components,
175 std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
176 *_mesh_subset_base_nodes};
177 _local_to_global_index_map_with_base_nodes =
178 std::make_unique<NumLib::LocalToGlobalIndexMap>(
179 std::move(all_mesh_subsets_base_nodes),
184 *_local_to_global_index_map_with_base_nodes, _mesh);
186 assert(_local_to_global_index_map);
187 assert(_local_to_global_index_map_with_base_nodes);
191 template <
int DisplacementDim>
195 unsigned const integration_order)
197 using nlohmann::json;
199 const int mechanical_process_id = _use_monolithic_scheme ? 0 : 1;
200 const int deformation_variable_id = _use_monolithic_scheme ? 1 : 0;
205 getProcessVariables(mechanical_process_id)[deformation_variable_id]
207 .getShapeFunctionOrder(),
211 auto add_secondary_variable = [&](std::string
const&
name,
212 int const num_components,
213 auto get_ip_values_function)
215 _secondary_variables.addSecondaryVariable(
219 std::move(get_ip_values_function)));
222 add_secondary_variable(
"sigma",
224 DisplacementDim>::RowsAtCompileTime,
225 &LocalAssemblerIF::getIntPtSigma);
227 add_secondary_variable(
"swelling_stress",
229 DisplacementDim>::RowsAtCompileTime,
230 &LocalAssemblerIF::getIntPtSwellingStress);
232 add_secondary_variable(
"epsilon",
234 DisplacementDim>::RowsAtCompileTime,
235 &LocalAssemblerIF::getIntPtEpsilon);
237 add_secondary_variable(
"velocity", DisplacementDim,
238 &LocalAssemblerIF::getIntPtDarcyVelocity);
240 add_secondary_variable(
"saturation", 1,
241 &LocalAssemblerIF::getIntPtSaturation);
243 add_secondary_variable(
"micro_saturation", 1,
244 &LocalAssemblerIF::getIntPtMicroSaturation);
246 add_secondary_variable(
"micro_pressure", 1,
247 &LocalAssemblerIF::getIntPtMicroPressure);
249 add_secondary_variable(
"porosity", 1, &LocalAssemblerIF::getIntPtPorosity);
251 add_secondary_variable(
"transport_porosity", 1,
252 &LocalAssemblerIF::getIntPtTransportPorosity);
254 add_secondary_variable(
"dry_density_solid", 1,
255 &LocalAssemblerIF::getIntPtDryDensitySolid);
262 add_secondary_variable);
267 _process_data.solid_materials, _local_assemblers,
268 _integration_point_writer, integration_order);
270 _process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
274 _process_data.element_porosity = MeshLib::getOrCreateMeshProperty<double>(
278 _process_data.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
282 DisplacementDim>::RowsAtCompileTime);
284 _process_data.pressure_interpolated =
285 MeshLib::getOrCreateMeshProperty<double>(
290 for (
auto const& ip_writer : _integration_point_writer)
293 auto const&
name = ip_writer->name();
298 auto const& mesh_property =
302 if (mesh_property.getMeshItemType() !=
311 if (ip_meta_data.n_components !=
312 mesh_property.getNumberOfGlobalComponents())
315 "Different number of components in meta data ({:d}) than in "
316 "the integration point field data for '{:s}': {:d}.",
317 ip_meta_data.n_components,
name,
318 mesh_property.getNumberOfGlobalComponents());
323 std::size_t position = 0;
324 for (
auto& local_asm : _local_assemblers)
326 std::size_t
const integration_points_read =
327 local_asm->setIPDataInitialConditions(
328 name, &mesh_property[position],
329 ip_meta_data.integration_order);
330 if (integration_points_read == 0)
333 "No integration points read in the integration point "
334 "initial conditions set function for {:s} variable.",
337 position += integration_points_read * ip_meta_data.n_components;
344 *_local_to_global_index_map);
347 template <
int DisplacementDim>
350 if (_use_monolithic_scheme)
352 const int monolithic_process_id = 0;
353 initializeProcessBoundaryConditionsAndSourceTerms(
354 *_local_to_global_index_map, monolithic_process_id);
360 const int hydraulic_process_id = 0;
361 initializeProcessBoundaryConditionsAndSourceTerms(
362 *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
365 const int mechanical_process_id = 1;
366 initializeProcessBoundaryConditionsAndSourceTerms(
367 *_local_to_global_index_map, mechanical_process_id);
370 template <
int DisplacementDim>
374 int const process_id)
381 DBUG(
"SetInitialConditions RichardsMechanicsProcess.");
388 _use_monolithic_scheme, process_id);
391 template <
int DisplacementDim>
393 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
394 std::vector<GlobalVector*>
const& xdot,
int const process_id,
397 DBUG(
"Assemble the equations for RichardsMechanics");
399 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
400 dof_table = {std::ref(*_local_to_global_index_map)};
410 template <
int DisplacementDim>
413 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
414 std::vector<GlobalVector*>
const& xdot,
const double dxdot_dx,
415 const double dx_dx,
int const process_id,
GlobalMatrix& M,
418 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
421 if (_use_monolithic_scheme)
424 "Assemble the Jacobian of RichardsMechanics for the monolithic"
426 dof_tables.emplace_back(*_local_to_global_index_map);
434 "Assemble the Jacobian equations of liquid fluid process in "
435 "RichardsMechanics for the staggered scheme.");
440 "Assemble the Jacobian equations of mechanical process in "
441 "RichardsMechanics for the staggered scheme.");
443 dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
444 dof_tables.emplace_back(*_local_to_global_index_map);
452 dxdot_dx, dx_dx, process_id, M, K, b, Jac);
454 auto copyRhs = [&](
int const variable_id,
auto& output_vector)
456 if (_use_monolithic_scheme)
460 std::negate<double>());
466 std::negate<double>());
469 if (_use_monolithic_scheme || process_id == 0)
471 copyRhs(0, *_hydraulic_flow);
473 if (_use_monolithic_scheme || process_id == 1)
475 copyRhs(1, *_nodal_forces);
479 template <
int DisplacementDim>
481 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
482 const int process_id)
484 if (hasMechanicalProcess(process_id))
486 DBUG(
"PostTimestep RichardsMechanicsProcess.");
487 auto const dof_tables = getDOFTables(x.size());
490 getProcessVariables(process_id)[0];
492 &LocalAssemblerIF::postTimestep, _local_assemblers,
497 template <
int DisplacementDim>
500 std::vector<GlobalVector*>
const& x,
502 int const process_id)
509 DBUG(
"Compute the secondary variables for RichardsMechanicsProcess.");
510 auto const dof_tables = getDOFTables(x.size());
514 &LocalAssemblerIF::computeSecondaryVariable, _local_assemblers,
518 template <
int DisplacementDim>
520 DisplacementDim>::getDOFTableForExtrapolatorData()
const
522 const bool manage_storage =
false;
523 return std::make_tuple(_local_to_global_index_map_single_component.get(),
527 template <
int DisplacementDim>
530 const int process_id)
const
532 if (hasMechanicalProcess(process_id))
534 return *_local_to_global_index_map;
538 return *_local_to_global_index_map_with_base_nodes;
541 template <
int DisplacementDim>
542 std::vector<NumLib::LocalToGlobalIndexMap const*>
544 int const number_of_processes)
const
546 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
547 dof_tables.reserve(number_of_processes);
548 std::generate_n(std::back_inserter(dof_tables), number_of_processes,
549 [&]() {
return &getDOFTable(dof_tables.size()); });
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
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::unique_ptr< IntegrationPointWriter > > _integration_point_writer
std::vector< std::unique_ptr< LocalAssemblerIF > > _local_assemblers
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
void initializeBoundaryConditions() override
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
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
void constructDofTable() 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_dot, int const process_id) override
std::vector< NumLib::LocalToGlobalIndexMap const * > getDOFTables(const int number_of_processes) const
MeshLib::PropertyVector< double > * _hydraulic_flow
RichardsMechanicsProcess(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, RichardsMechanicsProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
bool isLinear() const override
void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &x, double const t, int const process_id) override
MeshLib::PropertyVector< double > * _nodal_forces
void postTimestepConcreteProcess(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
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
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 > getSwellingStress() const =0
virtual std::vector< double > getSaturation() const =0
virtual std::vector< double > getSigma() const =0
virtual std::vector< double > getPorosity() const =0
virtual std::vector< double > getTransportPorosity() const =0
virtual std::vector< double > getEpsilon() const =0