25 namespace ThermoRichardsMechanics
27 template <
int DisplacementDim>
31 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
32 std::vector<std::unique_ptr<ParameterLib::ParameterBase>>
const& parameters,
33 unsigned const integration_order,
34 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
38 bool const use_monolithic_scheme)
39 :
Process(std::move(
name), mesh, std::move(jacobian_assembler), parameters,
40 integration_order, std::move(process_variables),
41 std::move(secondary_variables), use_monolithic_scheme),
42 process_data_(std::move(process_data))
50 heat_flux_ = MeshLib::getOrCreateMeshProperty<double>(
57 std::make_unique<IntegrationPointWriter>(
63 std::make_unique<IntegrationPointWriter>(
64 "saturation_ip", 1 , integration_order,
68 std::make_unique<IntegrationPointWriter>(
69 "porosity_ip", 1 , integration_order,
73 std::make_unique<IntegrationPointWriter>(
74 "transport_porosity_ip", 1 , integration_order,
78 std::make_unique<IntegrationPointWriter>(
85 std::make_unique<IntegrationPointWriter>(
92 template <
int DisplacementDim>
98 template <
int DisplacementDim>
103 auto const& l = *_local_to_global_index_map;
104 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
105 &l.getGhostIndices(), &this->_sparsity_pattern};
108 template <
int DisplacementDim>
112 _mesh_subset_all_nodes =
113 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
116 mesh_subset_base_nodes_ =
117 std::make_unique<MeshLib::MeshSubset>(_mesh, base_nodes_);
121 std::vector<MeshLib::MeshSubset> all__meshsubsets_single_component{
122 *_mesh_subset_all_nodes};
123 local_to_global_index_map_single_component_ =
124 std::make_unique<NumLib::LocalToGlobalIndexMap>(
125 std::move(all__meshsubsets_single_component),
130 std::vector<MeshLib::MeshSubset> all__meshsubsets{*mesh_subset_base_nodes_};
133 all__meshsubsets.push_back(*mesh_subset_base_nodes_);
136 const int monolithic_process_id = 0;
137 std::generate_n(std::back_inserter(all__meshsubsets),
138 getProcessVariables(monolithic_process_id)[2]
140 .getNumberOfGlobalComponents(),
141 [&]() {
return *_mesh_subset_all_nodes; });
143 std::vector<int>
const vec_n_components{1, 1, DisplacementDim};
144 _local_to_global_index_map =
145 std::make_unique<NumLib::LocalToGlobalIndexMap>(
146 std::move(all__meshsubsets), vec_n_components,
148 assert(_local_to_global_index_map);
151 template <
int DisplacementDim>
155 unsigned const integration_order)
157 using nlohmann::json;
159 const int mechanical_process_id = 0;
160 const int deformation_variable_id = 2;
165 getProcessVariables(mechanical_process_id)[deformation_variable_id]
167 .getShapeFunctionOrder(),
171 auto add_secondary_variable = [&](std::string
const&
name,
172 int const num_components,
173 auto get_ip_values_function)
175 _secondary_variables.addSecondaryVariable(
179 std::move(get_ip_values_function)));
182 add_secondary_variable(
"sigma",
184 DisplacementDim>::RowsAtCompileTime,
185 &LocalAssemblerIF::getIntPtSigma);
187 add_secondary_variable(
"swelling_stress",
189 DisplacementDim>::RowsAtCompileTime,
190 &LocalAssemblerIF::getIntPtSwellingStress);
192 add_secondary_variable(
"epsilon",
194 DisplacementDim>::RowsAtCompileTime,
195 &LocalAssemblerIF::getIntPtEpsilon);
197 add_secondary_variable(
"velocity", DisplacementDim,
198 &LocalAssemblerIF::getIntPtDarcyVelocity);
200 add_secondary_variable(
"saturation", 1,
201 &LocalAssemblerIF::getIntPtSaturation);
203 add_secondary_variable(
"porosity", 1, &LocalAssemblerIF::getIntPtPorosity);
205 add_secondary_variable(
"transport_porosity", 1,
206 &LocalAssemblerIF::getIntPtTransportPorosity);
208 add_secondary_variable(
"dry_density_solid", 1,
209 &LocalAssemblerIF::getIntPtDryDensitySolid);
216 add_secondary_variable);
218 process_data_.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
222 process_data_.element_porosity = MeshLib::getOrCreateMeshProperty<double>(
226 process_data_.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
230 DisplacementDim>::RowsAtCompileTime);
232 process_data_.pressure_interpolated =
233 MeshLib::getOrCreateMeshProperty<double>(
236 process_data_.temperature_interpolated =
237 MeshLib::getOrCreateMeshProperty<double>(
238 const_cast<MeshLib::Mesh&
>(mesh),
"temperature_interpolated",
242 for (
auto const& ip_writer : _integration_point_writer)
245 auto const&
name = ip_writer->name();
250 auto const& _meshproperty =
254 if (_meshproperty.getMeshItemType() !=
263 if (ip_meta_data.n_components !=
264 _meshproperty.getNumberOfGlobalComponents())
267 "Different number of components in meta data ({:d}) than in "
268 "the integration point field data for '{:s}': {:d}.",
269 ip_meta_data.n_components,
name,
270 _meshproperty.getNumberOfGlobalComponents());
275 std::size_t position = 0;
276 for (
auto& local_asm : local_assemblers_)
278 std::size_t
const integration_points_read =
279 local_asm->setIPDataInitialConditions(
280 name, &_meshproperty[position],
281 ip_meta_data.integration_order);
282 if (integration_points_read == 0)
285 "No integration points read in the integration point "
286 "initial conditions set function.");
288 position += integration_points_read * ip_meta_data.n_components;
295 *_local_to_global_index_map);
298 template <
int DisplacementDim>
300 DisplacementDim>::initializeBoundaryConditions()
302 const int process_id = 0;
303 initializeProcessBoundaryConditionsAndSourceTerms(
304 *_local_to_global_index_map, process_id);
307 template <
int DisplacementDim>
311 int const process_id)
313 DBUG(
"SetInitialConditions ThermoRichardsMechanicsProcess.");
317 *_local_to_global_index_map, *x[process_id], t, _use_monolithic_scheme,
321 template <
int DisplacementDim>
323 const double ,
double const ,
324 std::vector<GlobalVector*>
const& ,
325 std::vector<GlobalVector*>
const& ,
int const ,
329 "The Picard method or the Newton-Raphson method with numerical "
330 "Jacobian is not implemented for ThermoRichardsMechanics with the full "
331 "monolithic coupling scheme");
334 template <
int DisplacementDim>
337 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
338 std::vector<GlobalVector*>
const& xdot,
const double dxdot_dx,
339 const double dx_dx,
int const process_id,
GlobalMatrix& M,
342 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
346 "Assemble the Jacobian of ThermoRichardsMechanics for the monolithic "
348 dof_tables.emplace_back(*_local_to_global_index_map);
355 dxdot_dx, dx_dx, process_id, M, K, b, Jac);
357 auto copyRhs = [&](
int const variable_id,
auto& output_vector)
360 output_vector, std::negate<double>());
363 copyRhs(0, *heat_flux_);
364 copyRhs(1, *hydraulic_flow_);
365 copyRhs(2, *nodal_forces_);
368 template <
int DisplacementDim>
371 double const t,
double const dt,
372 const int process_id)
374 DBUG(
"PostTimestep ThermoRichardsMechanicsProcess.");
376 auto const dof_tables = getDOFTables(x.size());
380 &LocalAssemblerIF::postTimestep, local_assemblers_,
384 template <
int DisplacementDim>
387 std::vector<GlobalVector*>
const& x,
389 int const process_id)
391 DBUG(
"Compute the secondary variables for ThermoRichardsMechanicsProcess.");
393 auto const dof_tables = getDOFTables(x.size());
398 &LocalAssemblerIF::computeSecondaryVariable, local_assemblers_,
402 template <
int DisplacementDim>
404 DisplacementDim>::getDOFTableForExtrapolatorData()
const
406 const bool manage_storage =
false;
407 return std::make_tuple(local_to_global_index_map_single_component_.get(),
411 template <
int DisplacementDim>
416 return *_local_to_global_index_map;
419 template <
int DisplacementDim>
420 std::vector<NumLib::LocalToGlobalIndexMap const*>
422 int const number_of_processes)
const
424 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
425 dof_tables.reserve(number_of_processes);
426 std::generate_n(std::back_inserter(dof_tables), number_of_processes,
427 [&]() {
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
Handles configuration of several secondary variables from the project file.
Global assembler for the monolithic scheme of the non-isothermal Richards flow coupled with mechanics...
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, int const process_id) override
bool isLinear() const override
void constructDofTable() 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
MeshLib::PropertyVector< double > * nodal_forces_
ThermoRichardsMechanicsProcess(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, ThermoRichardsMechanicsProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
MeshLib::PropertyVector< double > * heat_flux_
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
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::vector< NumLib::LocalToGlobalIndexMap const * > getDOFTables(const int number_of_processes) const
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &x, double const t, int const process_id) override
std::vector< std::unique_ptr< LocalAssemblerIF > > local_assemblers_
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
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)
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.
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 > getTransportPorosity() const =0
virtual std::vector< double > getPorosity() const =0
virtual std::vector< double > getEpsilon() const =0
virtual std::vector< double > getSaturation() const =0
virtual std::vector< double > getSigma() const =0
virtual std::vector< double > getSwellingStress() const =0