Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Geo/12964 u quasistatic p euler #13022

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Geo/12964 u quasistatic p euler #13022

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
c415f6b
do not touch vector variables ( for velocity and acceleration ) in qu…
WPK4FEM Jan 20, 2025
30ac1f0
adapt unit tests for quasistatic euler scheme
WPK4FEM Jan 20, 2025
acf0f2e
min_iterations 10 to scale up all steps as before
WPK4FEM Jan 20, 2025
a03e61f
remove some not used parameters from project parameters and material …
WPK4FEM Jan 20, 2025
e1b1b8d
naming convention and (hopefully) Sonar Cloud codesmell indications i…
WPK4FEM Jan 20, 2025
8b17922
tell the derivation by this->
WPK4FEM Jan 20, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 17 additions & 1 deletion ...MechanicsApplication/custom_strategies/schemes/backward_euler_quasistatic_U_Pw_scheme.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@
#include "backward_euler_scheme.hpp"
#include "includes/define.h"
#include "includes/model_part.h"
#include "solving_strategies/schemes/scheme.h"
#include "utilities/parallel_utilities.h"

// Application includes
Expand All @@ -39,6 +38,23 @@ class BackwardEulerQuasistaticUPwScheme : public BackwardEulerScheme<TSparseSpac
}

protected:
inline void UpdateVariablesDerivatives(ModelPart& rModelPart) override
{
KRATOS_TRY

block_for_each(rModelPart.Nodes(), [this](Node& rNode) {
// static here means no velocities and accelerations for the displacement/rotation D.O.F.
for (const auto& r_first_order_scalar_variable : this->GetFirstOrderScalarVariables()) {
if (rNode.IsFixed(r_first_order_scalar_variable.first_time_derivative)) continue;

rNode.FastGetSolutionStepValue(r_first_order_scalar_variable.first_time_derivative) =
this->CalculateDerivative(r_first_order_scalar_variable.instance, rNode);
}
});

KRATOS_CATCH("")
}

std::string Info() const override { return "BackwardEulerQuasistaticUPwScheme"; }
}; // Class BackwardEulerQuasistaticUPwScheme

Expand Down
4 changes: 2 additions & 2 deletions ...eoMechanicsApplication/custom_strategies/schemes/geomechanics_time_integration_scheme.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -150,9 +150,9 @@ class GeoMechanicsTimeIntegrationScheme : public Scheme<TSparseSpace, TDenseSpac
return !(rComponent.IsDefined(ACTIVE)) || rComponent.Is(ACTIVE);
}

void FinalizeSolutionStep(ModelPart& rModelPart, TSystemMatrixType& A, TSystemVectorType& Dx, TSystemVectorType& b) override
void FinalizeSolutionStep(ModelPart& rModelPart, TSystemMatrixType& rA, TSystemVectorType& rDx, TSystemVectorType& rb) override
{
FinalizeSolutionStepActiveEntities(rModelPart, A, Dx, b);
FinalizeSolutionStepActiveEntities(rModelPart, rA, rDx, rb);
}

void CalculateSystemContributions(Element& rCurrentElement,
Expand Down
144 changes: 72 additions & 72 deletions ...cations/GeoMechanicsApplication/custom_strategies/schemes/newmark_dynamic_U_Pw_scheme.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,76 +121,76 @@ class NewmarkDynamicUPwScheme : public GeneralizedNewmarkScheme<TSparseSpace, TD
}

