24 namespace ThermoHydroMechanics
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))
49 _heat_flux = MeshLib::getOrCreateMeshProperty<double>(
53 std::make_unique<IntegrationPointWriter>(
60 std::make_unique<IntegrationPointWriter>(
67 template <
int DisplacementDim>
73 template <
int DisplacementDim>
76 const int process_id)
const
80 if (_use_monolithic_scheme || process_id == 2)
82 auto const& l = *_local_to_global_index_map;
83 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
84 &l.getGhostIndices(), &this->_sparsity_pattern};
88 auto const& l = *_local_to_global_index_map_with_base_nodes;
89 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
90 &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
93 template <
int DisplacementDim>
97 _mesh_subset_all_nodes =
98 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
101 _mesh_subset_base_nodes =
102 std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
106 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
107 *_mesh_subset_all_nodes};
108 _local_to_global_index_map_single_component =
109 std::make_unique<NumLib::LocalToGlobalIndexMap>(
110 std::move(all_mesh_subsets_single_component),
114 if (_use_monolithic_scheme)
117 std::vector<MeshLib::MeshSubset> all_mesh_subsets{
118 *_mesh_subset_base_nodes};
121 all_mesh_subsets.push_back(*_mesh_subset_base_nodes);
124 const int monolithic_process_id = 0;
125 std::generate_n(std::back_inserter(all_mesh_subsets),
126 getProcessVariables(monolithic_process_id)[2]
128 .getNumberOfGlobalComponents(),
129 [&]() {
return *_mesh_subset_all_nodes; });
131 std::vector<int>
const vec_n_components{1, 1, DisplacementDim};
132 _local_to_global_index_map =
133 std::make_unique<NumLib::LocalToGlobalIndexMap>(
134 std::move(all_mesh_subsets), vec_n_components,
136 assert(_local_to_global_index_map);
141 const int process_id = 2;
142 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
143 std::generate_n(std::back_inserter(all_mesh_subsets),
144 getProcessVariables(process_id)[0]
146 .getNumberOfGlobalComponents(),
147 [&]() {
return *_mesh_subset_all_nodes; });
149 std::vector<int>
const vec_n_components{DisplacementDim};
150 _local_to_global_index_map =
151 std::make_unique<NumLib::LocalToGlobalIndexMap>(
152 std::move(all_mesh_subsets), vec_n_components,
157 std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
158 *_mesh_subset_base_nodes};
159 _local_to_global_index_map_with_base_nodes =
160 std::make_unique<NumLib::LocalToGlobalIndexMap>(
161 std::move(all_mesh_subsets_base_nodes),
166 *_local_to_global_index_map_with_base_nodes, _mesh);
168 assert(_local_to_global_index_map);
169 assert(_local_to_global_index_map_with_base_nodes);
173 template <
int DisplacementDim>
177 unsigned const integration_order)
179 const int mechanical_process_id = _use_monolithic_scheme ? 0 : 2;
180 const int deformation_variable_id = _use_monolithic_scheme ? 2 : 0;
185 getProcessVariables(mechanical_process_id)[deformation_variable_id]
187 .getShapeFunctionOrder(),
191 _secondary_variables.addSecondaryVariable(
194 DisplacementDim>::RowsAtCompileTime,
195 getExtrapolator(), _local_assemblers,
198 _secondary_variables.addSecondaryVariable(
201 DisplacementDim>::RowsAtCompileTime,
202 getExtrapolator(), _local_assemblers,
205 _secondary_variables.addSecondaryVariable(
211 _process_data.pressure_interpolated =
212 MeshLib::getOrCreateMeshProperty<double>(
216 _process_data.temperature_interpolated =
217 MeshLib::getOrCreateMeshProperty<double>(
218 const_cast<MeshLib::Mesh&
>(mesh),
"temperature_interpolated",
222 for (
auto const& ip_writer : _integration_point_writer)
225 auto const&
name = ip_writer->name();
230 auto const& mesh_property =
234 if (mesh_property.getMeshItemType() !=
243 if (ip_meta_data.n_components !=
244 mesh_property.getNumberOfGlobalComponents())
247 "Different number of components in meta data ({:d}) than in "
248 "the integration point field data for '{:s}': {:d}.",
249 ip_meta_data.n_components,
name,
250 mesh_property.getNumberOfGlobalComponents());
255 std::size_t position = 0;
256 for (
auto& local_asm : _local_assemblers)
258 std::size_t
const integration_points_read =
259 local_asm->setIPDataInitialConditions(
260 name, &mesh_property[position],
261 ip_meta_data.integration_order);
262 if (integration_points_read == 0)
265 "No integration points read in the integration point "
266 "initial conditions set function.");
268 position += integration_points_read * ip_meta_data.n_components;
275 *_local_to_global_index_map);
278 template <
int DisplacementDim>
280 DisplacementDim>::initializeBoundaryConditions()
282 if (_use_monolithic_scheme)
284 const int process_id_of_thermohydromechancs = 0;
285 initializeProcessBoundaryConditionsAndSourceTerms(
286 *_local_to_global_index_map, process_id_of_thermohydromechancs);
292 const int thermal_process_id = 0;
293 initializeProcessBoundaryConditionsAndSourceTerms(
294 *_local_to_global_index_map_with_base_nodes, thermal_process_id);
297 const int hydraulic_process_id = 1;
298 initializeProcessBoundaryConditionsAndSourceTerms(
299 *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
302 const int mechanical_process_id = 2;
303 initializeProcessBoundaryConditionsAndSourceTerms(
304 *_local_to_global_index_map, mechanical_process_id);
307 template <
int DisplacementDim>
309 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
310 std::vector<GlobalVector*>
const& xdot,
int const process_id,
313 DBUG(
"Assemble the equations for ThermoHydroMechanics");
315 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
316 dof_table = {std::ref(*_local_to_global_index_map)};
326 template <
int DisplacementDim>
329 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
330 std::vector<GlobalVector*>
const& xdot,
const double dxdot_dx,
331 const double dx_dx,
int const process_id,
GlobalMatrix& M,
334 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
337 if (_use_monolithic_scheme)
340 "Assemble the Jacobian of ThermoHydroMechanics for the monolithic "
342 dof_tables.emplace_back(*_local_to_global_index_map);
350 "Assemble the Jacobian equations of heat transport process in "
351 "ThermoHydroMechanics for the staggered scheme.");
353 else if (process_id == 1)
356 "Assemble the Jacobian equations of liquid fluid process in "
357 "ThermoHydroMechanics for the staggered scheme.");
362 "Assemble the Jacobian equations of mechanical process in "
363 "ThermoHydroMechanics for the staggered scheme.");
365 dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
366 dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
367 dof_tables.emplace_back(*_local_to_global_index_map);
375 dxdot_dx, dx_dx, process_id, M, K, b, Jac);
377 auto copyRhs = [&](
int const variable_id,
auto& output_vector)
379 if (_use_monolithic_scheme)
383 std::negate<double>());
389 std::negate<double>());
392 if (_use_monolithic_scheme || process_id == 0)
394 copyRhs(0, *_heat_flux);
396 if (_use_monolithic_scheme || process_id == 1)
398 copyRhs(1, *_hydraulic_flow);
400 if (_use_monolithic_scheme || process_id == 2)
402 copyRhs(2, *_nodal_forces);
406 template <
int DisplacementDim>
408 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
409 const int process_id)
411 DBUG(
"PreTimestep ThermoHydroMechanicsProcess.");
413 if (hasMechanicalProcess(process_id))
416 getProcessVariables(process_id)[0];
421 *x[process_id], t, dt);
425 template <
int DisplacementDim>
427 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
428 const int process_id)
435 DBUG(
"PostTimestep ThermoHydroMechanicsProcess.");
436 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
437 auto const n_processes = x.size();
438 dof_tables.reserve(n_processes);
439 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
441 dof_tables.push_back(&getDOFTable(process_id));
451 template <
int DisplacementDim>
455 double const dt,
const int process_id)
457 if (!hasMechanicalProcess(process_id))
462 DBUG(
"PostNonLinearSolver ThermoHydroMechanicsProcess.");
470 _use_monolithic_scheme, process_id);
473 template <
int DisplacementDim>
476 std::vector<GlobalVector*>
const& x,
478 const int process_id)
485 DBUG(
"Compute the secondary variables for ThermoHydroMechanicsProcess.");
486 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
487 auto const n_processes = x.size();
488 dof_tables.reserve(n_processes);
489 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
491 dof_tables.push_back(&getDOFTable(process_id));
500 template <
int DisplacementDim>
502 DisplacementDim>::getDOFTableForExtrapolatorData()
const
504 const bool manage_storage =
false;
505 return std::make_tuple(_local_to_global_index_map_single_component.get(),
509 template <
int DisplacementDim>
512 const int process_id)
const
514 if (hasMechanicalProcess(process_id))
516 return *_local_to_global_index_map;
520 return *_local_to_global_index_map_with_base_nodes;
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
virtual void postTimestep(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
virtual void preTimestep(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void initialize(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void computeSecondaryVariable(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, int const process_id)
void postNonLinearSolver(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, GlobalVector const &xdot, double const t, double const dt, bool const use_monolithic_scheme, int const process_id)
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.
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
void postNonLinearSolverConcreteProcess(GlobalVector const &x, GlobalVector const &xdot, const double t, double const dt, int const process_id) override
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
MeshLib::PropertyVector< double > * _hydraulic_flow
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double dt, int const process_id) 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 > * _heat_flux
MeshLib::PropertyVector< double > * _nodal_forces
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, const int process_id) override
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
bool isLinear() const override
ThermoHydroMechanicsProcess(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, ThermoHydroMechanicsProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
std::vector< std::unique_ptr< LocalAssemblerInterface > > _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 constructDofTable() 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)
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)
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 > getSigma() const =0
virtual std::vector< double > getEpsilon() const =0
virtual std::vector< double > const & getIntPtSigma(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual std::vector< double > const & getIntPtEpsilon(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0