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>(
60 std::make_unique<IntegrationPointWriter>(
61 "saturation_ip", 1 , integration_order,
65 std::make_unique<IntegrationPointWriter>(
72 template <
int DisplacementDim>
78 template <
int DisplacementDim>
81 const int process_id)
const
85 if (_use_monolithic_scheme || process_id == deformation_process_id)
87 auto const& l = *_local_to_global_index_map;
88 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
89 &l.getGhostIndices(), &this->_sparsity_pattern};
93 auto const& l = *_local_to_global_index_map_with_base_nodes;
94 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
95 &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
98 template <
int DisplacementDim>
102 _mesh_subset_all_nodes =
103 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
106 _mesh_subset_base_nodes =
107 std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
111 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
112 *_mesh_subset_all_nodes};
113 _local_to_global_index_map_single_component =
114 std::make_unique<NumLib::LocalToGlobalIndexMap>(
115 std::move(all_mesh_subsets_single_component),
119 if (_use_monolithic_scheme)
122 std::vector<MeshLib::MeshSubset> all_mesh_subsets{
123 *_mesh_subset_base_nodes};
126 all_mesh_subsets.push_back(*_mesh_subset_base_nodes);
129 all_mesh_subsets.push_back(*_mesh_subset_base_nodes);
133 std::back_inserter(all_mesh_subsets),
134 getProcessVariables(monolithic_process_id)[deformation_process_id]
136 .getNumberOfGlobalComponents(),
137 [&]() {
return *_mesh_subset_all_nodes; });
139 std::vector<int>
const vec_n_components{
140 n_gas_pressure_components, n_capillary_pressure_components,
141 n_temperature_components, n_displacement_components};
143 _local_to_global_index_map =
144 std::make_unique<NumLib::LocalToGlobalIndexMap>(
145 std::move(all_mesh_subsets), vec_n_components,
147 assert(_local_to_global_index_map);
151 OGS_FATAL(
"A Staggered version of TH2M is not implemented.");
155 template <
int DisplacementDim>
159 unsigned const integration_order)
161 const int mechanical_process_id =
162 _use_monolithic_scheme ? 0 : deformation_process_id;
163 const int deformation_variable_id =
164 _use_monolithic_scheme ? deformation_process_id : 0;
169 getProcessVariables(mechanical_process_id)[deformation_variable_id]
171 .getShapeFunctionOrder(),
175 auto add_secondary_variable = [&](std::string
const&
name,
176 int const num_components,
177 auto get_ip_values_function)
179 _secondary_variables.addSecondaryVariable(
183 std::move(get_ip_values_function)));
186 add_secondary_variable(
"sigma",
188 DisplacementDim>::RowsAtCompileTime,
190 add_secondary_variable(
"epsilon",
192 DisplacementDim>::RowsAtCompileTime,
194 add_secondary_variable(
"velocity_gas", mesh.
getDimension(),
196 add_secondary_variable(
199 add_secondary_variable(
"saturation", 1,
201 add_secondary_variable(
"vapour_pressure", 1,
203 add_secondary_variable(
"porosity", 1,
205 add_secondary_variable(
"gas_density", 1,
207 add_secondary_variable(
"solid_density", 1,
209 add_secondary_variable(
"liquid_density", 1,
211 add_secondary_variable(
"mole_fraction_gas", 1,
213 add_secondary_variable(
"mass_fraction_gas", 1,
215 add_secondary_variable(
216 "mass_fraction_liquid", 1,
219 add_secondary_variable(
220 "relative_permeability_gas", 1,
222 add_secondary_variable(
223 "relative_permeability_liquid", 1,
226 add_secondary_variable(
"enthalpy_gas", 1,
229 add_secondary_variable(
"enthalpy_liquid", 1,
232 add_secondary_variable(
"enthalpy_solid", 1,
235 _process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
239 _process_data.gas_pressure_interpolated =
240 MeshLib::getOrCreateMeshProperty<double>(
241 const_cast<MeshLib::Mesh&
>(mesh),
"gas_pressure_interpolated",
244 _process_data.capillary_pressure_interpolated =
245 MeshLib::getOrCreateMeshProperty<double>(
246 const_cast<MeshLib::Mesh&
>(mesh),
"capillary_pressure_interpolated",
249 _process_data.liquid_pressure_interpolated =
250 MeshLib::getOrCreateMeshProperty<double>(
251 const_cast<MeshLib::Mesh&
>(mesh),
"liquid_pressure_interpolated",
254 _process_data.temperature_interpolated =
255 MeshLib::getOrCreateMeshProperty<double>(
256 const_cast<MeshLib::Mesh&
>(mesh),
"temperature_interpolated",
260 for (
auto const& ip_writer : _integration_point_writer)
263 auto const&
name = ip_writer->name();
268 auto const& _meshproperty =
272 if (_meshproperty.getMeshItemType() !=
281 if (ip_meta_data.n_components !=
282 _meshproperty.getNumberOfGlobalComponents())
285 "Different number of components in meta data ({:d}) than in "
286 "the integration point field data for '{:s}': {:d}.",
287 ip_meta_data.n_components,
name,
288 _meshproperty.getNumberOfGlobalComponents());
293 std::size_t position = 0;
294 for (
auto& local_asm : _local_assemblers)
296 std::size_t
const integration_points_read =
297 local_asm->setIPDataInitialConditions(
298 name, &_meshproperty[position],
299 ip_meta_data.integration_order);
300 if (integration_points_read == 0)
303 "No integration points read in the integration point "
304 "initial conditions set function.");
306 position += integration_points_read * ip_meta_data.n_components;
313 *_local_to_global_index_map);
316 template <
int DisplacementDim>
319 if (_use_monolithic_scheme)
321 initializeProcessBoundaryConditionsAndSourceTerms(
322 *_local_to_global_index_map, monolithic_process_id);
327 OGS_FATAL(
"A Staggered version of TH2M is not implemented.");
330 template <
int DisplacementDim>
332 std::vector<GlobalVector*>& x,
double const t,
int const process_id)
339 DBUG(
"Set initial conditions of TH2MProcess.");
343 getDOFTable(process_id), *x[process_id], t, _use_monolithic_scheme,
347 template <
int DisplacementDim>
349 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
350 std::vector<GlobalVector*>
const& xdot,
int const process_id,
353 DBUG(
"Assemble the equations for TH2M");
355 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
356 dof_table = {std::ref(*_local_to_global_index_map)};
366 template <
int DisplacementDim>
368 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
369 std::vector<GlobalVector*>
const& xdot,
const double dxdot_dx,
373 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
376 if (_use_monolithic_scheme)
378 DBUG(
"Assemble the Jacobian of TH2M for the monolithic scheme.");
379 dof_tables.emplace_back(*_local_to_global_index_map);
384 OGS_FATAL(
"A Staggered version of TH2M is not implemented.");
392 dxdot_dx, dx_dx, process_id, M, K, b, Jac);
394 auto copyRhs = [&](
int const variable_id,
auto& output_vector)
396 if (_use_monolithic_scheme)
400 std::negate<double>());
406 std::negate<double>());
409 if (_use_monolithic_scheme || process_id == 1)
411 copyRhs(0, *_hydraulic_flow);
413 if (_use_monolithic_scheme || process_id == 2)
415 copyRhs(1, *_nodal_forces);
419 template <
int DisplacementDim>
421 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
422 const int process_id)
424 DBUG(
"PreTimestep TH2MProcess.");
426 if (hasMechanicalProcess(process_id))
429 getProcessVariables(process_id)[0];
434 *x[process_id], t, dt);
438 template <
int DisplacementDim>
440 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
441 const int process_id)
443 DBUG(
"PostTimestep TH2MProcess.");
444 auto const dof_tables = getDOFTables(x.size());
451 template <
int DisplacementDim>
453 double const t,
double const dt, std::vector<GlobalVector*>
const& x,
461 DBUG(
"Compute the secondary variables for TH2MProcess.");
462 auto const dof_tables = getDOFTables(x.size());
470 template <
int DisplacementDim>
471 std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
474 const bool manage_storage =
false;
475 return std::make_tuple(_local_to_global_index_map_single_component.get(),
479 template <
int DisplacementDim>
481 const int process_id)
const
483 if (hasMechanicalProcess(process_id))
485 return *_local_to_global_index_map;
489 return *_local_to_global_index_map_with_base_nodes;
492 template <
int DisplacementDim>
493 std::vector<NumLib::LocalToGlobalIndexMap const*>
496 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
497 dof_tables.reserve(number_of_processes);
498 std::generate_n(std::back_inserter(dof_tables), number_of_processes,
499 [&]() {
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
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)
void setInitialConditions(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, bool const use_monolithic_scheme, int const process_id)
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)
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.
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
bool isLinear() const override
std::vector< std::unique_ptr< LocalAssemblerInterface > > _local_assemblers
std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData() const override
TH2MProcess(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, TH2MProcessData< 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
MeshLib::PropertyVector< double > * _nodal_forces
std::vector< NumLib::LocalToGlobalIndexMap const * > getDOFTables(int const number_of_processes) const
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double dt, int const) override
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
MeshLib::PropertyVector< double > * _hydraulic_flow
void constructDofTable() override
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, const int process_id) override
NumLib::LocalToGlobalIndexMap const & getDOFTable(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 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
void initializeBoundaryConditions() 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.
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 > const & getIntPtEnthalpyLiquid(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 & getIntPtMassFractionLiquid(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 & getIntPtLiquidDensity(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 & getIntPtMoleFractionGas(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 & getIntPtSolidDensity(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 & getIntPtPorosity(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 & getIntPtGasDensity(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 > getSaturation() const =0
virtual std::vector< double > getEpsilon() const =0
virtual std::vector< double > const & getIntPtRelativePermeabilityGas(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 & getIntPtDarcyVelocityLiquid(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 & getIntPtEnthalpyGas(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 & getIntPtRelativePermeabilityLiquid(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 > getSigma() const =0
virtual std::vector< double > const & getIntPtEnthalpySolid(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 & 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 & getIntPtVapourPressure(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 & getIntPtMassFractionGas(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 & getIntPtDarcyVelocityGas(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
virtual std::vector< double > const & getIntPtSaturation(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0