void CalculateSystemContributions(Condition& rCurrentCondition,
LocalSystemMatrixType& LHS_Contribution,
LocalSystemVectorType& RHS_Contribution,
Element::EquationIdVectorType& EquationId,
const ProcessInfo& CurrentProcessInfo) override
LocalSystemMatrixType& rLHS_Contribution,
LocalSystemVectorType& rRHS_Contribution,
Element::EquationIdVectorType& rEquationIds,
const ProcessInfo& rCurrentProcessInfo) override
{
KRATOS_TRY

int thread = OpenMPUtils::ThisThread();

rCurrentCondition.CalculateLocalSystem(LHS_Contribution, RHS_Contribution, CurrentProcessInfo);
rCurrentCondition.CalculateLocalSystem(rLHS_Contribution, rRHS_Contribution, rCurrentProcessInfo);

rCurrentCondition.CalculateMassMatrix(mMassMatrix[thread], CurrentProcessInfo);
rCurrentCondition.CalculateMassMatrix(mMassMatrix[thread], rCurrentProcessInfo);

rCurrentCondition.CalculateDampingMatrix(mDampingMatrix[thread], CurrentProcessInfo);
rCurrentCondition.CalculateDampingMatrix(mDampingMatrix[thread], rCurrentProcessInfo);

this->AddDynamicsToLHS(LHS_Contribution, mMassMatrix[thread], mDampingMatrix[thread], CurrentProcessInfo);
this->AddDynamicsToLHS(rLHS_Contribution, mMassMatrix[thread], mDampingMatrix[thread], rCurrentProcessInfo);

this->AddDynamicsToRHS(rCurrentCondition, RHS_Contribution, mMassMatrix[thread],
mDampingMatrix[thread], CurrentProcessInfo);
this->AddDynamicsToRHS(rCurrentCondition, rRHS_Contribution, mMassMatrix[thread],
mDampingMatrix[thread], rCurrentProcessInfo);

rCurrentCondition.EquationIdVector(EquationId, CurrentProcessInfo);
rCurrentCondition.EquationIdVector(rEquationIds, rCurrentProcessInfo);

KRATOS_CATCH("")
}

void CalculateSystemContributions(Element& rCurrentElement,
LocalSystemMatrixType& LHS_Contribution,
LocalSystemVectorType& RHS_Contribution,
Element::EquationIdVectorType& EquationId,
const ProcessInfo& CurrentProcessInfo) override
LocalSystemMatrixType& rLHS_Contribution,
LocalSystemVectorType& rRHS_Contribution,
Element::EquationIdVectorType& rEquationIds,
const ProcessInfo& rCurrentProcessInfo) override
{
KRATOS_TRY

int thread = OpenMPUtils::ThisThread();

rCurrentElement.CalculateLocalSystem(LHS_Contribution, RHS_Contribution, CurrentProcessInfo);
rCurrentElement.CalculateLocalSystem(rLHS_Contribution, rRHS_Contribution, rCurrentProcessInfo);

rCurrentElement.CalculateMassMatrix(mMassMatrix[thread], CurrentProcessInfo);
rCurrentElement.CalculateMassMatrix(mMassMatrix[thread], rCurrentProcessInfo);

rCurrentElement.CalculateDampingMatrix(mDampingMatrix[thread], CurrentProcessInfo);
rCurrentElement.CalculateDampingMatrix(mDampingMatrix[thread], rCurrentProcessInfo);

this->AddDynamicsToLHS(LHS_Contribution, mMassMatrix[thread], mDampingMatrix[thread], CurrentProcessInfo);
this->AddDynamicsToLHS(rLHS_Contribution, mMassMatrix[thread], mDampingMatrix[thread], rCurrentProcessInfo);

this->AddDynamicsToRHS(rCurrentElement, RHS_Contribution, mMassMatrix[thread],
mDampingMatrix[thread], CurrentProcessInfo);
this->AddDynamicsToRHS(rCurrentElement, rRHS_Contribution, mMassMatrix[thread],
mDampingMatrix[thread], rCurrentProcessInfo);

rCurrentElement.EquationIdVector(EquationId, CurrentProcessInfo);
rCurrentElement.EquationIdVector(rEquationIds, rCurrentProcessInfo);

KRATOS_CATCH("")
}

