OGS
HTFEM.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <Eigen/Core>
14#include <vector>
15
17#include "HTProcessData.h"
28
29namespace ProcessLib
30{
31namespace HT
32{
33template <typename ShapeFunction, int GlobalDim>
35{
38
41
46
47public:
48 HTFEM(MeshLib::Element const& element,
49 std::size_t const local_matrix_size,
50 NumLib::GenericIntegrationMethod const& integration_method,
51 bool const is_axially_symmetric,
52 HTProcessData const& process_data,
53 const unsigned dof_per_node)
55 _element(element),
56 _process_data(process_data),
57 _integration_method(integration_method)
58 {
59 // This assertion is valid only if all nodal d.o.f. use the same shape
60 // matrices.
61 assert(local_matrix_size == ShapeFunction::NPOINTS * dof_per_node);
62 (void)local_matrix_size;
63 (void)dof_per_node;
64
65 unsigned const n_integration_points =
67 _ip_data.reserve(n_integration_points);
68
71
72 double const aperture_size = _process_data.aperture_size(0.0, pos)[0];
73
74 auto const shape_matrices =
76 GlobalDim>(element, is_axially_symmetric,
78
79 for (unsigned ip = 0; ip < n_integration_points; ip++)
80 {
81 _ip_data.emplace_back(
82 shape_matrices[ip].N, shape_matrices[ip].dNdx,
84 shape_matrices[ip].integralMeasure *
85 shape_matrices[ip].detJ * aperture_size);
86 }
87 }
88
89 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
90 const unsigned integration_point) const override
91 {
92 auto const& N = _ip_data[integration_point].N;
93
94 // assumes N is stored contiguously in memory
95 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
96 }
97
100 Eigen::Vector3d getFlux(MathLib::Point3d const& pnt_local_coords,
101 double const t,
102 std::vector<double> const& local_x) const override
103 {
104 // Eval shape matrices at given point
105 // Note: Axial symmetry is set to false here, because we only need dNdx
106 // here, which is not affected by axial symmetry.
107 auto const shape_matrices =
109 GlobalDim>(
110 _element, false /*is_axially_symmetric*/,
111 std::array{pnt_local_coords})[0];
112
114 pos.setElementID(this->_element.getID());
115
117
118 // local_x contains the local temperature and pressure values
119 double T_int_pt = 0.0;
120 double p_int_pt = 0.0;
121 NumLib::shapeFunctionInterpolate(local_x, shape_matrices.N, T_int_pt,
122 p_int_pt);
123
124 vars.temperature = T_int_pt;
125 vars.phase_pressure = p_int_pt;
126
127 auto const& medium =
129 auto const& liquid_phase = medium.phase("AqueousLiquid");
130
131 // TODO (naumov) Temporary value not used by current material models.
132 // Need extension of secondary variables interface.
133 double const dt = std::numeric_limits<double>::quiet_NaN();
134 // fetch permeability, viscosity, density
135 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
137 .value(vars, pos, t, dt));
138
139 auto const mu =
141 .template value<double>(vars, pos, t, dt);
142 GlobalDimMatrixType const K_over_mu = K / mu;
143
144 auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
145 &local_x[local_x.size() / 2], ShapeFunction::NPOINTS);
147 -K_over_mu * shape_matrices.dNdx * p_nodal_values;
148
149 if (this->_process_data.has_gravity)
150 {
151 auto const rho_w =
152 liquid_phase
154 .template value<double>(vars, pos, t, dt);
155 auto const b =
157 [this->_element.getID()];
158 q += K_over_mu * rho_w * b;
159 }
160
161 Eigen::Vector3d flux;
162 flux.head<GlobalDim>() = q;
163 return flux;
164 }
165
166protected:
169
171 std::vector<
173 Eigen::aligned_allocator<
176
178 MaterialPropertyLib::VariableArray const& vars, const double porosity,
179 const double fluid_density, const double specific_heat_capacity_fluid,
180 ParameterLib::SpatialPosition const& pos, double const t,
181 double const dt)
182 {
183 auto const& medium =
184 *_process_data.media_map.getMedium(this->_element.getID());
185 auto const& solid_phase = medium.phase("Solid");
186
187 auto const specific_heat_capacity_solid =
188 solid_phase
189 .property(
191 .template value<double>(vars, pos, t, dt);
192
193 auto const solid_density =
195 .template value<double>(vars, pos, t, dt);
196
197 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
198 fluid_density * specific_heat_capacity_fluid * porosity;
199 }
200
203 const double fluid_density, const double specific_heat_capacity_fluid,
204 const GlobalDimVectorType& velocity,
205 ParameterLib::SpatialPosition const& pos, double const t,
206 double const dt)
207 {
208 auto const& medium =
210
211 auto thermal_conductivity =
212 MaterialPropertyLib::formEigenTensor<GlobalDim>(
213 medium
214 .property(
216 .value(vars, pos, t, dt));
217
218 auto const thermal_dispersivity_transversal =
219 medium
221 thermal_transversal_dispersivity)
222 .template value<double>();
223
224 auto const thermal_dispersivity_longitudinal =
225 medium
227 thermal_longitudinal_dispersivity)
228 .template value<double>();
229
230 // Thermal conductivity is moved outside and zero matrix is passed
231 // instead due to multiplication with fluid's density times specific
232 // heat capacity.
233 return thermal_conductivity +
234 fluid_density * specific_heat_capacity_fluid *
237 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
238 velocity, 0 /* phi */, thermal_dispersivity_transversal,
239 thermal_dispersivity_longitudinal);
240 }
241
242 std::vector<double> const& getIntPtDarcyVelocityLocal(
243 const double t, std::vector<double> const& local_x,
244 std::vector<double>& cache) const
245 {
246 std::vector<double> local_p{
247 local_x.data() + pressure_index,
248 local_x.data() + pressure_index + pressure_size};
249 std::vector<double> local_T{
250 local_x.data() + temperature_index,
251 local_x.data() + temperature_index + temperature_size};
252
253 auto const n_integration_points =
255
256 cache.clear();
257 auto cache_mat = MathLib::createZeroedMatrix<
258 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
259 cache, GlobalDim, n_integration_points);
260
263
265
266 auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
267 &local_p[0], ShapeFunction::NPOINTS);
268
269 auto const& medium =
271 auto const& liquid_phase = medium.phase("AqueousLiquid");
272
273 for (unsigned ip = 0; ip < n_integration_points; ++ip)
274 {
275 auto const& ip_data = _ip_data[ip];
276 auto const& N = ip_data.N;
277 auto const& dNdx = ip_data.dNdx;
278
279 pos.setIntegrationPoint(ip);
280
281 double T_int_pt = 0.0;
282 double p_int_pt = 0.0;
283 NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
284 NumLib::shapeFunctionInterpolate(local_T, N, T_int_pt);
285
286 vars.temperature = T_int_pt;
287 vars.phase_pressure = p_int_pt;
288
289 // TODO (naumov) Temporary value not used by current material
290 // models. Need extension of secondary variables interface.
291 double const dt = std::numeric_limits<double>::quiet_NaN();
292 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
294 .value(vars, pos, t, dt));
295
296 auto const mu =
297 liquid_phase
299 .template value<double>(vars, pos, t, dt);
300 GlobalDimMatrixType const K_over_mu = K / mu;
301
302 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
303
305 {
306 auto const rho_w =
307 liquid_phase
309 .template value<double>(vars, pos, t, dt);
310 auto const b =
312 [_element.getID()];
313 // here it is assumed that the vector b is directed 'downwards'
314 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
315 }
316 }
317
318 return cache;
319 }
320
321protected:
322 static const int pressure_index = ShapeFunction::NPOINTS;
323 static const int pressure_size = ShapeFunction::NPOINTS;
324 static const int temperature_index = 0;
325 static const int temperature_size = ShapeFunction::NPOINTS;
326};
327
328} // namespace HT
329} // namespace ProcessLib
Phase const & phase(std::size_t index) const
Definition: Medium.cpp:33
Property const & property(PropertyType const &p) const
Definition: Phase.cpp:52
double getWeight() const
Definition: WeightedPoint.h:80
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void setElementID(std::size_t element_id)
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
Definition: HTFEM.h:45
std::vector< double > const & getIntPtDarcyVelocityLocal(const double t, std::vector< double > const &local_x, std::vector< double > &cache) const
Definition: HTFEM.h:242
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
Definition: HTFEM.h:37
NumLib::GenericIntegrationMethod const & _integration_method
Definition: HTFEM.h:170
static const int temperature_index
Definition: HTFEM.h:324
typename ShapeMatricesType::NodalVectorType NodalVectorType
Definition: HTFEM.h:39
static const int temperature_size
Definition: HTFEM.h:325
double getHeatEnergyCoefficient(MaterialPropertyLib::VariableArray const &vars, const double porosity, const double fluid_density, const double specific_heat_capacity_fluid, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
Definition: HTFEM.h:177
HTFEM(MeshLib::Element const &element, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, HTProcessData const &process_data, const unsigned dof_per_node)
Definition: HTFEM.h:48
static const int pressure_size
Definition: HTFEM.h:323
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
Definition: HTFEM.h:36
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
Definition: HTFEM.h:42
GlobalDimMatrixType getThermalConductivityDispersivity(MaterialPropertyLib::VariableArray const &vars, const double fluid_density, const double specific_heat_capacity_fluid, const GlobalDimVectorType &velocity, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
Definition: HTFEM.h:201
Eigen::Vector3d getFlux(MathLib::Point3d const &pnt_local_coords, double const t, std::vector< double > const &local_x) const override
Definition: HTFEM.h:100
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
Definition: HTFEM.h:89
MeshLib::Element const & _element
Definition: HTFEM.h:167
HTProcessData const & _process_data
Definition: HTFEM.h:168
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
Definition: HTFEM.h:44
static const int pressure_index
Definition: HTFEM.h:322
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType > > > _ip_data
Definition: HTFEM.h:175
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
Definition: HTFEM.h:40
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition: Interpolation.h:27
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
Eigen::MatrixXd computeHydrodynamicDispersion(NumericalStabilization const &stabilizer, std::size_t const element_id, Eigen::MatrixXd const &pore_diffusion_coefficient, Eigen::VectorXd const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
MatrixType< GlobalDim, ShapeFunction::NPOINTS > GlobalDimNodalMatrixType
MatrixType< GlobalDim, GlobalDim > GlobalDimMatrixType
VectorType< GlobalDim > GlobalDimVectorType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
ParameterLib::Parameter< double > const & aperture_size
Definition: HTProcessData.h:44
MaterialPropertyLib::MaterialSpatialDistributionMap media_map
Definition: HTProcessData.h:27
NumLib::NumericalStabilization stabilizer
Definition: HTProcessData.h:36
std::vector< Eigen::VectorXd > const projected_specific_body_force_vectors
Projected specific body force vector: R * R^T * b.
Definition: HTProcessData.h:39