OGS
TH2MFEM.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <memory>
14#include <vector>
15
33
34namespace ProcessLib
35{
36namespace TH2M
37{
40template <typename ShapeMatrixType>
42{
43 std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N_u;
44};
45
46template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
47 int DisplacementDim>
48class TH2MLocalAssembler : public LocalAssemblerInterface<DisplacementDim>
49{
50public:
53
56
57 template <int N>
58 using VectorType =
59 typename ShapeMatricesTypePressure::template VectorType<N>;
60
61 template <int M, int N>
62 using MatrixType =
63 typename ShapeMatricesTypePressure::template MatrixType<M, N>;
64
67
70
71 static int const KelvinVectorSize =
73 using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
74
75 static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
76 DisplacementDim,
78
80
83
85 MeshLib::Element const& e,
86 std::size_t const /*local_matrix_size*/,
87 NumLib::GenericIntegrationMethod const& integration_method,
88 bool const is_axially_symmetric,
90
91private:
94 std::string_view const name,
95 double const* values,
96 int const integration_order) override;
97
98 void setInitialConditionsConcrete(Eigen::VectorXd const local_x,
99 double const t,
100 int const process_id) override;
101
102 void assemble(double const /*t*/, double const /*dt*/,
103 std::vector<double> const& /*local_x*/,
104 std::vector<double> const& /*local_x_prev*/,
105 std::vector<double>& /*local_M_data*/,
106 std::vector<double>& /*local_K_data*/,
107 std::vector<double>& /*local_rhs_data*/) override;
108
109 void assembleWithJacobian(double const t, double const dt,
110 std::vector<double> const& local_x,
111 std::vector<double> const& local_x_prev,
112 std::vector<double>& /*local_M_data*/,
113 std::vector<double>& /*local_K_data*/,
114 std::vector<double>& local_rhs_data,
115 std::vector<double>& local_Jac_data) override;
116
117 void initializeConcrete() override
118 {
119 unsigned const n_integration_points =
121
122 for (unsigned ip = 0; ip < n_integration_points; ip++)
123 {
124 auto& ip_data = _ip_data[ip];
125
126 ParameterLib::SpatialPosition const x_position{
127 std::nullopt, this->element_.getID(), ip,
129 ShapeFunctionDisplacement,
131 ip_data.N_u))};
132
134 if (this->process_data_.initial_stress.value)
135 {
136 this->current_states_[ip].eff_stress_data.sigma.noalias() =
138 DisplacementDim>(
139 (*this->process_data_.initial_stress.value)(
140 std::numeric_limits<
141 double>::quiet_NaN() /* time independent */,
142 x_position));
143 }
144
145 double const t = 0; // TODO (naumov) pass t from top
146 auto& material_state = this->material_states_[ip];
147 this->solid_material_.initializeInternalStateVariables(
148 t, x_position, *material_state.material_state_variables);
149
150 ip_data.pushBackState();
151 material_state.pushBackState();
152 }
153
154 for (unsigned ip = 0; ip < n_integration_points; ip++)
155 {
156 this->prev_states_[ip] = this->current_states_[ip];
157 }
158 }
159
160 void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
161 Eigen::VectorXd const& /*local_x_prev*/,
162 double const /*t*/, double const /*dt*/,
163 int const /*process_id*/) override
164 {
165 unsigned const n_integration_points =
167
168 for (unsigned ip = 0; ip < n_integration_points; ip++)
169 {
170 _ip_data[ip].pushBackState();
171 this->material_states_[ip].pushBackState();
172 }
173
174 for (unsigned ip = 0; ip < n_integration_points; ip++)
175 {
176 this->prev_states_[ip] = this->current_states_[ip];
177 }
178 }
179
181 double const t, double const dt, Eigen::VectorXd const& local_x,
182 Eigen::VectorXd const& local_x_prev) override;
183
184 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
185 const unsigned integration_point) const override
186 {
187 auto const& N_u = _secondary_data.N_u[integration_point];
188
189 // assumes N is stored contiguously in memory
190 return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
191 }
192
193 std::vector<double> const& getIntPtDarcyVelocityGas(
194 const double t,
195 std::vector<GlobalVector*> const& x,
196 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
197 std::vector<double>& cache) const override;
198 std::vector<double> const& getIntPtDarcyVelocityLiquid(
199 const double t,
200 std::vector<GlobalVector*> const& x,
201 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
202 std::vector<double>& cache) const override;
203 std::vector<double> const& getIntPtDiffusionVelocityVapourGas(
204 const double t,
205 std::vector<GlobalVector*> const& x,
206 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
207 std::vector<double>& cache) const override;
208 std::vector<double> const& getIntPtDiffusionVelocityGasGas(
209 const double t,
210 std::vector<GlobalVector*> const& x,
211 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
212 std::vector<double>& cache) const override;
213 std::vector<double> const& getIntPtDiffusionVelocitySoluteLiquid(
214 const double t,
215 std::vector<GlobalVector*> const& x,
216 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
217 std::vector<double>& cache) const override;
218 std::vector<double> const& getIntPtDiffusionVelocityLiquidLiquid(
219 const double t,
220 std::vector<GlobalVector*> const& x,
221 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
222 std::vector<double>& cache) const override;
223
224 std::tuple<
225 std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>,
226 std::vector<
229 Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
230 double const t, double const dt,
232 models);
233
234 // TODO: Here is some refactoring potential. All secondary variables could
235 // be stored in some container to avoid defining one method for each
236 // variable.
237
238 virtual std::vector<double> const& getIntPtLiquidDensity(
239 const double /*t*/,
240 std::vector<GlobalVector*> const& /*x*/,
241 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
242 std::vector<double>& cache) const override
243 {
245 &IpData::rhoLR, cache);
246 }
247
248 virtual std::vector<double> const& getIntPtGasDensity(
249 const double /*t*/,
250 std::vector<GlobalVector*> const& /*x*/,
251 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
252 std::vector<double>& cache) const override
253 {
255 &IpData::rhoGR, cache);
256 }
257
258 virtual std::vector<double> const& getIntPtSolidDensity(
259 const double /*t*/,
260 std::vector<GlobalVector*> const& /*x*/,
261 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
262 std::vector<double>& cache) const override
263 {
265 &IpData::rhoSR, cache);
266 }
267
268 virtual std::vector<double> const& getIntPtVapourPressure(
269 const double /*t*/,
270 std::vector<GlobalVector*> const& /*x*/,
271 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
272 std::vector<double>& cache) const override
273 {
275 &IpData::pWGR, cache);
276 }
277
278 virtual std::vector<double> const& getIntPtPorosity(
279 const double /*t*/,
280 std::vector<GlobalVector*> const& /*x*/,
281 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
282 std::vector<double>& cache) const override
283 {
285 cache);
286 }
287
288 virtual std::vector<double> const& getIntPtMoleFractionGas(
289 const double /*t*/,
290 std::vector<GlobalVector*> const& /*x*/,
291 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
292 std::vector<double>& cache) const override
293 {
295 &IpData::xnCG, cache);
296 }
297 virtual std::vector<double> const& getIntPtMassFractionGas(
298 const double /*t*/,
299 std::vector<GlobalVector*> const& /*x*/,
300 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
301 std::vector<double>& cache) const override
302 {
304 &IpData::xmCG, cache);
305 }
306 virtual std::vector<double> const& getIntPtMassFractionLiquid(
307 const double /*t*/,
308 std::vector<GlobalVector*> const& /*x*/,
309 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
310 std::vector<double>& cache) const override
311 {
313 &IpData::xmWL, cache);
314 }
315
316 virtual std::vector<double> const& getIntPtEnthalpyGas(
317 const double /*t*/,
318 std::vector<GlobalVector*> const& /*x*/,
319 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
320 std::vector<double>& cache) const override
321 {
323 cache);
324 }
325 virtual std::vector<double> const& getIntPtEnthalpyLiquid(
326 const double /*t*/,
327 std::vector<GlobalVector*> const& /*x*/,
328 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
329 std::vector<double>& cache) const override
330 {
332 cache);
333 }
334 virtual std::vector<double> const& getIntPtEnthalpySolid(
335 const double /*t*/,
336 std::vector<GlobalVector*> const& /*x*/,
337 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
338 std::vector<double>& cache) const override
339 {
341 cache);
342 }
343
344private:
347 using IpData =
349 ShapeMatricesTypePressure, DisplacementDim,
350 ShapeFunctionDisplacement::NPOINTS>;
351 std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
352
354 typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
356
357 // The shape function of pressure has the same form with the shape function
358 // of temperature
359 static const int gas_pressure_index = 0;
360 static const int gas_pressure_size = ShapeFunctionPressure::NPOINTS;
361 static const int capillary_pressure_index = ShapeFunctionPressure::NPOINTS;
362 static const int capillary_pressure_size = ShapeFunctionPressure::NPOINTS;
363 static const int temperature_index = 2 * ShapeFunctionPressure::NPOINTS;
364 static const int temperature_size = ShapeFunctionPressure::NPOINTS;
365 static const int displacement_index = ShapeFunctionPressure::NPOINTS * 3;
366 static const int displacement_size =
367 ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
368
369 static const int C_index = 0;
370 static const int C_size = ShapeFunctionPressure::NPOINTS;
371 static const int W_index = ShapeFunctionPressure::NPOINTS;
372 static const int W_size = ShapeFunctionPressure::NPOINTS;
373};
374
375} // namespace TH2M
376} // namespace ProcessLib
377
378#include "TH2MFEM-impl.h"
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
virtual std::vector< double > const & getIntPtEnthalpyGas(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:316
void assemble(double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
std::size_t setIPDataInitialConditions(std::string_view const name, double const *values, int const integration_order) override
static constexpr auto & N_u_op
Definition TH2MFEM.h:75
virtual std::vector< double > const & getIntPtMoleFractionGas(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:288
static const int capillary_pressure_index
Definition TH2MFEM.h:361
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition TH2MFEM.h:51
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
Definition TH2MFEM.h:73
std::vector< double > const & getIntPtDiffusionVelocitySoluteLiquid(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
virtual std::vector< double > const & getIntPtMassFractionGas(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:297
virtual std::vector< double > const & getIntPtEnthalpyLiquid(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:325
typename ShapeMatricesTypePressure::template VectorType< N > VectorType
Definition TH2MFEM.h:58
virtual std::vector< double > const & getIntPtSolidDensity(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:258
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition TH2MFEM.h:355
virtual std::vector< double > const & getIntPtVapourPressure(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:268
void assembleWithJacobian(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
std::vector< double > const & getIntPtDiffusionVelocityGasGas(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
void postTimestepConcrete(Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
Definition TH2MFEM.h:160
virtual std::vector< double > const & getIntPtGasDensity(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:248
static const int capillary_pressure_size
Definition TH2MFEM.h:362
TH2MLocalAssembler(TH2MLocalAssembler &&)=delete
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 override
std::vector< double > const & getIntPtDiffusionVelocityLiquidLiquid(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
typename ShapeMatricesTypePressure::template MatrixType< M, N > MatrixType
Definition TH2MFEM.h:62
virtual std::vector< double > const & getIntPtLiquidDensity(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:238
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
Definition TH2MFEM.h:184
virtual std::vector< double > const & getIntPtMassFractionLiquid(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:306
static const int gas_pressure_index
Definition TH2MFEM.h:359
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition TH2MFEM.h:351
static int const KelvinVectorSize
Definition TH2MFEM.h:71
std::vector< double > const & getIntPtDiffusionVelocityVapourGas(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > BMatricesType
Definition TH2MFEM.h:345
virtual std::vector< double > const & getIntPtPorosity(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:278
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition TH2MFEM.h:54
virtual std::vector< double > const & getIntPtEnthalpySolid(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Definition TH2MFEM.h:334
TH2MLocalAssembler(TH2MLocalAssembler const &)=delete
std::tuple< std::vector< ConstitutiveRelations::ConstitutiveData< DisplacementDim > >, std::vector< ConstitutiveRelations::ConstitutiveTempData< DisplacementDim > > > updateConstitutiveVariables(Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, ConstitutiveRelations::ConstitutiveModels< DisplacementDim > const &models)
void computeSecondaryVariableConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev) override
typename ShapeMatricesTypePressure::GlobalDimVectorType GlobalDimVectorType
Definition TH2MFEM.h:68
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 override
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
Definition TH2MFEM.h:65
void setInitialConditionsConcrete(Eigen::VectorXd const local_x, double const t, int const process_id) override
static const int displacement_index
Definition TH2MFEM.h:365
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
constexpr Eigen::CwiseNullaryOp< EigenBlockMatrixViewFunctor< D, M >, typename EigenBlockMatrixViewFunctor< D, M >::Matrix > eigenBlockMatrixView(const Eigen::MatrixBase< M > &matrix)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
MatrixType< GlobalDim, GlobalDim > GlobalDimMatrixType
VectorType< GlobalDim > GlobalDimVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
NumLib::GenericIntegrationMethod const & integration_method_
TH2MProcessData< DisplacementDim > & process_data_
ConstitutiveRelations::SolidConstitutiveRelation< DisplacementDim > const & solid_material_
std::vector< typename ConstitutiveRelations::StatefulData< DisplacementDim > > current_states_
std::vector< typename ConstitutiveRelations::StatefulDataPrev< DisplacementDim > > prev_states_
std::vector< ConstitutiveRelations::MaterialStateData< DisplacementDim > > material_states_
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
Definition TH2MFEM.h:43