void CalculateRHSContribution(Element& rCurrentElement,
LocalSystemVectorType& RHS_Contribution,
Element::EquationIdVectorType& EquationId,
const ProcessInfo& CurrentProcessInfo) override
LocalSystemVectorType& rRHS_Contribution,
Element::EquationIdVectorType& rEquationIds,
const ProcessInfo& rCurrentProcessInfo) override
{
KRATOS_TRY

int thread = OpenMPUtils::ThisThread();

rCurrentElement.CalculateRightHandSide(RHS_Contribution, CurrentProcessInfo);
rCurrentElement.CalculateRightHandSide(rRHS_Contribution, rCurrentProcessInfo);

rCurrentElement.CalculateMassMatrix(mMassMatrix[thread], CurrentProcessInfo);
rCurrentElement.CalculateMassMatrix(mMassMatrix[thread], rCurrentProcessInfo);

rCurrentElement.CalculateDampingMatrix(mDampingMatrix[thread], CurrentProcessInfo);
rCurrentElement.CalculateDampingMatrix(mDampingMatrix[thread], rCurrentProcessInfo);

this->AddDynamicsToRHS(rCurrentElement, RHS_Contribution, mMassMatrix[thread],
mDampingMatrix[thread], CurrentProcessInfo);
this->AddDynamicsToRHS(rCurrentElement, rRHS_Contribution, mMassMatrix[thread],
mDampingMatrix[thread], rCurrentProcessInfo);

rCurrentElement.EquationIdVector(EquationId, CurrentProcessInfo);
rCurrentElement.EquationIdVector(rEquationIds, rCurrentProcessInfo);

KRATOS_CATCH("")
}
Expand Down Expand Up @@ -219,114 +219,114 @@ class NewmarkDynamicUPwScheme : public GeneralizedNewmarkScheme<TSparseSpace, TD
}

void CalculateLHSContribution(Condition& rCurrentCondition,
LocalSystemMatrixType& LHS_Contribution,
Element::EquationIdVectorType& EquationId,
const ProcessInfo& CurrentProcessInfo) override
LocalSystemMatrixType& rLHS_Contribution,
Element::EquationIdVectorType& rEquationIds,
const ProcessInfo& rCurrentProcessInfo) override
{
KRATOS_TRY

int thread = OpenMPUtils::ThisThread();

rCurrentCondition.CalculateLeftHandSide(LHS_Contribution, CurrentProcessInfo);
rCurrentCondition.CalculateLeftHandSide(rLHS_Contribution, rCurrentProcessInfo);

rCurrentCondition.CalculateMassMatrix(mMassMatrix[thread], CurrentProcessInfo);
rCurrentCondition.CalculateMassMatrix(mMassMatrix[thread], rCurrentProcessInfo);

rCurrentCondition.CalculateDampingMatrix(mDampingMatrix[thread], CurrentProcessInfo);
rCurrentCondition.CalculateDampingMatrix(mDampingMatrix[thread], rCurrentProcessInfo);

this->AddDynamicsToLHS(LHS_Contribution, mMassMatrix[thread], mDampingMatrix[thread], CurrentProcessInfo);
this->AddDynamicsToLHS(rLHS_Contribution, mMassMatrix[thread], mDampingMatrix[thread], rCurrentProcessInfo);

rCurrentCondition.EquationIdVector(EquationId, CurrentProcessInfo);
rCurrentCondition.EquationIdVector(rEquationIds, rCurrentProcessInfo);

KRATOS_CATCH("")
}

void CalculateLHSContribution(Element& rCurrentElement,
LocalSystemMatrixType& LHS_Contribution,
Element::EquationIdVectorType& EquationId,
const ProcessInfo& CurrentProcessInfo) override
LocalSystemMatrixType& rLHS_Contribution,
Element::EquationIdVectorType& rEquationIds,
const ProcessInfo& rCurrentProcessInfo) override
{
KRATOS_TRY

int thread = OpenMPUtils::ThisThread();

rCurrentElement.CalculateLeftHandSide(LHS_Contribution, CurrentProcessInfo);
rCurrentElement.CalculateLeftHandSide(rLHS_Contribution, rCurrentProcessInfo);

rCurrentElement.CalculateMassMatrix(mMassMatrix[thread], CurrentProcessInfo);
rCurrentElement.CalculateMassMatrix(mMassMatrix[thread], rCurrentProcessInfo);

rCurrentElement.CalculateDampingMatrix(mDampingMatrix[thread], CurrentProcessInfo);
rCurrentElement.CalculateDampingMatrix(mDampingMatrix[thread], rCurrentProcessInfo);

this->AddDynamicsToLHS(LHS_Contribution, mMassMatrix[thread], mDampingMatrix[thread], CurrentProcessInfo);
this->AddDynamicsToLHS(rLHS_Contribution, mMassMatrix[thread], mDampingMatrix[thread], rCurrentProcessInfo);

rCurrentElement.EquationIdVector(EquationId, CurrentProcessInfo);
rCurrentElement.EquationIdVector(rEquationIds, rCurrentProcessInfo);

KRATOS_CATCH("")
}

