OGS
RichardsFlowFEM.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <vector>
14
32
33namespace ProcessLib
34{
35namespace RichardsFlow
36{
37template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType,
38 typename NodalMatrixType>
40{
41 IntegrationPointData(NodalRowVectorType const& N_,
42 GlobalDimNodalMatrixType const& dNdx_,
43 double const& integration_weight_,
44 NodalMatrixType const mass_operator_)
45 : N(N_),
46 dNdx(dNdx_),
47 integration_weight(integration_weight_),
48 mass_operator(mass_operator_)
49 {
50 }
51
52 NodalRowVectorType const N;
53 GlobalDimNodalMatrixType const dNdx;
54 double const integration_weight;
55 NodalMatrixType const mass_operator;
56
58};
59const unsigned NUM_NODAL_DOF = 1;
60
64{
65public:
66 virtual std::vector<double> const& getIntPtSaturation(
67 const double t,
68 std::vector<GlobalVector*> const& x,
69 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
70 std::vector<double>& cache) const = 0;
71
72 virtual std::vector<double> const& getIntPtDarcyVelocity(
73 const double t,
74 std::vector<GlobalVector*> const& x,
75 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
76 std::vector<double>& cache) const = 0;
77};
78
79template <typename ShapeFunction, int GlobalDim>
81{
84
86 ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
87
95
96public:
98 MeshLib::Element const& element,
99 std::size_t const local_matrix_size,
100 NumLib::GenericIntegrationMethod const& integration_method,
101 bool is_axially_symmetric,
102 RichardsFlowProcessData const& process_data)
103 : _element(element),
104 _process_data(process_data),
105 _integration_method(integration_method),
107 std::vector<double>(_integration_method.getNumberOfPoints()))
108 {
109 // This assertion is valid only if all nodal d.o.f. use the same shape
110 // matrices.
111 assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
112 (void)local_matrix_size;
113
114 unsigned const n_integration_points =
116 _ip_data.reserve(n_integration_points);
117
118 auto const shape_matrices =
120 GlobalDim>(element, is_axially_symmetric,
122
123 for (unsigned ip = 0; ip < n_integration_points; ip++)
124 {
125 auto const& sm = shape_matrices[ip];
126 const double integration_factor = sm.integralMeasure * sm.detJ;
127 _ip_data.emplace_back(
128 sm.N, sm.dNdx,
130 integration_factor,
131 sm.N.transpose() * sm.N * integration_factor *
133 }
134 }
135
136 void assemble(double const t, double const dt,
137 std::vector<double> const& local_x,
138 std::vector<double> const& /*local_x_prev*/,
139 std::vector<double>& local_M_data,
140 std::vector<double>& local_K_data,
141 std::vector<double>& local_b_data) override
142 {
143 auto const local_matrix_size = local_x.size();
144 // This assertion is valid only if all nodal d.o.f. use the same shape
145 // matrices.
146 assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
147
148 auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
149 local_M_data, local_matrix_size, local_matrix_size);
150 auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
151 local_K_data, local_matrix_size, local_matrix_size);
152 auto local_b = MathLib::createZeroedVector<NodalVectorType>(
153 local_b_data, local_matrix_size);
154
155 unsigned const n_integration_points =
159
160 auto const& medium =
162 auto const& liquid_phase = medium.phase("AqueousLiquid");
164 vars.temperature =
165 medium
166 .property(
168 .template value<double>(vars, pos, t, dt);
169
170 for (unsigned ip = 0; ip < n_integration_points; ip++)
171 {
172 pos.setIntegrationPoint(ip);
173 double p_int_pt = 0.0;
174 NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, p_int_pt);
175
176 vars.liquid_phase_pressure = p_int_pt;
177 vars.capillary_pressure = -p_int_pt;
178
179 auto const permeability =
180 MaterialPropertyLib::formEigenTensor<GlobalDim>(
181 medium.property(MaterialPropertyLib::permeability)
182 .value(vars, pos, t, dt));
183
184 auto const porosity =
186 .template value<double>(vars, pos, t, dt);
187
188 double const Sw =
190 .template value<double>(vars, pos, t, dt);
191 _saturation[ip] = Sw;
192 vars.liquid_saturation = Sw;
193
194 double const dSw_dpc =
196 .template dValue<double>(
198 pos, t, dt);
199
200 auto const drhow_dp =
201 liquid_phase
203 .template dValue<double>(
204 vars,
206 pos, t, dt);
207 auto const storage =
209 .template value<double>(vars, pos, t, dt);
210 double const mass_mat_coeff =
211 storage * Sw + porosity * Sw * drhow_dp - porosity * dSw_dpc;
212
213 local_M.noalias() += mass_mat_coeff * _ip_data[ip].mass_operator;
214
215 double const k_rel =
216 medium
218 relative_permeability)
219 .template value<double>(vars, pos, t, dt);
220
221 auto const mu =
222 liquid_phase.property(MaterialPropertyLib::viscosity)
223 .template value<double>(vars, pos, t, dt);
224 local_K.noalias() += _ip_data[ip].dNdx.transpose() * permeability *
225 _ip_data[ip].dNdx *
226 _ip_data[ip].integration_weight * (k_rel / mu);
227
229 {
230 auto const rho_w =
231 liquid_phase
233 .template value<double>(vars, pos, t, dt);
234 auto const& body_force = _process_data.specific_body_force;
235 assert(body_force.size() == GlobalDim);
236 NodalVectorType gravity_operator =
237 _ip_data[ip].dNdx.transpose() * permeability * body_force *
238 _ip_data[ip].integration_weight;
239 local_b.noalias() += (k_rel / mu) * rho_w * gravity_operator;
240 }
241 }
243 {
244 for (int idx_ml = 0; idx_ml < local_M.cols(); idx_ml++)
245 {
246 double const mass_lump_val = local_M.col(idx_ml).sum();
247 local_M.col(idx_ml).setZero();
248 local_M(idx_ml, idx_ml) = mass_lump_val;
249 }
250 } // end of mass lumping
251 }
252
253 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
254 const unsigned integration_point) const override
255 {
256 auto const& N = _ip_data[integration_point].N;
257
258 // assumes N is stored contiguously in memory
259 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
260 }
261
262 std::vector<double> const& getIntPtSaturation(
263 const double /*t*/,
264 std::vector<GlobalVector*> const& /*x*/,
265 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
266 std::vector<double>& /*cache*/) const override
267 {
268 assert(!_saturation.empty());
269 return _saturation;
270 }
271
272 std::vector<double> const& getIntPtDarcyVelocity(
273 const double t,
274 std::vector<GlobalVector*> const& x,
275 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
276 std::vector<double>& cache) const override
277 {
278 // TODO (tf) Temporary value not used by current material models. Need
279 // extension of secondary variable interface
280 double const dt = std::numeric_limits<double>::quiet_NaN();
281
282 constexpr int process_id = 0; // monolithic scheme.
283 auto const indices =
284 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
285 assert(!indices.empty());
286 auto const local_x = x[process_id]->get(indices);
287
290
291 auto const& medium =
293 auto const& liquid_phase = medium.phase("AqueousLiquid");
294
296 vars.temperature =
297 medium
298 .property(
300 .template value<double>(vars, pos, t, dt);
301
302 auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
303 &local_x[0], ShapeFunction::NPOINTS);
304
305 unsigned const n_integration_points =
307
308 cache.clear();
309 auto cache_vec = MathLib::createZeroedMatrix<
310 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
311 cache, GlobalDim, n_integration_points);
312
313 for (unsigned ip = 0; ip < n_integration_points; ++ip)
314 {
315 double p_int_pt = 0.0;
316 NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, p_int_pt);
317 vars.liquid_phase_pressure = p_int_pt;
318 vars.capillary_pressure = -p_int_pt;
319
320 double const Sw =
322 .template value<double>(vars, pos, t, dt);
323 vars.liquid_saturation = Sw;
324
325 auto const permeability =
326 MaterialPropertyLib::formEigenTensor<GlobalDim>(
327 medium.property(MaterialPropertyLib::permeability)
328 .value(vars, pos, t, dt));
329
330 double const k_rel =
331 medium
333 relative_permeability)
334 .template value<double>(vars, pos, t, dt);
335 auto const mu =
336 liquid_phase.property(MaterialPropertyLib::viscosity)
337 .template value<double>(vars, pos, t, dt);
338 auto const K_mat_coeff = permeability * (k_rel / mu);
339 cache_vec.col(ip).noalias() =
340 -K_mat_coeff * _ip_data[ip].dNdx * p_nodal_values;
342 {
343 auto const rho_w =
344 liquid_phase
346 .template value<double>(vars, pos, t, dt);
347 auto const& body_force = _process_data.specific_body_force;
348 assert(body_force.size() == GlobalDim);
349 // here it is assumed that the vector body_force is directed
350 // 'downwards'
351 cache_vec.col(ip).noalias() += K_mat_coeff * rho_w * body_force;
352 }
353 }
354
355 return cache;
356 }
357
358private:
361
363 std::vector<
366 Eigen::aligned_allocator<IntegrationPointData<
369 std::vector<double> _saturation;
370};
371
372} // namespace RichardsFlow
373} // namespace ProcessLib
Definition of the PiecewiseLinearInterpolation class.
Phase const & phase(std::size_t index) const
Definition: Medium.cpp:33
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)
void setIntegrationPoint(unsigned integration_point)
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
NumLib::GenericIntegrationMethod const & _integration_method
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType > > > _ip_data
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
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 override
std::vector< double > const & getIntPtSaturation(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const override
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
RichardsFlowProcessData const & _process_data
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool is_axially_symmetric, RichardsFlowProcessData const &process_data)
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
typename ShapeMatricesType::NodalVectorType NodalVectorType
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
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
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< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
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)
const unsigned NUM_NODAL_DOF
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
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
IntegrationPointData(NodalRowVectorType const &N_, GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_, NodalMatrixType const mass_operator_)
MaterialPropertyLib::MaterialSpatialDistributionMap media_map