Skip to content

Commit

Permalink
Merge pull request #12658 from KratosMultiphysics/structural/hotfix-m…
Browse files Browse the repository at this point in the history 
…ixed-u-eps-vol-rhs

[FastPR][Structural] Hotfixes in u-eps vol elements RHS
  • Loading branch information
@rubenzorrilla
rubenzorrilla authored Sep 9, 2024
2 parents 4ef6e05 + 6626fe0 commit 255ec37
Show file tree
Hide file tree
Showing 6 changed files with 68 additions and 40 deletions.
4 changes: 2 additions & 2 deletions ...MechanicsApplication/custom_elements/small_displacement_mixed_volumetric_strain_element.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -458,6 +458,8 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) SmallDisplacementMixedVolumet
///@name Protected member Variables
///@{

bool mIsDynamic; // Bool variable to indicate if the problem is dynamic

std::vector<Vector> mDisplacementSubscale1; /// Gauss points displacement subscale at previous time step

std::vector<Vector> mDisplacementSubscale2; /// Gauss points displacement subscale at two previous time step
Expand Down Expand Up @@ -698,8 +700,6 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) SmallDisplacementMixedVolumet
///@name Member Variables
///@{

bool mIsDynamic; // Bool variable to indicate if the problem is dynamic

Matrix mAnisotropyTensor; // The anisotropy transformation tensor

Matrix mInverseAnisotropyTensor; // The inverse of the anisotropy transformation tensor
Expand Down
26 changes: 22 additions & 4 deletions ...icsApplication/custom_elements/small_displacement_mixed_volumetric_strain_oss_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,23 @@ Element::Pointer SmallDisplacementMixedVolumetricStrainOssElement::Clone (
/***********************************************************************************/
/***********************************************************************************/

void SmallDisplacementMixedVolumetricStrainOssElement::Initialize(const ProcessInfo &rCurrentProcessInfo)
{
KRATOS_TRY

BaseType::Initialize(rCurrentProcessInfo);

// Deactivate the dynamic subscales as these are not implemented in the OSS elements yet
mIsDynamic = 0;
mDisplacementSubscale1.resize(0);
mDisplacementSubscale2.resize(0);

KRATOS_CATCH( "" )
}

/***********************************************************************************/
/***********************************************************************************/

void SmallDisplacementMixedVolumetricStrainOssElement::CalculateRightHandSide(
VectorType& rRightHandSideVector,
const ProcessInfo& rCurrentProcessInfo)
Expand Down Expand Up @@ -149,14 +166,14 @@ void SmallDisplacementMixedVolumetricStrainOssElement::CalculateRightHandSide(
gauss_point_auxiliary_variables);

// Assemble the current Gauss point OSS projection operator
// Note that this calculates the OSS stabilization operator to be applied to the RHS
CalculateOssStabilizationOperatorGaussPointContribution(
oss_proj_op,
kinematic_variables,
gauss_point_auxiliary_variables);
}

// Substract the OSS projections to the current residual
// Note that the OSS stabilization projection operator already includes the minus sign
VectorType proj_vect(matrix_size);
for (IndexType i_node = 0; i_node < n_nodes; ++i_node) {
const auto& r_us_proj = r_geometry[i_node].FastGetSolutionStepValue(DISPLACEMENT_PROJECTION);
Expand Down Expand Up @@ -246,6 +263,7 @@ void SmallDisplacementMixedVolumetricStrainOssElement::CalculateLocalSystem(
gauss_point_auxiliary_variables);

// Assemble the current Gauss point OSS projection operator
// Note that this calculates the OSS stabilization operator to be applied to the RHS
CalculateOssStabilizationOperatorGaussPointContribution(
oss_proj_op,
kinematic_variables,
Expand Down Expand Up @@ -467,7 +485,7 @@ void SmallDisplacementMixedVolumetricStrainOssElement::Calculate(
const auto& r_body_force = gauss_point_auxiliary_variables.BodyForce;
const auto& r_grad_eps = gauss_point_auxiliary_variables.VolumetricStrainGradient;
for (IndexType d = 0; d < dim; ++d) {
aux_proj[d] = -r_body_force[d] - gauss_point_auxiliary_variables.BulkModulus * r_grad_eps[d];
aux_proj[d] = r_body_force[d] + gauss_point_auxiliary_variables.BulkModulus * r_grad_eps[d];
}

// Nodal assembly of the projection contributions
Expand Down Expand Up @@ -570,10 +588,10 @@ void SmallDisplacementMixedVolumetricStrainOssElement::CalculateOssStabilization

for (IndexType j = 0; j < n_nodes; ++j) {
const double N_j = rThisKinematicVariables.N[j];
rOrthogonalSubScalesOperator(i*block_size + dim, j*block_size + dim) += aux_w_kappa_tau_2 * N_i * N_j;
rOrthogonalSubScalesOperator(i*block_size + dim, j*block_size + dim) -= aux_w_kappa_tau_2 * N_i * N_j;
for (IndexType d = 0; d < dim; ++d) {
rOrthogonalSubScalesOperator(i*block_size + dim, j*block_size + d) -= aux_w_kappa_tau_1 * G_i[d] * N_j;
rOrthogonalSubScalesOperator(i*block_size + d, j*block_size + dim) -= aux_w_kappa_tau_2 * psi_i[d] * N_j;
rOrthogonalSubScalesOperator(i*block_size + d, j*block_size + dim) += aux_w_kappa_tau_2 * psi_i[d] * N_j;
}
}
}
Expand Down
2 changes: 2 additions & 0 deletions ...anicsApplication/custom_elements/small_displacement_mixed_volumetric_strain_oss_element.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,6 +135,8 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) SmallDisplacementMixedVolumet
IndexType NewId,
NodesArrayType const& rThisNodes) const override;

