Overview
Here we add the corresponding dynamic contributions to MyLaplacianElement
.
“Mass” matrix
There is no second term derivative for the temperature, so the corresponding “mass” matrix will be a zero matrix.
We need to add to my_laplacian_element.h
:
/**
* @brief This is called during the assembling process in order to calculate the elemental mass matrix
* @param rMassMatrix The elemental mass matrix
* @param rCurrentProcessInfo The current process info instance
*/
void CalculateMassMatrix(
MatrixType& rMassMatrix,
ProcessInfo& rCurrentProcessInfo
) override;
Then we will add the following to our my_laplacian_element.cpp
, this is a trivial code with the code already computed on the base class element.h
it will be enough, we are doing this for consistency.
void MyLaplacianElement::CalculateMassMatrix(
MatrixType& rMassMatrix,
ProcessInfo& rCurrentProcessInfo
)
{
KRATOS_TRY;
SizeType dimension = r_geom.WorkingSpaceDimension();
SizeType number_of_nodes = r_geom.size();
SizeType mat_size = dimension * number_of_nodes;
rMassMatrix = ZeroMatrix( mat_size, mat_size );
KRATOS_CATCH("");
}
“Damping” matrix
This “damping matrix” will be computed similar to a regular mass matrix, but using the specific heat instead of the density.
We need to add to my_laplacian_element.h
:
/**
* @brief This is called during the assembling process in order to calculate the elemental damping matrix
* @param rDampingMatrix The elemental damping matrix
* @param rCurrentProcessInfo The current process info instance
*/
void CalculateDampingMatrix(
MatrixType& rDampingMatrix,
ProcessInfo& rCurrentProcessInfo
) override;
Now we can add the following to my_laplacian_element.cpp
:
void MyLaplacianElement::CalculateDampingMatrix(
MatrixType& rDampingMatrix,
ProcessInfo& rCurrentProcessInfo
)
{
KRATOS_TRY;
const auto& r_geom = GetGeometry();
const auto& r_prop = GetProperties();
SizeType dimension = r_geom.WorkingSpaceDimension();
SizeType number_of_nodes = r_geom.size();
SizeType mat_size = dimension * number_of_nodes;
rDampingMatrix = ZeroMatrix( mat_size, mat_size );
KRATOS_ERROR_IF_NOT(r_prop.Has( SPECIFIC_HEAT )) << "SPECIFIC_HEAT has to be provided for the calculation of the DampingMatrix!" << std::endl;
const double specific_heat = r_prop[SPECIFIC_HEAT];
const double thickness = (dimension == 2 && r_prop.Has(THICKNESS)) ? r_prop[THICKNESS] : 1.0;
const bool compute_lumped_damped_matrix = rCurrentProcessInfo.Has(COMPUTE_LUMPED_MASS_MATRIX) ? rCurrentProcessInfo[COMPUTE_LUMPED_MASS_MATRIX] : false;
// LUMPED DAMPING MATRIX
if (compute_lumped_damped_matrix == true) {
const double total_specific_heat = GetGeometry().Volume() * specific_heat * thickness;
Vector lumping_factors;
lumping_factors = GetGeometry().LumpingFactors( lumping_factors );
for ( IndexType i = 0; i < number_of_nodes; ++i ) {
const double temp = lumping_factors[i] * total_specific_heat;
for ( IndexType j = 0; j < dimension; ++j ) {
IndexType index = i * dimension + j;
rDampingMatrix( index, index ) = temp;
}
}
} else { // CONSISTENT DAMPING
Matrix J0(dimension, dimension);
IntegrationMethod integration_method = IntegrationUtilities::GetIntegrationMethodForExactMassMatrixEvaluation(r_geom);
const GeometryType::IntegrationPointsArrayType& integration_points = r_geom.IntegrationPoints( integration_method );
const Matrix& Ncontainer = r_geom.ShapeFunctionsValues(integration_method);
for ( IndexType point_number = 0; point_number < integration_points.size(); ++point_number ) {
GeometryUtils::JacobianOnInitialConfiguration(r_geom, integration_points[point_number], J0);
const double detJ0 = MathUtils<double>::DetMat(J0);
const double integration_weight = GetIntegrationWeight(integration_points, point_number, detJ0) * thickness;
const Vector& rN = row(Ncontainer,point_number);
for ( IndexType i = 0; i < number_of_nodes; ++i ) {
const SizeType index_i = i * dimension;
for ( IndexType j = 0; j < number_of_nodes; ++j ) {
const SizeType index_j = j * dimension;
const double NiNj_weight = rN[i] * rN[j] * integration_weight * specific_heat;
for ( IndexType k = 0; k < dimension; ++k )
rDampingMatrix( index_i + k, index_j + k ) += NiNj_weight;
}
}
}
}
KRATOS_CATCH("");
}