26 template <
unsigned int TDim,
unsigned int TNumNodes>
40 :
Element(NewId, pGeometry, pProperties)
46 return make_intrusive<TransientThermalElement>(NewId, GetGeometry().
Create(rThisNodes), pProperties);
49 Element::Pointer
Create(
IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties)
const override
51 return make_intrusive<TransientThermalElement>(NewId, pGeom, pProperties);
58 rElementalDofList.resize(TNumNodes);
59 for (
unsigned int i = 0;
i < TNumNodes; ++
i) {
60 rElementalDofList[
i] = GetGeometry()[
i].pGetDof(TEMPERATURE);
70 rResult.resize(TNumNodes,
false);
71 for (
unsigned int i = 0;
i < TNumNodes; ++
i) {
72 rResult[
i] = GetGeometry()[
i].GetDof(TEMPERATURE).EquationId();
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);
93 AddContributionsToLhsMatrix(rLeftHandSideMatrix, conductivity_matrix, capacity_matrix,
94 rCurrentProcessInfo[DT_TEMPERATURE_COEFFICIENT]);
95 AddContributionsToRhsVector(rRightHandSideVector, conductivity_matrix, capacity_matrix);
116 KRATOS_ERROR <<
"Can't return integration method: unexpected "
118 << TNumNodes << std::endl;
127 CheckHasSolutionStepsDataFor(TEMPERATURE);
128 CheckHasSolutionStepsDataFor(DT_TEMPERATURE);
129 CheckHasDofsFor(TEMPERATURE);
131 CheckForNonZeroZCoordinateIn2D();
139 void CheckDomainSize()
const
141 constexpr
auto min_domain_size = 1.0e-15;
143 <<
"DomainSize smaller than " << min_domain_size <<
" for element " << Id() << std::endl;
146 void CheckHasSolutionStepsDataFor(
const Variable<double>& rVariable)
const
148 for (
const auto&
node : GetGeometry()) {
150 <<
"Missing variable " << rVariable.Name() <<
" on node " <<
node.Id() << std::endl;
154 void CheckHasDofsFor(
const Variable<double>& rVariable)
const
156 for (
const auto&
node : GetGeometry()) {
158 <<
"Missing degree of freedom for " << rVariable.Name() <<
" on node " <<
node.Id()
163 void CheckProperties()
const
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);
177 if constexpr (TDim == 3) {
178 CheckProperty(THERMAL_CONDUCTIVITY_SOLID_ZZ);
179 CheckProperty(THERMAL_CONDUCTIVITY_SOLID_YZ);
180 CheckProperty(THERMAL_CONDUCTIVITY_SOLID_XZ);
187 << rVariable.
Name() <<
" does not exist in the thermal element's properties" << std::endl;
189 << rVariable.
Name() <<
" has an invalid value at element " << Id() << std::endl;
195 << rVariable.
Name() <<
" does not exist in the thermal element's properties" << std::endl;
197 << rVariable.
Name() <<
" has a value of (" << GetProperties()[rVariable]
198 <<
") instead of (" << rName <<
") at element " << Id() << std::endl;
201 void CheckForNonZeroZCoordinateIn2D()
const
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; });
208 <<
" Node with non-zero Z coordinate found. Id: " << pos->Id() << std::endl;
212 static void AddContributionsToLhsMatrix(
MatrixType& rLeftHandSideMatrix,
213 const BoundedMatrix<double, TNumNodes, TNumNodes>& rConductivityMatrix,
214 const BoundedMatrix<double, TNumNodes, TNumNodes>& rCapacityMatrix,
215 double DtTemperatureCoefficient)
217 rLeftHandSideMatrix = rConductivityMatrix;
218 rLeftHandSideMatrix += (DtTemperatureCoefficient * rCapacityMatrix);
221 void AddContributionsToRhsVector(
VectorType& rRightHandSideVector,
222 const BoundedMatrix<double, TNumNodes, TNumNodes>& rConductivityMatrix,
223 const BoundedMatrix<double, TNumNodes, TNumNodes>& rCapacityMatrix)
const
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;
233 Vector CalculateIntegrationCoefficients(
const Vector& rDetJContainer)
const
235 const auto& r_integration_points = GetGeometry().IntegrationPoints(GetIntegrationMethod());
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];
247 BoundedMatrix<double, TNumNodes, TNumNodes> CalculateConductivityMatrix(
249 const Vector& rIntegrationCoefficients,
250 const ProcessInfo& rCurrentProcessInfo)
const
252 GeoThermalDispersionLaw law{TDim};
253 const auto constitutive_matrix =
254 law.CalculateThermalDispersionMatrix(GetProperties(), rCurrentProcessInfo);
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];
269 BoundedMatrix<double, TNumNodes, TNumNodes> CalculateCapacityMatrix(
const Vector& rIntegrationCoefficients,
270 const ProcessInfo& rCurrentProcessInfo)
const
272 const auto& r_properties = GetProperties();
273 RetentionLaw::Parameters
parameters(r_properties, rCurrentProcessInfo);
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];
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];
293 array_1d<double, TNumNodes> GetNodalValuesOf(
const Variable<double>& rNodalVariable)
const
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);
305 void save(
Serializer& rSerializer)
const override
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
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