void Initialize(const ProcessInfo &rCurrentProcessInfo) override;

void CalculateLocalSystem(
MatrixType& rLeftHandSideMatrix,
VectorType& rRightHandSideVector,
Expand Down
46 changes: 26 additions & 20 deletions ...ion/custom_elements/small_displacement_mixed_volumetric_strain_oss_non_linear_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,31 +240,34 @@ void SmallDisplacementMixedVolumetricStrainOssNonLinearElement::CalculateLeftHan
rCurrentProcessInfo);

// Call the base element OSS operator
// Note that this calculates the LHS of the OSS stabilization operator
MatrixType aux_oss_stab_operator(matrix_size, matrix_size);
CalculateOrthogonalSubScalesStabilizationOperator(
aux_oss_stab_operator,
rCurrentProcessInfo);

// Call the base element OSS operator
// Call the element OSS operator
// Note that this calculates the OSS operator to be applied to the RHS (transpose to the stabilization operator)
MatrixType aux_oss_operator(matrix_size, matrix_size);
CalculateOrthogonalSubScalesOperator(
aux_oss_operator,
rCurrentProcessInfo);

// Call the base element OSS lumped projection operator
// Call the element OSS lumped projection operator
// Note that this calculates the lumped projection operator to be applied to the RHS
MatrixType aux_lumped_mass_operator(matrix_size, matrix_size);
CalculateOrthogonalSubScalesLumpedProjectionOperator(
aux_lumped_mass_operator,
rCurrentProcessInfo);

// Assemble the extended modal analysis LHS
// Assemble the extended LHS
for (IndexType i = 0; i < matrix_size; ++i) {
for (IndexType j = 0; j < matrix_size; ++j) {
rLeftHandSideMatrix(i,j) = aux_stiffnes(i,j);
rLeftHandSideMatrix(i, matrix_size + j) = aux_oss_stab_operator(i,j);
rLeftHandSideMatrix(matrix_size + i, j) = aux_oss_operator(i,j);
rLeftHandSideMatrix(i, matrix_size + j) = -aux_oss_stab_operator(i,j);
rLeftHandSideMatrix(matrix_size + i, j) = -aux_oss_operator(i,j);
// rLeftHandSideMatrix(matrix_size + j, i) = aux_oss_stab_operator(i,j);
rLeftHandSideMatrix(matrix_size + i, matrix_size + j) = aux_lumped_mass_operator(i,j);
rLeftHandSideMatrix(matrix_size + i, matrix_size + j) = -aux_lumped_mass_operator(i,j);
}
}
}
Expand All @@ -288,20 +291,24 @@ void SmallDisplacementMixedVolumetricStrainOssNonLinearElement::CalculateRightHa
if (rRightHandSideVector.size() != matrix_size_full) {
rRightHandSideVector.resize(matrix_size_full, false);
}
rRightHandSideVector.clear();

// Call the base element to get the standard residual matrix
// Call the base element to get the standard residual vector
// Note that this already includes the substraction of the OSS projections
VectorType aux_rhs(block_size * n_nodes);
BaseType::CalculateRightHandSide(
aux_rhs,
rCurrentProcessInfo);

