26namespace RichardsMechanics
28template <
int DisplacementDim>
32 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
33 std::vector<std::unique_ptr<ParameterLib::ParameterBase>>
const& parameters,
34 unsigned const integration_order,
35 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
39 bool const use_monolithic_scheme)
40 :
Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
41 integration_order, std::move(process_variables),
42 std::move(secondary_variables), use_monolithic_scheme),
43 _process_data(std::move(process_data))
55 std::make_unique<MeshLib::IntegrationPointWriter>(
61 std::make_unique<MeshLib::IntegrationPointWriter>(
62 "saturation_ip", 1 , integration_order,
66 std::make_unique<MeshLib::IntegrationPointWriter>(
67 "porosity_ip", 1 , integration_order,
71 std::make_unique<MeshLib::IntegrationPointWriter>(
72 "transport_porosity_ip", 1 , integration_order,
76 std::make_unique<MeshLib::IntegrationPointWriter>(
83 std::make_unique<MeshLib::IntegrationPointWriter>(
90template <
int DisplacementDim>
96template <
int DisplacementDim>
99 const int process_id)
const
103 if (_use_monolithic_scheme || process_id == 1)
105 auto const& l = *_local_to_global_index_map;
106 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
107 &l.getGhostIndices(), &this->_sparsity_pattern};
111 auto const& l = *_local_to_global_index_map_with_base_nodes;
112 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
113 &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
116template <
int DisplacementDim>
120 _mesh_subset_all_nodes =
121 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
124 _mesh_subset_base_nodes =
125 std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
129 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
130 *_mesh_subset_all_nodes};
131 _local_to_global_index_map_single_component =
132 std::make_unique<NumLib::LocalToGlobalIndexMap>(
133 std::move(all_mesh_subsets_single_component),
137 if (_use_monolithic_scheme)
140 std::vector<MeshLib::MeshSubset> all_mesh_subsets{
141 *_mesh_subset_base_nodes};
144 const int monolithic_process_id = 0;
145 std::generate_n(std::back_inserter(all_mesh_subsets),
146 getProcessVariables(monolithic_process_id)[1]
148 .getNumberOfGlobalComponents(),
149 [&]() {
return *_mesh_subset_all_nodes; });
151 std::vector<int>
const vec_n_components{1, DisplacementDim};
152 _local_to_global_index_map =
153 std::make_unique<NumLib::LocalToGlobalIndexMap>(
154 std::move(all_mesh_subsets), vec_n_components,
156 assert(_local_to_global_index_map);
161 const int process_id = 1;
162 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
163 std::generate_n(std::back_inserter(all_mesh_subsets),
164 getProcessVariables(process_id)[0]
166 .getNumberOfGlobalComponents(),
167 [&]() {
return *_mesh_subset_all_nodes; });
169 std::vector<int>
const vec_n_components{DisplacementDim};
170 _local_to_global_index_map =
171 std::make_unique<NumLib::LocalToGlobalIndexMap>(
172 std::move(all_mesh_subsets), vec_n_components,
177 std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
178 *_mesh_subset_base_nodes};
179 _local_to_global_index_map_with_base_nodes =
180 std::make_unique<NumLib::LocalToGlobalIndexMap>(
181 std::move(all_mesh_subsets_base_nodes),
186 *_local_to_global_index_map_with_base_nodes, _mesh);
188 assert(_local_to_global_index_map);
189 assert(_local_to_global_index_map_with_base_nodes);
193template <
int DisplacementDim>
197 unsigned const integration_order)
205 auto add_secondary_variable = [&](std::string
const& name,
206 int const num_components,
207 auto get_ip_values_function)
209 _secondary_variables.addSecondaryVariable(
213 std::move(get_ip_values_function)));
216 add_secondary_variable(
"sigma",
218 DisplacementDim>::RowsAtCompileTime,
219 &LocalAssemblerIF::getIntPtSigma);
221 add_secondary_variable(
"swelling_stress",
223 DisplacementDim>::RowsAtCompileTime,
224 &LocalAssemblerIF::getIntPtSwellingStress);
226 add_secondary_variable(
"epsilon",
228 DisplacementDim>::RowsAtCompileTime,
229 &LocalAssemblerIF::getIntPtEpsilon);
231 add_secondary_variable(
"velocity", DisplacementDim,
232 &LocalAssemblerIF::getIntPtDarcyVelocity);
234 add_secondary_variable(
"saturation", 1,
235 &LocalAssemblerIF::getIntPtSaturation);
237 add_secondary_variable(
"micro_saturation", 1,
238 &LocalAssemblerIF::getIntPtMicroSaturation);
240 add_secondary_variable(
"micro_pressure", 1,
241 &LocalAssemblerIF::getIntPtMicroPressure);
243 add_secondary_variable(
"porosity", 1, &LocalAssemblerIF::getIntPtPorosity);
245 add_secondary_variable(
"transport_porosity", 1,
246 &LocalAssemblerIF::getIntPtTransportPorosity);
248 add_secondary_variable(
"dry_density_solid", 1,
249 &LocalAssemblerIF::getIntPtDryDensitySolid);
256 add_secondary_variable);
260 _process_data.solid_materials, _local_assemblers,
261 _integration_point_writer, integration_order);
263 _process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
267 _process_data.element_porosity = MeshLib::getOrCreateMeshProperty<double>(
271 _process_data.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
275 DisplacementDim>::RowsAtCompileTime);
277 _process_data.pressure_interpolated =
278 MeshLib::getOrCreateMeshProperty<double>(
288 *_local_to_global_index_map);
291template <
int DisplacementDim>
293 std::map<
int, std::shared_ptr<MaterialPropertyLib::Medium>>
const& media)
295 if (_use_monolithic_scheme)
297 const int monolithic_process_id = 0;
298 initializeProcessBoundaryConditionsAndSourceTerms(
299 *_local_to_global_index_map, monolithic_process_id, media);
305 const int hydraulic_process_id = 0;
306 initializeProcessBoundaryConditionsAndSourceTerms(
307 *_local_to_global_index_map_with_base_nodes, hydraulic_process_id,
311 const int mechanical_process_id = 1;
312 initializeProcessBoundaryConditionsAndSourceTerms(
313 *_local_to_global_index_map, mechanical_process_id, media);
316template <
int DisplacementDim>
320 int const process_id)
327 DBUG(
"SetInitialConditions RichardsMechanicsProcess.");
332 &LocalAssemblerIF::setInitialConditions, _local_assemblers,
336template <
int DisplacementDim>
338 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
339 std::vector<GlobalVector*>
const& x_prev,
int const process_id,
342 DBUG(
"Assemble the equations for RichardsMechanics");
344 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
345 _local_to_global_index_map.get()};
355template <
int DisplacementDim>
358 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
359 std::vector<GlobalVector*>
const& x_prev,
int const process_id,
363 if (_use_monolithic_scheme)
366 "Assemble the Jacobian of RichardsMechanics for the monolithic"
375 "Assemble the Jacobian equations of liquid fluid process in "
376 "RichardsMechanics for the staggered scheme.");
381 "Assemble the Jacobian equations of mechanical process in "
382 "RichardsMechanics for the staggered scheme.");
388 auto const dof_tables = getDOFTables(x.size());
392 x_prev, process_id, M, K, b, Jac);
394 auto copyRhs = [&](
int const variable_id,
auto& output_vector)
396 if (_use_monolithic_scheme)
398 transformVariableFromGlobalVector(b, variable_id, *dof_tables[0],
400 std::negate<double>());
404 transformVariableFromGlobalVector(b, 0, *dof_tables[process_id],
406 std::negate<double>());
409 if (_use_monolithic_scheme || process_id == 0)
411 copyRhs(0, *_hydraulic_flow);
413 if (_use_monolithic_scheme || process_id == 1)
415 copyRhs(1, *_nodal_forces);
419template <
int DisplacementDim>
421 std::vector<GlobalVector*>
const& x,
422 std::vector<GlobalVector*>
const& x_prev,
double const t,
double const dt,
423 const int process_id)
425 if (hasMechanicalProcess(process_id))
427 DBUG(
"PostTimestep RichardsMechanicsProcess.");
430 getProcessVariables(process_id)[0];
432 &LocalAssemblerIF::postTimestep, _local_assemblers,
438template <
int DisplacementDim>
441 std::vector<GlobalVector*>
const& x,
443 int const process_id)
450 DBUG(
"Compute the secondary variables for RichardsMechanicsProcess.");
454 &LocalAssemblerIF::computeSecondaryVariable, _local_assemblers,
459template <
int DisplacementDim>
461 DisplacementDim>::getDOFTableForExtrapolatorData()
const
463 const bool manage_storage =
false;
464 return std::make_tuple(_local_to_global_index_map_single_component.get(),
468template <
int DisplacementDim>
471 const int process_id)
const
473 if (hasMechanicalProcess(process_id))
475 return *_local_to_global_index_map;
479 return *_local_to_global_index_map_with_base_nodes;
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Global vector based on Eigen vector.
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).
Properties & getProperties()
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::unique_ptr< MeshLib::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 &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) 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_prev, int const process_id) override
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 assembleWithJacobianConcreteProcess(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) 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
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< 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)
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.
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 > getSaturation() const =0
virtual std::vector< double > getSwellingStress() const =0
virtual std::vector< double > getPorosity() const =0
virtual std::vector< double > getEpsilon() const =0
virtual std::vector< double > getSigma() const =0
virtual std::vector< double > getTransportPorosity() const =0