protected:
void AddDynamicsToLHS(LocalSystemMatrixType& LHS_Contribution,
LocalSystemMatrixType& M,
LocalSystemMatrixType& C,
const ProcessInfo& CurrentProcessInfo)
void AddDynamicsToLHS(LocalSystemMatrixType& rLHS_Contribution,
LocalSystemMatrixType& rM,
LocalSystemMatrixType& rC,
const ProcessInfo& rCurrentProcessInfo)
{
KRATOS_TRY

// adding mass contribution
if (M.size1() != 0)
noalias(LHS_Contribution) +=
(1.0 / (this->GetBeta() * this->GetDeltaTime() * this->GetDeltaTime())) * M;
if (rM.size1() != 0)
noalias(rLHS_Contribution) +=
(1.0 / (this->GetBeta() * this->GetDeltaTime() * this->GetDeltaTime())) * rM;

// adding damping contribution
if (C.size1() != 0)
noalias(LHS_Contribution) += (this->GetGamma() / (this->GetBeta() * this->GetDeltaTime())) * C;
if (rC.size1() != 0)
noalias(rLHS_Contribution) += (this->GetGamma() / (this->GetBeta() * this->GetDeltaTime())) * rC;

KRATOS_CATCH("")
}

void AddDynamicsToRHS(Condition& rCurrentCondition,
LocalSystemVectorType& RHS_Contribution,
LocalSystemMatrixType& M,
LocalSystemMatrixType& C,
const ProcessInfo& CurrentProcessInfo)
LocalSystemVectorType& rRHS_Contribution,
LocalSystemMatrixType& rM,
LocalSystemMatrixType& rC,
const ProcessInfo& rCurrentProcessInfo)
{
KRATOS_TRY

int thread = OpenMPUtils::ThisThread();

// adding inertia contribution
if (M.size1() != 0) {
if (rM.size1() != 0) {
rCurrentCondition.GetSecondDerivativesVector(mAccelerationVector[thread], 0);
noalias(RHS_Contribution) -= prod(M, mAccelerationVector[thread]);
noalias(rRHS_Contribution) -= prod(rM, mAccelerationVector[thread]);
}

// adding damping contribution
if (C.size1() != 0) {
if (rC.size1() != 0) {
rCurrentCondition.GetFirstDerivativesVector(mVelocityVector[thread], 0);
noalias(RHS_Contribution) -= prod(C, mVelocityVector[thread]);
noalias(rRHS_Contribution) -= prod(rC, mVelocityVector[thread]);
}

KRATOS_CATCH("")
}

void AddDynamicsToRHS(Element& rCurrentElement,
LocalSystemVectorType& RHS_Contribution,
LocalSystemMatrixType& M,
LocalSystemMatrixType& C,
const ProcessInfo& CurrentProcessInfo)
LocalSystemVectorType& rRHS_Contribution,
LocalSystemMatrixType& rM,
LocalSystemMatrixType& rC,
const ProcessInfo& rCurrentProcessInfo)
{
KRATOS_TRY

int thread = OpenMPUtils::ThisThread();

// adding inertia contribution
if (M.size1() != 0) {
if (rM.size1() != 0) {
rCurrentElement.GetSecondDerivativesVector(mAccelerationVector[thread], 0);
noalias(RHS_Contribution) -= prod(M, mAccelerationVector[thread]);
noalias(rRHS_Contribution) -= prod(rM, mAccelerationVector[thread]);
}

// adding damping contribution
if (C.size1() != 0) {
if (rC.size1() != 0) {
rCurrentElement.GetFirstDerivativesVector(mVelocityVector[thread], 0);
noalias(RHS_Contribution) -= prod(C, mVelocityVector[thread]);
noalias(rRHS_Contribution) -= prod(rC, mVelocityVector[thread]);
}

KRATOS_CATCH("")
Expand Down
Loading
Loading