// Call the base element OSS operator
MatrixType aux_oss_stab_operator(matrix_size, matrix_size);
CalculateOrthogonalSubScalesStabilizationOperator(
aux_oss_stab_operator,
// Call the element OSS operator
// Note that this calculates the OSS operator to be applied to the RHS (transpose to the stabilization operator)
MatrixType aux_oss_operator(matrix_size, matrix_size);
CalculateOrthogonalSubScalesOperator(
aux_oss_operator,
rCurrentProcessInfo);

// Call the base element OSS lumped projection operator
// Call the element OSS lumped projection operator
// Note that this calculates the lumped projection operator to be applied to the RHS
MatrixType aux_lumped_mass_operator(matrix_size, matrix_size);
CalculateOrthogonalSubScalesLumpedProjectionOperator(
aux_lumped_mass_operator,
Expand Down Expand Up @@ -339,14 +346,13 @@ void SmallDisplacementMixedVolumetricStrainOssNonLinearElement::CalculateRightHa
}

// Calculate auxiliary arrays
const VectorType oss_stab_times_unk = prod(trans(aux_oss_stab_operator), aux_unk);
const VectorType oss_stab_times_proj = prod(aux_oss_stab_operator, aux_projection);
const VectorType oss_stab_times_unk = prod(aux_oss_operator, aux_unk);
const VectorType lumped_mass_times_proj = prod(aux_lumped_mass_operator, aux_projection);

// Assemble the extended modal analysis RHS
for (IndexType i = 0; i < matrix_size; ++i) {
rRightHandSideVector(i) = aux_rhs(i) - oss_stab_times_proj(i);
rRightHandSideVector(matrix_size + i) = - oss_stab_times_unk(i) - lumped_mass_times_proj(i);
rRightHandSideVector(i) = aux_rhs(i); // Note that this already includes the projections substraction
rRightHandSideVector(matrix_size + i) = lumped_mass_times_proj(i) + oss_stab_times_unk(i);
}
}

Expand Down Expand Up @@ -426,7 +432,7 @@ void SmallDisplacementMixedVolumetricStrainOssNonLinearElement::GetSecondDerivat
rValues.clear();

