OGS
HeatConductionFEM.h
Go to the documentation of this file.
1 
11 #pragma once
12 
13 #include <vector>
14 
16 #include "MaterialLib/MPL/Medium.h"
27 
28 namespace ProcessLib
29 {
30 namespace HeatConduction
31 {
32 const unsigned NUM_NODAL_DOF = 1;
33 
37 {
38 public:
39  virtual std::vector<double> const& getIntPtHeatFluxX(
40  const double t,
41  std::vector<GlobalVector*> const& x,
42  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
43  std::vector<double>& cache) const = 0;
44 
45  virtual std::vector<double> const& getIntPtHeatFluxY(
46  const double t,
47  std::vector<GlobalVector*> const& x,
48  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
49  std::vector<double>& cache) const = 0;
50 
51  virtual std::vector<double> const& getIntPtHeatFluxZ(
52  const double t,
53  std::vector<GlobalVector*> const& x,
54  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
55  std::vector<double>& cache) const = 0;
56 };
57 
58 template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
60 {
63 
65  ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
66 
70 
71 public:
75  std::size_t const local_matrix_size,
76  bool is_axially_symmetric,
77  unsigned const integration_order,
78  HeatConductionProcessData const& process_data)
79  : _element(element),
80  _process_data(process_data),
81  _integration_method(integration_order),
84  GlobalDim>(
85  element, is_axially_symmetric, _integration_method)),
87  GlobalDim,
88  std::vector<double>(_integration_method.getNumberOfPoints()))
89  {
90  // This assertion is valid only if all nodal d.o.f. use the same shape
91  // matrices.
92  assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
93  (void)local_matrix_size;
94  }
95 
96  void assemble(double const t, double const dt,
97  std::vector<double> const& local_x,
98  std::vector<double> const& /*local_xdot*/,
99  std::vector<double>& local_M_data,
100  std::vector<double>& local_K_data,
101  std::vector<double>& /*local_b_data*/) override
102  {
103  auto const local_matrix_size = local_x.size();
104  // This assertion is valid only if all nodal d.o.f. use the same shape
105  // matrices.
106  assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
107 
108  auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
109  local_M_data, local_matrix_size, local_matrix_size);
110  auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
111  local_K_data, local_matrix_size, local_matrix_size);
112 
113  unsigned const n_integration_points =
114  _integration_method.getNumberOfPoints();
115 
117  pos.setElementID(_element.getID());
118 
119  auto const& medium =
120  *_process_data.media_map->getMedium(_element.getID());
122 
123  for (unsigned ip = 0; ip < n_integration_points; ip++)
124  {
125  pos.setIntegrationPoint(ip);
126  auto const& sm = _shape_matrices[ip];
127  auto const& wp = _integration_method.getWeightedPoint(ip);
128 
129  // get the local temperature and put it in the variable array for
130  // access in MPL
131  double T_int_pt = 0.0;
132  NumLib::shapeFunctionInterpolate(local_x, sm.N, T_int_pt);
133  vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
134  T_int_pt;
135 
136  auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
137  medium
138  .property(
140  .value(vars, pos, t, dt));
141  auto const specific_heat_capacity =
142  medium
145  .template value<double>(vars, pos, t, dt);
146  auto const density =
147  medium
148  .property(
150  .template value<double>(vars, pos, t, dt);
151 
152  local_K.noalias() += sm.dNdx.transpose() * k * sm.dNdx * sm.detJ *
153  wp.getWeight() * sm.integralMeasure;
154  local_M.noalias() += sm.N.transpose() * density *
155  specific_heat_capacity * sm.N * sm.detJ *
156  wp.getWeight() * sm.integralMeasure;
157  }
159  {
160  local_M = local_M.colwise().sum().eval().asDiagonal();
161  }
162  }
163 
164  void assembleWithJacobian(double const t, double const dt,
165  std::vector<double> const& local_x,
166  std::vector<double> const& local_xdot,
167  const double /*dxdot_dx*/, const double /*dx_dx*/,
168  std::vector<double>& /*local_M_data*/,
169  std::vector<double>& /*local_K_data*/,
170  std::vector<double>& local_rhs_data,
171  std::vector<double>& local_Jac_data) override
172  {
173  auto const local_matrix_size = local_x.size();
174  // This assertion is valid only if all nodal d.o.f. use the same shape
175  // matrices.
176  assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
177 
178  auto x = Eigen::Map<NodalVectorType const>(local_x.data(),
179  local_matrix_size);
180 
181  auto x_dot = Eigen::Map<NodalVectorType const>(local_xdot.data(),
182  local_matrix_size);
183 
184  auto local_Jac = MathLib::createZeroedMatrix<NodalMatrixType>(
185  local_Jac_data, local_matrix_size, local_matrix_size);
186  auto local_rhs = MathLib::createZeroedVector<NodalVectorType>(
187  local_rhs_data, local_matrix_size);
188 
189  NodalMatrixType laplace =
190  NodalMatrixType::Zero(local_matrix_size, local_matrix_size);
192  NodalMatrixType::Zero(local_matrix_size, local_matrix_size);
193 
194  unsigned const n_integration_points =
195  _integration_method.getNumberOfPoints();
196 
198  pos.setElementID(_element.getID());
199 
200  auto const& medium =
201  *_process_data.media_map->getMedium(_element.getID());
203 
204  for (unsigned ip = 0; ip < n_integration_points; ip++)
205  {
206  pos.setIntegrationPoint(ip);
207  auto const& sm = _shape_matrices[ip];
208  double const w =
209  _integration_method.getWeightedPoint(ip).getWeight() * sm.detJ *
210  sm.integralMeasure;
211 
212  // get the local temperature and put it in the variable array for
213  // access in MPL
214  double T_int_pt = 0.0;
215  NumLib::shapeFunctionInterpolate(local_x, sm.N, T_int_pt);
216  vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
217  T_int_pt;
218 
219  auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
220  medium
221  .property(
223  .value(vars, pos, t, dt));
224  auto const specific_heat_capacity =
225  medium
228  .template value<double>(vars, pos, t, dt);
229  auto const density =
230  medium
231  .property(
233  .template value<double>(vars, pos, t, dt);
234 
235  laplace.noalias() += sm.dNdx.transpose() * k * sm.dNdx * w;
236  storage.noalias() +=
237  sm.N.transpose() * density * specific_heat_capacity * sm.N * w;
238  }
240  {
241  storage = storage.colwise().sum().eval().asDiagonal();
242  }
243 
244  local_Jac.noalias() += laplace + storage / dt;
245  local_rhs.noalias() -= laplace * x + storage * x_dot;
246  }
247 
249  double const t, double const dt, Eigen::VectorXd const& local_x,
250  Eigen::VectorXd const& /*local_x_dot*/) override
251  {
252  unsigned const n_integration_points =
253  _integration_method.getNumberOfPoints();
254 
256  pos.setElementID(_element.getID());
257 
258  auto const& medium =
259  *_process_data.media_map->getMedium(_element.getID());
261 
262  for (unsigned ip = 0; ip < n_integration_points; ip++)
263  {
264  pos.setIntegrationPoint(ip);
265  auto const& sm = _shape_matrices[ip];
266  // get the local temperature and put it in the variable array for
267  // access in MPL
268  double T_int_pt = 0.0;
269  NumLib::shapeFunctionInterpolate(local_x, sm.N, T_int_pt);
270  vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
271  T_int_pt;
272 
273  auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
274  medium
275  .property(
277  .value(vars, pos, t, dt));
278  // heat flux only computed for output.
279  GlobalDimVectorType const heat_flux = -k * sm.dNdx * local_x;
280 
281  for (int d = 0; d < GlobalDim; ++d)
282  {
283  _heat_fluxes[d][ip] = heat_flux[d];
284  }
285  }
286  }
287 
288  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
289  const unsigned integration_point) const override
290  {
291  auto const& N = _shape_matrices[integration_point].N;
292 
293  // assumes N is stored contiguously in memory
294  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
295  }
296 
297  std::vector<double> const& getIntPtHeatFluxX(
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  {
303  assert(!_heat_fluxes.empty());
304  return _heat_fluxes[0];
305  }
306 
307  std::vector<double> const& getIntPtHeatFluxY(
308  const double /*t*/,
309  std::vector<GlobalVector*> const& /*x*/,
310  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
311  std::vector<double>& /*cache*/) const override
312  {
313  assert(_heat_fluxes.size() > 1);
314  return _heat_fluxes[1];
315  }
316 
317  std::vector<double> const& getIntPtHeatFluxZ(
318  const double /*t*/,
319  std::vector<GlobalVector*> const& /*x*/,
320  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
321  std::vector<double>& /*cache*/) const override
322  {
323  assert(_heat_fluxes.size() > 2);
324  return _heat_fluxes[2];
325  }
326 
327 private:
330 
331  IntegrationMethod const _integration_method;
332  std::vector<ShapeMatrices, Eigen::aligned_allocator<ShapeMatrices>>
334 
335  std::vector<std::vector<double>> _heat_fluxes;
336 };
337 
338 } // namespace HeatConduction
339 } // 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)
virtual std::vector< double > const & getIntPtHeatFluxY(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 & getIntPtHeatFluxX(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 & getIntPtHeatFluxZ(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
std::vector< std::vector< double > > _heat_fluxes
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, bool is_axially_symmetric, unsigned const integration_order, HeatConductionProcessData const &process_data)
std::vector< double > const & getIntPtHeatFluxY(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const override
typename LocalAssemblerTraits::LocalMatrix NodalMatrixType
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 > &) override
HeatConductionProcessData const & _process_data
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
std::vector< double > const & getIntPtHeatFluxZ(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const override
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
typename LocalAssemblerTraits::LocalVector NodalVectorType
void assembleWithJacobian(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_xdot, const double, const double, std::vector< double > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
void computeSecondaryVariableConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &) override
std::vector< ShapeMatrices, Eigen::aligned_allocator< ShapeMatrices > > _shape_matrices
std::vector< double > const & getIntPtHeatFluxX(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const override
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Definition: VariableType.h:108
static const double t
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Definition: Interpolation.h:79
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
bool const mass_lumping
If set mass lumping will be applied to the equation.
Matrix< NNodes *NodalDOF, NNodes *NodalDOF > LocalMatrix