OGS
SteadyStateDiffusionFEM.h
Go to the documentation of this file.
1 
11 #pragma once
12 
13 #include <vector>
14 
15 #include "MaterialLib/MPL/Medium.h"
25 #include "ParameterLib/Parameter.h"
29 
30 namespace ProcessLib
31 {
32 namespace SteadyStateDiffusion
33 {
34 const unsigned NUM_NODAL_DOF = 1;
35 
39 {
40 public:
41  virtual std::vector<double> const& getIntPtDarcyVelocity(
42  const double t,
43  std::vector<GlobalVector*> const& x,
44  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
45  std::vector<double>& cache) const = 0;
46 };
47 
48 template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
50 {
53 
55  ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
56 
60 
61 public:
65  std::size_t const /*local_matrix_size*/,
66  bool is_axially_symmetric,
67  unsigned const integration_order,
68  SteadyStateDiffusionData const& process_data)
69  : _element(element),
70  _process_data(process_data),
71  _integration_method(integration_order),
74  GlobalDim>(
75  element, is_axially_symmetric, _integration_method))
76  {
77  }
78 
79  void assemble(double const t, double const dt,
80  std::vector<double> const& local_x,
81  std::vector<double> const& /*local_xdot*/,
82  std::vector<double>& /*local_M_data*/,
83  std::vector<double>& local_K_data,
84  std::vector<double>& /*local_b_data*/) override
85  {
86  auto const local_matrix_size = local_x.size();
87  // This assertion is valid only if all nodal d.o.f. use the same shape
88  // matrices.
89  assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
90 
91  auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
92  local_K_data, local_matrix_size, local_matrix_size);
93 
94  unsigned const n_integration_points =
95  _integration_method.getNumberOfPoints();
96 
99 
100  auto const& medium =
101  *_process_data.media_map->getMedium(_element.getID());
103  vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
104  medium
105  .property(
107  .template value<double>(vars, pos, t, dt);
108 
109  for (unsigned ip = 0; ip < n_integration_points; ip++)
110  {
111  pos.setIntegrationPoint(ip);
112  auto const& sm = _shape_matrices[ip];
113  auto const& wp = _integration_method.getWeightedPoint(ip);
114 
115  double p_int_pt = 0.0;
116  NumLib::shapeFunctionInterpolate(local_x, sm.N, p_int_pt);
117  vars[static_cast<int>(
119  auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
121  .value(vars, pos, t, dt));
122 
123  local_K.noalias() += sm.dNdx.transpose() * k * sm.dNdx * sm.detJ *
124  sm.integralMeasure * wp.getWeight();
125  }
126  }
127 
130  Eigen::Vector3d getFlux(MathLib::Point3d const& p_local_coords,
131  double const t,
132  std::vector<double> const& local_x) const override
133  {
134  // TODO (tf) Temporary value not used by current material models. Need
135  // extension of getFlux interface
136  double const dt = std::numeric_limits<double>::quiet_NaN();
137 
138  // Eval shape matrices at given point
139  // Note: Axial symmetry is set to false here, because we only need dNdx
140  // here, which is not affected by axial symmetry.
141  auto const shape_matrices =
143  GlobalDim>(
144  _element, false /*is_axially_symmetric*/,
145  std::array{p_local_coords})[0];
146 
147  // fetch hydraulic conductivity
149  pos.setElementID(_element.getID());
150  auto const& medium =
151  *_process_data.media_map->getMedium(_element.getID());
152 
154  vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
155  medium
156  .property(
158  .template value<double>(vars, pos, t, dt);
159  double pressure = 0.0;
160  NumLib::shapeFunctionInterpolate(local_x, shape_matrices.N, pressure);
161  vars[static_cast<int>(MaterialPropertyLib::Variable::phase_pressure)] =
162  pressure;
163 
164  auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
166  .value(vars, pos, t, dt));
167 
168  Eigen::Vector3d flux(0.0, 0.0, 0.0);
169  flux.head<GlobalDim>() =
170  -k * shape_matrices.dNdx *
171  Eigen::Map<const NodalVectorType>(local_x.data(), local_x.size());
172 
173  return flux;
174  }
175 
176  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
177  const unsigned integration_point) const override
178  {
179  auto const& N = _shape_matrices[integration_point].N;
180 
181  // assumes N is stored contiguously in memory
182  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
183  }
184 
185  std::vector<double> const& getIntPtDarcyVelocity(
186  const double t,
187  std::vector<GlobalVector*> const& x,
188  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
189  std::vector<double>& cache) const override
190  {
191  // TODO (tf) Temporary value not used by current material models. Need
192  // extension of secondary variable interface.
193  double const dt = std::numeric_limits<double>::quiet_NaN();
194 
195  auto const n_integration_points =
196  _integration_method.getNumberOfPoints();
197 
198  int const process_id = 0; // monolithic scheme
199  auto const indices =
200  NumLib::getIndices(_element.getID(), *dof_table[process_id]);
201  assert(!indices.empty());
202  auto const local_x = x[process_id]->get(indices);
203  auto const local_x_vec =
204  MathLib::toVector<Eigen::Matrix<double, ShapeFunction::NPOINTS, 1>>(
205  local_x, ShapeFunction::NPOINTS);
206 
207  cache.clear();
208  auto cache_mat = MathLib::createZeroedMatrix<
209  Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
210  cache, GlobalDim, n_integration_points);
211 
213  pos.setElementID(_element.getID());
214 
215  auto const& medium =
216  *_process_data.media_map->getMedium(_element.getID());
217 
219  vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
220  medium
221  .property(
223  .template value<double>(vars, pos, t, dt);
224  double pressure = 0.0;
225  for (unsigned i = 0; i < n_integration_points; ++i)
226  {
227  pos.setIntegrationPoint(i);
229  pressure);
230  vars[static_cast<int>(
232 
233  auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
235  .value(vars, pos, t, dt));
236  // dimensions: (d x 1) = (d x n) * (n x 1)
237  cache_mat.col(i).noalias() =
238  -k * _shape_matrices[i].dNdx * local_x_vec;
239  }
240 
241  return cache;
242  }
243 
244 private:
247 
248  IntegrationMethod const _integration_method;
249  std::vector<ShapeMatrices, Eigen::aligned_allocator<ShapeMatrices>>
251 };
252 
253 } // namespace SteadyStateDiffusion
254 } // namespace ProcessLib
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:82
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
LocalAssemblerData(MeshLib::Element const &element, std::size_t const, bool is_axially_symmetric, unsigned const integration_order, SteadyStateDiffusionData const &process_data)
void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &, std::vector< double > &, std::vector< double > &local_K_data, std::vector< double > &) override
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
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< ShapeMatrices, Eigen::aligned_allocator< ShapeMatrices > > _shape_matrices
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
typename LocalAssemblerTraits::LocalVector NodalVectorType
typename LocalAssemblerTraits::LocalMatrix NodalMatrixType
Eigen::Vector3d getFlux(MathLib::Point3d const &p_local_coords, double const t, std::vector< double > const &local_x) const override
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
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
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Definition: VariableType.h:116
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
static const double t
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Definition: Interpolation.h:79
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)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
VectorType< GlobalDim > GlobalDimVectorType
std::unique_ptr< MaterialPropertyLib::MaterialSpatialDistributionMap > media_map
Matrix< NNodes *NodalDOF, NNodes *NodalDOF > LocalMatrix