KratosMultiphysics
KRATOS Multiphysics (Kratos) is a framework for building parallel, multi-disciplinary simulation software, aiming at modularity, extensibility, and high performance. Kratos is written in C++, and counts with an extensive Python interface.
transient_thermal_element.h
Go to the documentation of this file.
1 // KRATOS___
2 // // ) )
3 // // ___ ___
4 // // ____ //___) ) // ) )
5 // // / / // // / /
6 // ((____/ / ((____ ((___/ / MECHANICS
7 //
8 // License: geo_mechanics_application/license.txt
9 //
10 // Main authors: Mohamed Nabi
11 // John van Esch
12 // Gennady Markelov
13 //
14 
15 #pragma once
16 
20 #include "includes/element.h"
21 #include "includes/serializer.h"
22 
23 namespace Kratos
24 {
25 
26 template <unsigned int TDim, unsigned int TNumNodes>
27 class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientThermalElement : public Element
28 {
29 public:
31 
32  explicit TransientThermalElement(IndexType NewId = 0) : Element(NewId) {}
33 
34  TransientThermalElement(IndexType NewId, GeometryType::Pointer pGeometry)
35  : Element(NewId, pGeometry)
36  {
37  }
38 
39  TransientThermalElement(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties)
40  : Element(NewId, pGeometry, pProperties)
41  {
42  }
43 
44  Element::Pointer Create(IndexType NewId, const NodesArrayType& rThisNodes, PropertiesType::Pointer pProperties) const override
45  {
46  return make_intrusive<TransientThermalElement>(NewId, GetGeometry().Create(rThisNodes), pProperties);
47  }
48 
49  Element::Pointer Create(IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const override
50  {
51  return make_intrusive<TransientThermalElement>(NewId, pGeom, pProperties);
52  }
53 
54  void GetDofList(DofsVectorType& rElementalDofList, const ProcessInfo&) const override
55  {
57 
58  rElementalDofList.resize(TNumNodes);
59  for (unsigned int i = 0; i < TNumNodes; ++i) {
60  rElementalDofList[i] = GetGeometry()[i].pGetDof(TEMPERATURE);
61  }
62 
63  KRATOS_CATCH("")
64  }
65 
66  void EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo&) const override
67  {
69 
70  rResult.resize(TNumNodes, false);
71  for (unsigned int i = 0; i < TNumNodes; ++i) {
72  rResult[i] = GetGeometry()[i].GetDof(TEMPERATURE).EquationId();
73  }
74 
75  KRATOS_CATCH("")
76  }
77 
78  void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix,
79  VectorType& rRightHandSideVector,
80  const ProcessInfo& rCurrentProcessInfo) override
81  {
83 
85  Vector det_J_container;
86  GetGeometry().ShapeFunctionsIntegrationPointsGradients(dN_dX_container, det_J_container,
87  GetIntegrationMethod());
88  const auto integration_coefficients = CalculateIntegrationCoefficients(det_J_container);
89  const auto conductivity_matrix =
90  CalculateConductivityMatrix(dN_dX_container, integration_coefficients, rCurrentProcessInfo);
91  const auto capacity_matrix = CalculateCapacityMatrix(integration_coefficients, rCurrentProcessInfo);
92 
93  AddContributionsToLhsMatrix(rLeftHandSideMatrix, conductivity_matrix, capacity_matrix,
94  rCurrentProcessInfo[DT_TEMPERATURE_COEFFICIENT]);
95  AddContributionsToRhsVector(rRightHandSideVector, conductivity_matrix, capacity_matrix);
96 
97  KRATOS_CATCH("")
98  }
99 
101  {
102  switch (TNumNodes) {
103  case 3:
104  case 4:
105  case 6:
106  case 8:
107  case 9:
108  case 20:
109  case 27:
111  case 10:
113  case 15:
115  default:
116  KRATOS_ERROR << "Can't return integration method: unexpected "
117  "number of nodes: "
118  << TNumNodes << std::endl;
119  }
120  }
121 
122  int Check(const ProcessInfo&) const override
123  {
124  KRATOS_TRY
125 
126  CheckDomainSize();
127  CheckHasSolutionStepsDataFor(TEMPERATURE);
128  CheckHasSolutionStepsDataFor(DT_TEMPERATURE);
129  CheckHasDofsFor(TEMPERATURE);
130  CheckProperties();
131  CheckForNonZeroZCoordinateIn2D();
132 
133  KRATOS_CATCH("")
134 
135  return 0;
136  }
137 
138 private:
139  void CheckDomainSize() const
140  {
141  constexpr auto min_domain_size = 1.0e-15;
142  KRATOS_ERROR_IF(GetGeometry().DomainSize() < min_domain_size)
143  << "DomainSize smaller than " << min_domain_size << " for element " << Id() << std::endl;
144  }
145 
146  void CheckHasSolutionStepsDataFor(const Variable<double>& rVariable) const
147  {
148  for (const auto& node : GetGeometry()) {
149  KRATOS_ERROR_IF_NOT(node.SolutionStepsDataHas(rVariable))
150  << "Missing variable " << rVariable.Name() << " on node " << node.Id() << std::endl;
151  }
152  }
153 
154  void CheckHasDofsFor(const Variable<double>& rVariable) const
155  {
156  for (const auto& node : GetGeometry()) {
157  KRATOS_ERROR_IF_NOT(node.HasDofFor(rVariable))
158  << "Missing degree of freedom for " << rVariable.Name() << " on node " << node.Id()
159  << std::endl;
160  }
161  }
162 
163  void CheckProperties() const
164  {
165  CheckProperty(DENSITY_WATER);
166  CheckProperty(POROSITY);
167  CheckProperty(RETENTION_LAW, "SaturatedLaw");
168  CheckProperty(SATURATED_SATURATION);
169  CheckProperty(DENSITY_SOLID);
170  CheckProperty(SPECIFIC_HEAT_CAPACITY_WATER);
171  CheckProperty(SPECIFIC_HEAT_CAPACITY_SOLID);
172  CheckProperty(THERMAL_CONDUCTIVITY_WATER);
173  CheckProperty(THERMAL_CONDUCTIVITY_SOLID_XX);
174  CheckProperty(THERMAL_CONDUCTIVITY_SOLID_YY);
175  CheckProperty(THERMAL_CONDUCTIVITY_SOLID_XY);
176 
177  if constexpr (TDim == 3) {
178  CheckProperty(THERMAL_CONDUCTIVITY_SOLID_ZZ);
179  CheckProperty(THERMAL_CONDUCTIVITY_SOLID_YZ);
180  CheckProperty(THERMAL_CONDUCTIVITY_SOLID_XZ);
181  }
182  }
183 
184  void CheckProperty(const Kratos::Variable<double>& rVariable) const
185  {
186  KRATOS_ERROR_IF_NOT(GetProperties().Has(rVariable))
187  << rVariable.Name() << " does not exist in the thermal element's properties" << std::endl;
188  KRATOS_ERROR_IF(GetProperties()[rVariable] < 0.0)
189  << rVariable.Name() << " has an invalid value at element " << Id() << std::endl;
190  }
191 
192  void CheckProperty(const Kratos::Variable<std::string>& rVariable, const std::string& rName) const
193  {
194  KRATOS_ERROR_IF_NOT(GetProperties().Has(rVariable))
195  << rVariable.Name() << " does not exist in the thermal element's properties" << std::endl;
196  KRATOS_ERROR_IF_NOT(GetProperties()[rVariable] == rName)
197  << rVariable.Name() << " has a value of (" << GetProperties()[rVariable]
198  << ") instead of (" << rName << ") at element " << Id() << std::endl;
199  }
200 
201  void CheckForNonZeroZCoordinateIn2D() const
202  {
203  if constexpr (TDim == 2) {
204  const auto& r_geometry = GetGeometry();
205  auto pos = std::find_if(r_geometry.begin(), r_geometry.end(),
206  [](const auto& node) { return node.Z() != 0.0; });
207  KRATOS_ERROR_IF_NOT(pos == r_geometry.end())
208  << " Node with non-zero Z coordinate found. Id: " << pos->Id() << std::endl;
209  }
210  }
211 
212  static void AddContributionsToLhsMatrix(MatrixType& rLeftHandSideMatrix,
213  const BoundedMatrix<double, TNumNodes, TNumNodes>& rConductivityMatrix,
214  const BoundedMatrix<double, TNumNodes, TNumNodes>& rCapacityMatrix,
215  double DtTemperatureCoefficient)
216  {
217  rLeftHandSideMatrix = rConductivityMatrix;
218  rLeftHandSideMatrix += (DtTemperatureCoefficient * rCapacityMatrix);
219  }
220 
221  void AddContributionsToRhsVector(VectorType& rRightHandSideVector,
222  const BoundedMatrix<double, TNumNodes, TNumNodes>& rConductivityMatrix,
223  const BoundedMatrix<double, TNumNodes, TNumNodes>& rCapacityMatrix) const
224  {
225  const auto capacity_vector =
226  array_1d<double, TNumNodes>{-prod(rCapacityMatrix, GetNodalValuesOf(DT_TEMPERATURE))};
227  rRightHandSideVector = capacity_vector;
228  const auto conductivity_vector =
229  array_1d<double, TNumNodes>{-prod(rConductivityMatrix, GetNodalValuesOf(TEMPERATURE))};
230  rRightHandSideVector += conductivity_vector;
231  }
232 
233  Vector CalculateIntegrationCoefficients(const Vector& rDetJContainer) const
234  {
235  const auto& r_integration_points = GetGeometry().IntegrationPoints(GetIntegrationMethod());
236 
237  auto result = Vector{r_integration_points.size()};
238  for (unsigned int integration_point_index = 0;
239  integration_point_index < r_integration_points.size(); ++integration_point_index) {
240  result[integration_point_index] = r_integration_points[integration_point_index].Weight() *
241  rDetJContainer[integration_point_index];
242  }
243 
244  return result;
245  }
246 
247  BoundedMatrix<double, TNumNodes, TNumNodes> CalculateConductivityMatrix(
248  const GeometryType::ShapeFunctionsGradientsType& rShapeFunctionGradients,
249  const Vector& rIntegrationCoefficients,
250  const ProcessInfo& rCurrentProcessInfo) const
251  {
252  GeoThermalDispersionLaw law{TDim};
253  const auto constitutive_matrix =
254  law.CalculateThermalDispersionMatrix(GetProperties(), rCurrentProcessInfo);
255 
256  auto result = BoundedMatrix<double, TNumNodes, TNumNodes>{ZeroMatrix{TNumNodes, TNumNodes}};
257  for (unsigned int integration_point_index = 0;
258  integration_point_index < GetGeometry().IntegrationPointsNumber(GetIntegrationMethod());
259  ++integration_point_index) {
260  BoundedMatrix<double, TDim, TNumNodes> Temp =
261  prod(constitutive_matrix, trans(rShapeFunctionGradients[integration_point_index]));
262  result += prod(rShapeFunctionGradients[integration_point_index], Temp) *
263  rIntegrationCoefficients[integration_point_index];
264  }
265 
266  return result;
267  }
268 
269  BoundedMatrix<double, TNumNodes, TNumNodes> CalculateCapacityMatrix(const Vector& rIntegrationCoefficients,
270  const ProcessInfo& rCurrentProcessInfo) const
271  {
272  const auto& r_properties = GetProperties();
273  RetentionLaw::Parameters parameters(r_properties, rCurrentProcessInfo);
274  auto retention_law = RetentionLawFactory::Clone(r_properties);
275  const double saturation = retention_law->CalculateSaturation(parameters);
276  const auto c_water = r_properties[POROSITY] * saturation * r_properties[DENSITY_WATER] *
277  r_properties[SPECIFIC_HEAT_CAPACITY_WATER];
278  const auto c_solid = (1.0 - r_properties[POROSITY]) * r_properties[DENSITY_SOLID] *
279  r_properties[SPECIFIC_HEAT_CAPACITY_SOLID];
280 
281  auto result = BoundedMatrix<double, TNumNodes, TNumNodes>{ZeroMatrix{TNumNodes, TNumNodes}};
282  const auto& r_N_container = GetGeometry().ShapeFunctionsValues(GetIntegrationMethod());
283  for (unsigned int integration_point_index = 0;
284  integration_point_index < GetGeometry().IntegrationPointsNumber(GetIntegrationMethod());
285  ++integration_point_index) {
286  const auto N = Vector{row(r_N_container, integration_point_index)};
287  result += (c_water + c_solid) * outer_prod(N, N) * rIntegrationCoefficients[integration_point_index];
288  }
289 
290  return result;
291  }
292 
293  array_1d<double, TNumNodes> GetNodalValuesOf(const Variable<double>& rNodalVariable) const
294  {
295  auto result = array_1d<double, TNumNodes>{};
296  const auto& r_geometry = GetGeometry();
297  std::transform(r_geometry.begin(), r_geometry.end(), result.begin(), [&rNodalVariable](const auto& node) {
298  return node.FastGetSolutionStepValue(rNodalVariable);
299  });
300  return result;
301  }
302 
303  friend class Serializer;
304 
305  void save(Serializer& rSerializer) const override
306  {
308  }
309 
310  void load(Serializer& rSerializer) override
311  {
313  }
314 };
315 
316 } // namespace Kratos
Base class for all Elements.
Definition: element.h:60
std::vector< DofType::Pointer > DofsVectorType
Definition: element.h:100
std::vector< std::size_t > EquationIdVectorType
Definition: element.h:98
std::size_t IndexType
Definition: flags.h:74
IntegrationMethod
Definition: geometry_data.h:76
GeometryData::ShapeFunctionsGradientsType ShapeFunctionsGradientsType
Definition: geometry.h:189
PointerVector is a container like stl vector but using a vector to store pointers to its data.
Definition: pointer_vector.h:72
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
static unique_ptr< RetentionLaw > Clone(const Properties &rMaterialProperties)
Definition: retention_law_factory.h:45
The serialization consists in storing the state of an object into a storage format like data file or ...
Definition: serializer.h:123
Definition: transient_thermal_element.h:28
Element::Pointer Create(IndexType NewId, const NodesArrayType &rThisNodes, PropertiesType::Pointer pProperties) const override
It creates a new element pointer.
Definition: transient_thermal_element.h:44
Element::Pointer Create(IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const override
It creates a new element pointer.
Definition: transient_thermal_element.h:49
void CalculateLocalSystem(MatrixType &rLeftHandSideMatrix, VectorType &rRightHandSideVector, const ProcessInfo &rCurrentProcessInfo) override
Definition: transient_thermal_element.h:78
void EquationIdVector(EquationIdVectorType &rResult, const ProcessInfo &) const override
Definition: transient_thermal_element.h:66
GeometryData::IntegrationMethod GetIntegrationMethod() const override
Definition: transient_thermal_element.h:100
TransientThermalElement(IndexType NewId, GeometryType::Pointer pGeometry)
Definition: transient_thermal_element.h:34
TransientThermalElement(IndexType NewId=0)
Definition: transient_thermal_element.h:32
int Check(const ProcessInfo &) const override
Definition: transient_thermal_element.h:122
void GetDofList(DofsVectorType &rElementalDofList, const ProcessInfo &) const override
Definition: transient_thermal_element.h:54
TransientThermalElement(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties)
Definition: transient_thermal_element.h:39
KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(TransientThermalElement)
const std::string & Name() const
Definition: variable_data.h:201
#define KRATOS_SERIALIZE_SAVE_BASE_CLASS(Serializer, BaseType)
Definition: define.h:812
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_TRY
Definition: define.h:109
#define KRATOS_SERIALIZE_LOAD_BASE_CLASS(Serializer, BaseType)
Definition: define.h:815
#define KRATOS_ERROR
Definition: exception.h:161
#define KRATOS_ERROR_IF_NOT(conditional)
Definition: exception.h:163
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
Matrix MatrixType
Definition: geometrical_transformation_utilities.h:55
Modeler::Pointer Create(const std::string &ModelerName, Model &rModel, const Parameters ModelParameters)
Checks if the modeler is registered.
Definition: modeler_factory.cpp:30
bool Has(const std::string &ModelerName)
Checks if the modeler is registered.
Definition: modeler_factory.cpp:24
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Internals::Matrix< double, AMatrix::dynamic, 1 > Vector
Definition: amatrix_interface.h:472
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
array_1d< double, 3 > VectorType
Definition: solid_mechanics_application_variables.cpp:19
AMatrix::VectorOuterProductExpression< TExpression1Type, TExpression2Type > outer_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:582
AMatrix::MatrixRow< const TExpressionType > row(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression, std::size_t RowIndex)
Definition: amatrix_interface.h:649
AMatrix::TransposeMatrix< const T > trans(const T &TheMatrix)
Definition: amatrix_interface.h:486
Temp
Definition: PecletTest.py:105
parameters
Definition: fluid_chimera_analysis.py:35
def load(f)
Definition: ode_solve.py:307
N
Definition: sensitivityMatrix.py:29
integer i
Definition: TensorModule.f:17
Definition: mesh_converter.cpp:38