for (IndexType i_node = 0; i_node < n_nodes; ++i_node) {
const SizeType index = i_node * block_size_full;
const SizeType index = i_node * block_size;
const auto& r_acc = r_geometry[i_node].FastGetSolutionStepValue(ACCELERATION, Step);
for(IndexType d = 0; d < dim; ++d) {
rValues[index + d] = r_acc[d];
Expand Down Expand Up @@ -589,10 +595,10 @@ void SmallDisplacementMixedVolumetricStrainOssNonLinearElement::CalculateOrthogo

for (IndexType i = 0; i < n_nodes; ++i) {
const double N_i = kinematic_variables.N[i];
rOrthogonalSubScalesOperator(i * block_size + dim, j * block_size + dim) += aux_w_kappa_tau_2 * N_i * N_j; // Note that we multiply by kappa*tau_2 to symmetrize
rOrthogonalSubScalesOperator(i * block_size + dim, j * block_size + dim) -= aux_w_kappa_tau_2 * N_i * N_j; // Note that we multiply by kappa*tau_2 to symmetrize
for (IndexType d = 0; d < dim; ++d) {
rOrthogonalSubScalesOperator(i * block_size + d, j * block_size + dim) -= aux_w_kappa_tau_1 * N_i * G_j[d]; // Note that we multiply by tau_1 to symmetrize
rOrthogonalSubScalesOperator(i * block_size + dim, j * block_size + d) -= aux_w_kappa_tau_2 * N_i * psi_j[d]; // Note that we multiply by kappa*tau_2 to symmetrize
rOrthogonalSubScalesOperator(i * block_size + dim, j * block_size + d) += aux_w_kappa_tau_2 * N_i * psi_j[d]; // Note that we multiply by kappa*tau_2 to symmetrize
}
}
}
Expand Down Expand Up @@ -745,7 +751,7 @@ void SmallDisplacementMixedVolumetricStrainOssNonLinearElement::CalculateOrthogo
for (IndexType d = 0; d < dim; ++d) {
rOrthogonalSubScalesLumpedProjectionOperator(i * block_size + d, i * block_size + d) += aux_w_tau_1 * N_i * sum_N_j; // Note that we multiply by tau_1 to scale the projection with the symmetrization factors (see CalculateOrthogonalSubScalesOperator)
}
rOrthogonalSubScalesLumpedProjectionOperator(i * block_size + dim, i * block_size + dim) += aux_w_kappa_tau_2 * N_i * sum_N_j; // Note that we multiply by kappa*tau_2 to scale the projection with the symmetrization factors (see CalculateOrthogonalSubScalesOperator)
rOrthogonalSubScalesLumpedProjectionOperator(i * block_size + dim, i * block_size + dim) -= aux_w_kappa_tau_2 * N_i * sum_N_j; // Note that we multiply by kappa*tau_2 to scale the projection with the symmetrization factors (see CalculateOrthogonalSubScalesOperator)
}
}
}
Expand Down
13 changes: 7 additions & 6 deletions ...plication/tests/cpp_tests/test_small_displacement_mixed_volumetric_strain_oss_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ KRATOS_TEST_CASE_IN_SUITE(SmallDisplacementMixedVolumetricStrainOssElement2D3N,

// Check RHS and LHS results
const double tolerance = 1.0e-5;
const std::vector<double> expected_RHS({51173.0769231, 51173.0769231, -1816.54062688, -12711.5384615, -38461.5384615, 10555.3527785, -38461.5384615, -12711.5384615, 3972.72630992});
const std::vector<double> expected_RHS({51134.6153846, 51134.6153846, -1827.75857559, -12673.0769231, -38461.5384615, 10540.9297016, -38461.5384615, -12673.0769231, 3959.9057971});
const std::vector<double> expected_LHS_row_0({778846.153846, 9615.38461538, -317307.692308, -394230.769231, -384615.384615, -317307.692308, -384615.384615, 375000, -317307.692308});
KRATOS_CHECK_VECTOR_RELATIVE_NEAR(RHS, expected_RHS, tolerance)
KRATOS_CHECK_VECTOR_RELATIVE_NEAR(row(LHS,0), expected_LHS_row_0, tolerance)
Expand Down Expand Up @@ -162,8 +162,8 @@ KRATOS_TEST_CASE_IN_SUITE(SmallDisplacementMixedVolumetricStrainOssElement2D3NDy

// Check RHS and LHS results
const double tolerance = 1.0e-5;
const std::vector<double> expected_RHS({51181.5777242, 51173.076373, -1831.8931645, -12703.0376604, -38461.5390116, 10570.705396, -38453.0376604, -12711.5390116, 3972.72623009});
const std::vector<double> expected_LHS_row_0({778846.153846, 9615.38461538, -327036.58062, -394230.769231, -384615.384615, -307578.803995, -384615.384615, 375000, -317307.692308});
const std::vector<double> expected_RHS({51134.6153846, 51134.6153846, -1827.75857559, -12673.0769231, -38461.5384615, 10540.9297016, -38461.5384615, -12673.0769231, 3959.9057971});
const std::vector<double> expected_LHS_row_0({778846.153846, 9615.38461538, -317307.692308, -394230.769231, -384615.384615, -317307.692308, -384615.384615, 375000, -317307.692308});
const std::vector<double> expected_mass_row_0({83.3333333333, 0, 0, 41.6666666667, 0, 0, 41.6666666667, 0, 0});
const std::vector<double> expected_mass_row_2({0, 0, 0, 0, 0, 0, 0, 0, 0});
KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(RHS, expected_RHS, tolerance)
Expand Down Expand Up @@ -271,9 +271,10 @@ KRATOS_TEST_CASE_IN_SUITE(SmallDisplacementMixedVolumetricStrainOssElementZienki
// Check results
const double tolerance = 1.0e-6;
const double expected_vol_strain = 1.49650698603e-05;
const double expected_vol_strain_proj = -0.000531463341208;
const std::vector<double> expected_disp = {-0.000533062, 0.000562992, 0.0};
const std::vector<double> expected_disp_proj = {-2.20921422407, 2.20921422407, 0.0};
const double expected_vol_strain_proj = 9.66365079115e-11;
const std::vector<double> expected_disp = {-7.41543286634e-06, 3.73455725869e-05, 0.0};
const std::vector<double> expected_disp_proj = {2.98507494006, -2.98507494006, 0.0};

KRATOS_CHECK_VECTOR_NEAR(p_node_3->FastGetSolutionStepValue(DISPLACEMENT), expected_disp, tolerance);
KRATOS_CHECK_VECTOR_NEAR(p_node_3->FastGetSolutionStepValue(DISPLACEMENT_PROJECTION), expected_disp_proj, tolerance);
KRATOS_CHECK_NEAR(p_node_3->FastGetSolutionStepValue(VOLUMETRIC_STRAIN), expected_vol_strain, tolerance);
Expand Down
Loading

0 comments on commit 255ec37

Please sign in to comment.