From 2c07bcd0ad00359725b84f7c2f3b83791b59ffad Mon Sep 17 00:00:00 2001 From: WPK4FEM <129858139+WPK4FEM@users.noreply.github.com> Date: Thu, 30 Nov 2023 17:56:07 +0100 Subject: [PATCH] [GeoMechanicsApplication] Implement SetValueOnIntegrationPoints for Vectors in U_Pw_small_strain_element. (#11863) * Implemented Vector overload SetValuesOnIntegrationPoints for UPwSmallStrainElement to enable K0 Procedure process for these elements. * Some reformatting, code smells and empty statements removed * Do not hide baseclass SetValuesOnIntegrationPoints. Simplified Info and PrintInfo. * avoid function hiding. Adressed some codesmells. Simplified Info and PrintInfo. * Default destructor, as indicated on SonarCloud. Some more tassel removed. * Review comments of Anne taken up * More review comments by Anne taken up. * Still more of the review comments by Anne taken up. --- .../custom_elements/U_Pw_base_element.hpp | 96 +- .../U_Pw_small_strain_element.cpp | 1286 +++++++---------- .../U_Pw_small_strain_element.hpp | 241 ++- .../mesh_stage1.mdpa | 2 +- 4 files changed, 707 insertions(+), 918 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.hpp index b5a3dce08779..00555e1dec2d 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.hpp @@ -11,8 +11,7 @@ // Vahid Galavi // -#if !defined(KRATOS_GEO_U_PW_ELEMENT_H_INCLUDED ) -#define KRATOS_GEO_U_PW_ELEMENT_H_INCLUDED +#pragma once // Project includes #include "containers/array_1d.h" @@ -37,16 +36,11 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element { public: - - /// The definition of the sizetype using SizeType = std::size_t; KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION( UPwBaseElement ); -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - - /// Default Constructor - UPwBaseElement(IndexType NewId = 0) : Element( NewId ) {} + explicit UPwBaseElement(IndexType NewId = 0) : Element( NewId ) {} /// Constructor using an array of nodes UPwBaseElement(IndexType NewId, const NodesArrayType& ThisNodes) : Element(NewId, ThisNodes) {} @@ -63,10 +57,11 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element mThisIntegrationMethod = this->GetIntegrationMethod(); } - /// Destructor - virtual ~UPwBaseElement() {} - -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + ~UPwBaseElement() override = default; + UPwBaseElement(const UPwBaseElement&) = delete; + UPwBaseElement& operator=(const UPwBaseElement&) = delete; + UPwBaseElement(UPwBaseElement&&) noexcept = delete; + UPwBaseElement& operator=(UPwBaseElement&&) noexcept = delete; Element::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes, @@ -86,26 +81,25 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element void ResetConstitutiveLaw() override; GeometryData::IntegrationMethod GetIntegrationMethod() const override; -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - void CalculateLocalSystem( MatrixType& rLeftHandSideMatrix, - VectorType& rRightHandSideVector, - const ProcessInfo& rCurrentProcessInfo ) override; + void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, + VectorType& rRightHandSideVector, + const ProcessInfo& rCurrentProcessInfo) override; - void CalculateLeftHandSide( MatrixType& rLeftHandSideMatrix, - const ProcessInfo& rCurrentProcessInfo ) override; + void CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, + const ProcessInfo& rCurrentProcessInfo) override; - void CalculateRightHandSide( VectorType& rRightHandSideVector, - const ProcessInfo& rCurrentProcessInfo ) override; + void CalculateRightHandSide(VectorType& rRightHandSideVector, + const ProcessInfo& rCurrentProcessInfo) override; - void EquationIdVector( EquationIdVectorType& rResult, - const ProcessInfo& rCurrentProcessInfo ) const override; + void EquationIdVector(EquationIdVectorType& rResult, + const ProcessInfo& rCurrentProcessInfo) const override; - void CalculateMassMatrix( MatrixType& rMassMatrix, - const ProcessInfo& rCurrentProcessInfo ) override; + void CalculateMassMatrix(MatrixType& rMassMatrix, + const ProcessInfo& rCurrentProcessInfo) override; - void CalculateDampingMatrix( MatrixType& rDampingMatrix, - const ProcessInfo& rCurrentProcessInfo ) override; + void CalculateDampingMatrix(MatrixType& rDampingMatrix, + const ProcessInfo& rCurrentProcessInfo) override; void GetValuesVector(Vector& rValues, int Step = 0) const override; @@ -113,27 +107,27 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element void GetSecondDerivativesVector(Vector& rValues, int Step = 0) const override; -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - void CalculateOnIntegrationPoints( const Variable& rVariable, + void CalculateOnIntegrationPoints(const Variable& rVariable, std::vector& rValues, const ProcessInfo& rCurrentProcessInfo ) override; void CalculateOnIntegrationPoints(const Variable>& rVariable, - std::vector>& rValues, - const ProcessInfo& rCurrentProcessInfo) override; + std::vector>& rValues, + const ProcessInfo& rCurrentProcessInfo) override; void CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector& rValues, - const ProcessInfo& rCurrentProcessInfo) override; + std::vector& rValues, + const ProcessInfo& rCurrentProcessInfo) override; void CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector& rValues, - const ProcessInfo& rCurrentProcessInfo) override; + std::vector& rValues, + const ProcessInfo& rCurrentProcessInfo) override; void CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector& rValues, - const ProcessInfo& rCurrentProcessInfo) override; + std::vector& rValues, + const ProcessInfo& rCurrentProcessInfo) override; + using Element::CalculateOnIntegrationPoints; void SetValuesOnIntegrationPoints(const Variable& rVariable, const std::vector& rValues, @@ -147,22 +141,18 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element const std::vector& rValues, const ProcessInfo& rCurrentProcessInfo) override; - // Turn back information as a string. + using Element::SetValuesOnIntegrationPoints; + std::string Info() const override { - std::stringstream buffer; - buffer << "U-Pw Base class Element #" << this->Id() << "\nConstitutive law: " << mConstitutiveLawVector[0]->Info(); - return buffer.str(); + return "U-Pw Base class Element #" + std::to_string(Id()) + "\nConstitutive law: " + mConstitutiveLawVector[0]->Info(); } - // Print information about this object. void PrintInfo(std::ostream& rOStream) const override { - rOStream << "U-Pw Base class Element #" << this->Id() << "\nConstitutive law: " << mConstitutiveLawVector[0]->Info(); + rOStream << Info(); } -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - protected: /// Member Variables @@ -175,16 +165,14 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element std::vector mStateVariablesFinalized; bool mIsInitialised = false; -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - virtual void CalculateMaterialStiffnessMatrix( MatrixType& rStiffnessMatrix, const ProcessInfo& CurrentProcessInfo ); virtual void CalculateAll( MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, const ProcessInfo& CurrentProcessInfo, - const bool CalculateStiffnessMatrixFlag, - const bool CalculateResidualVectorFlag); + bool CalculateStiffnessMatrixFlag, + bool CalculateResidualVectorFlag); virtual double CalculateIntegrationCoefficient( const GeometryType::IntegrationPointsArrayType& IntegrationPoints, unsigned int PointNumber, @@ -233,12 +221,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element virtual unsigned int GetNumberOfDOF() const; -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - private: - /// Serialization - friend class Serializer; void save(Serializer& rSerializer) const override @@ -251,15 +235,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element KRATOS_SERIALIZE_LOAD_BASE_CLASS( rSerializer, Element ) } - /// Assignment operator. - UPwBaseElement & operator=(UPwBaseElement const& rOther); - - /// Copy constructor. - UPwBaseElement(UPwBaseElement const& rOther); - }; // Class UPwBaseElement } // namespace Kratos - -#endif // KRATOS_GEO_U_PW_ELEMENT_H_INCLUDED defined diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp index 28158aa8b7bd..13cc88091b78 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -17,8 +17,7 @@ namespace Kratos { -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > +template < unsigned int TDim, unsigned int TNumNodes > Element::Pointer UPwSmallStrainElement:: Create(IndexType NewId, NodesArrayType const& ThisNodes, @@ -27,8 +26,7 @@ Element::Pointer UPwSmallStrainElement:: return Element::Pointer( new UPwSmallStrainElement( NewId, this->GetGeometry().Create( ThisNodes ), pProperties ) ); } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > +template < unsigned int TDim, unsigned int TNumNodes > Element::Pointer UPwSmallStrainElement:: Create(IndexType NewId, GeometryType::Pointer pGeom, @@ -37,10 +35,8 @@ Element::Pointer UPwSmallStrainElement:: return Element::Pointer( new UPwSmallStrainElement( NewId, pGeom, pProperties ) ); } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -int UPwSmallStrainElement:: - Check( const ProcessInfo& rCurrentProcessInfo ) const +template +int UPwSmallStrainElement::Check(const ProcessInfo& rCurrentProcessInfo) const { KRATOS_TRY @@ -52,44 +48,40 @@ int UPwSmallStrainElement:: const PropertiesType& rProp = this->GetProperties(); const GeometryType& rGeom = this->GetGeometry(); - if (rGeom.DomainSize() < 1.0e-15) - KRATOS_ERROR << "DomainSize < 1.0e-15 for the element " << this->Id() << std::endl; + KRATOS_ERROR_IF(rGeom.DomainSize() < 1.0e-15) << "DomainSize < 1.0e-15 for the element " << this->Id() << std::endl; // Verify specific properties - if ( rProp.Has( IGNORE_UNDRAINED ) == false) - KRATOS_ERROR << "IGNORE_UNDRAINED does not exist in the parameter list" << this->Id() << std::endl; + KRATOS_ERROR_IF_NOT(rProp.Has(IGNORE_UNDRAINED)) + << "IGNORE_UNDRAINED does not exist in the parameter list" << this->Id() << std::endl; - bool IgnoreUndrained = rProp[IGNORE_UNDRAINED]; - - if (!IgnoreUndrained) { - if ( rProp.Has( BULK_MODULUS_FLUID ) == false || rProp[BULK_MODULUS_FLUID] < 0.0 ) - KRATOS_ERROR << "BULK_MODULUS_FLUID has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; + if (!rProp[IGNORE_UNDRAINED]) { + KRATOS_ERROR_IF(!rProp.Has( BULK_MODULUS_FLUID ) || rProp[BULK_MODULUS_FLUID] < 0.0) + << "BULK_MODULUS_FLUID has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; - if ( rProp.Has( DYNAMIC_VISCOSITY ) == false || rProp[DYNAMIC_VISCOSITY] < 0.0 ) - KRATOS_ERROR << "DYNAMIC_VISCOSITY has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; + KRATOS_ERROR_IF(!rProp.Has( DYNAMIC_VISCOSITY ) || rProp[DYNAMIC_VISCOSITY] < 0.0) + << "DYNAMIC_VISCOSITY has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; - if ( rProp.Has( PERMEABILITY_XX ) == false || rProp[PERMEABILITY_XX] < 0.0 ) - KRATOS_ERROR << "PERMEABILITY_XX has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; + KRATOS_ERROR_IF(!rProp.Has( PERMEABILITY_XX ) || rProp[PERMEABILITY_XX] < 0.0) + << "PERMEABILITY_XX has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; - if ( rProp.Has( PERMEABILITY_YY ) == false || rProp[PERMEABILITY_YY] < 0.0 ) - KRATOS_ERROR << "PERMEABILITY_YY has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; + KRATOS_ERROR_IF(!rProp.Has( PERMEABILITY_YY ) || rProp[PERMEABILITY_YY] < 0.0) + << "PERMEABILITY_YY has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; - if ( rProp.Has( PERMEABILITY_XY ) == false || rProp[PERMEABILITY_XY] < 0.0 ) - KRATOS_ERROR << "PERMEABILITY_XY has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; + KRATOS_ERROR_IF(!rProp.Has( PERMEABILITY_XY ) || rProp[PERMEABILITY_XY] < 0.0) + << "PERMEABILITY_XY has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; if constexpr (TDim > 2) { - if ( rProp.Has( PERMEABILITY_ZZ ) == false || rProp[PERMEABILITY_ZZ] < 0.0 ) - KRATOS_ERROR << "PERMEABILITY_ZZ has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; + KRATOS_ERROR_IF(!rProp.Has( PERMEABILITY_ZZ ) || rProp[PERMEABILITY_ZZ] < 0.0) + << "PERMEABILITY_ZZ has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; - if ( rProp.Has( PERMEABILITY_YZ ) == false || rProp[PERMEABILITY_YZ] < 0.0 ) - KRATOS_ERROR << "PERMEABILITY_YZ has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; + KRATOS_ERROR_IF(!rProp.Has( PERMEABILITY_YZ ) || rProp[PERMEABILITY_YZ] < 0.0) + << "PERMEABILITY_YZ has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; - if ( rProp.Has( PERMEABILITY_ZX ) == false || rProp[PERMEABILITY_ZX] < 0.0 ) - KRATOS_ERROR << "PERMEABILITY_ZX has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; + KRATOS_ERROR_IF(!rProp.Has( PERMEABILITY_ZX ) || rProp[PERMEABILITY_ZX] < 0.0) + << "PERMEABILITY_ZX has Key zero, is not defined or has an invalid value at element" << this->Id() << std::endl; } } - // Verify that the constitutive law exists KRATOS_ERROR_IF_NOT(this->GetProperties().Has( CONSTITUTIVE_LAW )) << "Constitutive law not provided for property " << this->GetProperties().Id() << std::endl; @@ -115,53 +107,52 @@ int UPwSmallStrainElement:: } // Check constitutive law - if ( mConstitutiveLawVector.size() > 0 ) { - return mConstitutiveLawVector[0]->Check( rProp, rGeom, rCurrentProcessInfo ); + if (mConstitutiveLawVector.size() > 0) { + return mConstitutiveLawVector[0]->Check(rProp, rGeom, rCurrentProcessInfo); } - // Check constitutive law - if ( mRetentionLawVector.size() > 0 ) { - return mRetentionLawVector[0]->Check( rProp, rCurrentProcessInfo ); + // Check retention law + if (mRetentionLawVector.size() > 0) { + return mRetentionLawVector[0]->Check(rProp, rCurrentProcessInfo); } return ierr; - KRATOS_CATCH( "" ); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement::InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) +template +void UPwSmallStrainElement::InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) { - KRATOS_TRY; + KRATOS_TRY if (!mIsInitialised) this->Initialize(rCurrentProcessInfo); - //Defining necessary variables + // Defining necessary variables const GeometryType& rGeom = this->GetGeometry(); - const IndexType NumGPoints = rGeom.IntegrationPointsNumber( mThisIntegrationMethod ); + const IndexType NumGPoints = rGeom.IntegrationPointsNumber(mThisIntegrationMethod); ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom, this->GetProperties(), rCurrentProcessInfo); ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); ConstitutiveParameters.Set(ConstitutiveLaw::INITIALIZE_MATERIAL_RESPONSE); //Note: this is for nonlocal damage - //Element variables + // Element variables ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // create general parameters of retention law + // Create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); - //Loop over integration points + // Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - //Compute Np, GradNpT, B and StrainVector + // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - //Compute infinitessimal strain + // Compute infinitesimal strain this->CalculateStrain(Variables, GPoint); - //set gauss points variables to constitutivelaw parameters + // Set Gauss points variables to constitutive law parameters this->SetConstitutiveParameters(Variables, ConstitutiveParameters); // Initialize constitutive law @@ -173,60 +164,55 @@ void UPwSmallStrainElement::InitializeSolutionStep(const Process mRetentionLawVector[GPoint]->InitializeSolutionStep(RetentionParameters); } - // reset hydraulic discharge + // Reset hydraulic discharge this->ResetHydraulicDischarge(); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - ResetHydraulicDischarge() +template +void UPwSmallStrainElement::ResetHydraulicDischarge() { - KRATOS_TRY; + KRATOS_TRY - // reset hydraulic discharge + // Reset hydraulic discharge GeometryType& rGeom = this->GetGeometry(); for (unsigned int i=0; i -void UPwSmallStrainElement:: - CalculateHydraulicDischarge(const ProcessInfo& rCurrentProcessInfo) +template +void UPwSmallStrainElement::CalculateHydraulicDischarge(const ProcessInfo& rCurrentProcessInfo) { - KRATOS_TRY; + KRATOS_TRY std::vector> FluidFlux; - this->CalculateOnIntegrationPoints( FLUID_FLUX_VECTOR, - FluidFlux, - rCurrentProcessInfo ); + this->CalculateOnIntegrationPoints(FLUID_FLUX_VECTOR, + FluidFlux, + rCurrentProcessInfo); const GeometryType& rGeom = this->GetGeometry(); - const IndexType NumGPoints = rGeom.IntegrationPointsNumber( mThisIntegrationMethod ); - const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints( mThisIntegrationMethod ); + const IndexType NumGPoints = rGeom.IntegrationPointsNumber(mThisIntegrationMethod); + const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(mThisIntegrationMethod); ElementVariables Variables; - // gradient of shape functions and determinant of Jacobian - + // Gradient of shape functions and determinant of Jacobian Variables.GradNpTInitialConfiguration.resize(TNumNodes, TDim, false); Variables.GradNpT.resize(TNumNodes, TDim, false); - (Variables.detJContainer).resize(NumGPoints, false); - rGeom.ShapeFunctionsIntegrationPointsGradients( Variables.DN_DXContainer, + Variables.detJContainer.resize(NumGPoints, false); + rGeom.ShapeFunctionsIntegrationPointsGradients(Variables.DN_DXContainer, Variables.detJContainer, - mThisIntegrationMethod ); + mThisIntegrationMethod); - //Loop over integration points + // Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { noalias(Variables.GradNpT) = Variables.DN_DXContainer[GPoint]; Variables.detJ = Variables.detJContainer[GPoint]; - //Compute weighting coefficient for integration + // Compute weighting coefficient for integration Variables.IntegrationCoefficient = this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, @@ -244,39 +230,37 @@ void UPwSmallStrainElement:: } } - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - InitializeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) +template +void UPwSmallStrainElement::InitializeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) { - KRATOS_TRY; + KRATOS_TRY - //Defining necessary variables + // Defining necessary variables const GeometryType& rGeom = this->GetGeometry(); - const IndexType NumGPoints = rGeom.IntegrationPointsNumber( mThisIntegrationMethod ); + const IndexType NumGPoints = rGeom.IntegrationPointsNumber(mThisIntegrationMethod); - //Create constitutive law parameters: + // Create constitutive law parameters: ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom,this->GetProperties(),rCurrentProcessInfo); ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_STRESS); ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); ConstitutiveParameters.Set(ConstitutiveLaw::INITIALIZE_MATERIAL_RESPONSE); //Note: this is for nonlocal damage - //Element variables + // Element variables ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - //Loop over integration points + // Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { this->CalculateKinematics(Variables, GPoint); - //Compute infinitessimal strain + // Compute infinitesimal strain this->CalculateStrain(Variables, GPoint); - //set gauss points variables to constitutivelaw parameters + // Set gauss points variables to constitutive law parameters this->SetConstitutiveParameters(Variables, ConstitutiveParameters); // Compute Stress @@ -284,63 +268,59 @@ void UPwSmallStrainElement:: mConstitutiveLawVector[GPoint]->CalculateMaterialResponseCauchy(ConstitutiveParameters); } - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------- - -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement::FinalizeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) +template +void UPwSmallStrainElement::FinalizeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) { - KRATOS_TRY; + KRATOS_TRY this->InitializeNonLinearIteration(rCurrentProcessInfo); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo ) +template +void UPwSmallStrainElement::FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) { KRATOS_TRY this->CalculateHydraulicDischarge(rCurrentProcessInfo); - //Defining necessary variables + // Defining necessary variables const GeometryType& rGeom = this->GetGeometry(); const IndexType NumGPoints = rGeom.IntegrationPointsNumber( mThisIntegrationMethod ); ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom,this->GetProperties(),rCurrentProcessInfo); ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); - //Element variables + // Element variables ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // create general parameters of retention law + // Create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); Matrix StressContainer(NumGPoints, mStressVector[0].size()); - //Loop over integration points + // Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { this->CalculateKinematics(Variables, GPoint); - //Compute infinitesimal strain + // Compute infinitesimal strain this->CalculateStrain(Variables, GPoint); - //set Gauss points variables to constitutive law parameters + // Set Gauss points variables to constitutive law parameters this->SetConstitutiveParameters(Variables, ConstitutiveParameters); - //compute constitutive tensor and/or stresses + // Compute constitutive tensor and/or stresses noalias(Variables.StressVector) = mStressVector[GPoint]; ConstitutiveParameters.SetStressVector(Variables.StressVector); mConstitutiveLawVector[GPoint]->FinalizeMaterialResponseCauchy(ConstitutiveParameters); mStateVariablesFinalized[GPoint] = - mConstitutiveLawVector[GPoint]->GetValue( STATE_VARIABLES, - mStateVariablesFinalized[GPoint] ); + mConstitutiveLawVector[GPoint]->GetValue(STATE_VARIABLES, + mStateVariablesFinalized[GPoint]); // retention law mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters); @@ -351,14 +331,13 @@ void UPwSmallStrainElement:: if (rCurrentProcessInfo[NODAL_SMOOTHING]) this->ExtrapolateGPValues(StressContainer); - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement::SaveGPStress(Matrix& rStressContainer, - const Vector& StressVector, - unsigned int GPoint) +template +void UPwSmallStrainElement::SaveGPStress(Matrix& rStressContainer, + const Vector& StressVector, + unsigned int GPoint) { KRATOS_TRY @@ -376,11 +355,10 @@ void UPwSmallStrainElement::SaveGPStress(Matrix& rStressContaine * S1-0 = S[1] at GP 0 */ - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > +template void UPwSmallStrainElement::ExtrapolateGPValues(const Matrix& StressContainer) { KRATOS_TRY @@ -389,7 +367,7 @@ void UPwSmallStrainElement::ExtrapolateGPValues(const Matrix& S for ( unsigned int iNode = 0; iNode < TNumNodes; ++iNode ) { DamageContainer[iNode] = 0.0; - DamageContainer[iNode] = mConstitutiveLawVector[iNode]->GetValue( DAMAGE_VARIABLE, DamageContainer[iNode] ); + DamageContainer[iNode] = mConstitutiveLawVector[iNode]->GetValue(DAMAGE_VARIABLE, DamageContainer[iNode]); } GeometryType& rGeom = this->GetGeometry(); @@ -432,77 +410,91 @@ void UPwSmallStrainElement::ExtrapolateGPValues(const Matrix& S rGeom[i].UnSetLock(); } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement::CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector& rOutput, - const ProcessInfo& rCurrentProcessInfo) +template +void UPwSmallStrainElement::SetValuesOnIntegrationPoints(const Variable& rVariable, + const std::vector& rValues, + const ProcessInfo& rCurrentProcessInfo) { KRATOS_TRY - //Defining necessary variables + + if (rVariable == CAUCHY_STRESS_VECTOR) { + KRATOS_ERROR_IF (rValues.size() != mStressVector.size()) << + "Unexpected number of values for UPwSmallStrainElement::SetValuesOnIntegrationPoints" << std::endl; + std::copy(rValues.begin(), rValues.end(), mStressVector.begin()); + } else { + KRATOS_ERROR_IF (rValues.size() < mConstitutiveLawVector.size()) << + "Insufficient number of values for UPwSmallStrainElement::SetValuesOnIntegrationPoints" << std::endl; + for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { + mConstitutiveLawVector[GPoint]->SetValue(rVariable, rValues[GPoint], rCurrentProcessInfo); + } + } + + KRATOS_CATCH("") +} + +template +void UPwSmallStrainElement::CalculateOnIntegrationPoints(const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) +{ + KRATOS_TRY + // Defining necessary variables const GeometryType& rGeom = this->GetGeometry(); - const IndexType NumGPoints = rGeom.IntegrationPointsNumber( mThisIntegrationMethod ); + const IndexType NumGPoints = rGeom.IntegrationPointsNumber(mThisIntegrationMethod); auto& r_prop = this->GetProperties(); if ( rOutput.size() != NumGPoints ) rOutput.resize(NumGPoints); if (rVariable == VON_MISES_STRESS) { - //Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { StressStrainUtilities EquivalentStress; rOutput[GPoint] = EquivalentStress.CalculateVonMisesStress(mStressVector[GPoint]); } } else if (rVariable == MEAN_EFFECTIVE_STRESS ) { - //Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { StressStrainUtilities EquivalentStress; rOutput[GPoint] = EquivalentStress.CalculateMeanStress(mStressVector[GPoint]); } } else if (rVariable == MEAN_STRESS ) { std::vector StressVector; - CalculateOnIntegrationPoints( TOTAL_STRESS_VECTOR, StressVector, rCurrentProcessInfo ); + CalculateOnIntegrationPoints(TOTAL_STRESS_VECTOR, StressVector, rCurrentProcessInfo); - //loop integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { StressStrainUtilities EquivalentStress; rOutput[GPoint] = EquivalentStress.CalculateMeanStress(StressVector[GPoint]); } } else if (rVariable == ENGINEERING_VON_MISES_STRAIN ) { std::vector StrainVector; - CalculateOnIntegrationPoints( ENGINEERING_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo ); + CalculateOnIntegrationPoints(ENGINEERING_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo); - //loop integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { StressStrainUtilities EquivalentStress; rOutput[GPoint] = EquivalentStress.CalculateVonMisesStrain(StrainVector[GPoint]); } } else if (rVariable == ENGINEERING_VOLUMETRIC_STRAIN ) { std::vector StrainVector; - CalculateOnIntegrationPoints( ENGINEERING_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo ); + CalculateOnIntegrationPoints(ENGINEERING_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo); - //loop integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { StressStrainUtilities EquivalentStress; rOutput[GPoint] = EquivalentStress.CalculateTrace(StrainVector[GPoint]); } } else if (rVariable == GREEN_LAGRANGE_VON_MISES_STRAIN ) { std::vector StrainVector; - CalculateOnIntegrationPoints( GREEN_LAGRANGE_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo ); + CalculateOnIntegrationPoints(GREEN_LAGRANGE_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo); - //loop integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { StressStrainUtilities EquivalentStress; rOutput[GPoint] = EquivalentStress.CalculateVonMisesStrain(StrainVector[GPoint]); } } else if (rVariable == GREEN_LAGRANGE_VOLUMETRIC_STRAIN ) { std::vector StrainVector; - CalculateOnIntegrationPoints( GREEN_LAGRANGE_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo ); + CalculateOnIntegrationPoints(GREEN_LAGRANGE_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo); - //loop integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { StressStrainUtilities EquivalentStress; rOutput[GPoint] = EquivalentStress.CalculateTrace(StrainVector[GPoint]); @@ -513,7 +505,6 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const V rVariable == DERIVATIVE_OF_SATURATION || rVariable == RELATIVE_PERMEABILITY ) { - //Element variables ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); @@ -521,29 +512,27 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const V // create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); - //Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { //Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint ); + this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); if (rVariable == DEGREE_OF_SATURATION) rOutput[GPoint] = Variables.DegreeOfSaturation; if (rVariable == EFFECTIVE_SATURATION) rOutput[GPoint] = Variables.EffectiveSaturation; if (rVariable == BISHOP_COEFFICIENT) rOutput[GPoint] = Variables.BishopCoefficient; if (rVariable == DERIVATIVE_OF_SATURATION) rOutput[GPoint] = Variables.DerivativeOfSaturation; - if (rVariable == RELATIVE_PERMEABILITY ) rOutput[GPoint] = Variables.RelativePermeability; + if (rVariable == RELATIVE_PERMEABILITY) rOutput[GPoint] = Variables.RelativePermeability; } } else if (rVariable == HYDRAULIC_HEAD) { const PropertiesType& rProp = this->GetProperties(); - //Defining the shape functions, the Jacobian and the shape functions local gradients Containers - const Matrix& NContainer = rGeom.ShapeFunctionsValues( mThisIntegrationMethod ); + // Defining the shape functions, the Jacobian and the shape functions local gradients containers + const Matrix& NContainer = rGeom.ShapeFunctionsValues(mThisIntegrationMethod); const auto NodalHydraulicHead = GeoElementUtilities::CalculateNodalHydraulicHeadFromWaterPressures(rGeom, rProp); - //Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { double HydraulicHead = 0.0; for (unsigned int node = 0; node < TNumNodes; ++node) @@ -554,40 +543,27 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const V } else if (rVariable == CONFINED_STIFFNESS || rVariable == SHEAR_STIFFNESS) { size_t variable_index; - if (rVariable == CONFINED_STIFFNESS) - { - if (TDim == 2) - { + if (rVariable == CONFINED_STIFFNESS) { + if (TDim == 2) { variable_index = INDEX_2D_PLANE_STRAIN_XX; - } - else if (TDim == 3) - { + } else if (TDim == 3) { variable_index = INDEX_3D_XX; - } - else - { + } else { KRATOS_ERROR << "CONFINED_STIFFNESS can not be retrieved for dim " << TDim << " in element: " << this->Id() << std::endl; } - } - else if (rVariable == SHEAR_STIFFNESS) - { - if (TDim == 2) - { + } else if (rVariable == SHEAR_STIFFNESS) { + if (TDim == 2) { variable_index = INDEX_2D_PLANE_STRAIN_XY; - } - else if (TDim == 3) - { + } else if (TDim == 3) { variable_index = INDEX_3D_XZ; - } - else - { + } else { KRATOS_ERROR << "SHEAR_STIFFNESS can not be retrieved for dim " << TDim << " in element: " << this->Id() << std::endl; } } ElementVariables Variables; this->InitializeElementVariables(Variables, - rCurrentProcessInfo); + rCurrentProcessInfo); if (rOutput.size() != mConstitutiveLawVector.size()) rOutput.resize(mConstitutiveLawVector.size()); @@ -599,7 +575,6 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const V ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR); - //Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); @@ -620,7 +595,7 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const V } else if (r_prop.Has(rVariable)) { - // map initial material property to gauss points, as required for the output + // Map initial material property to gauss points, as required for the output rOutput.clear(); std::fill_n(std::back_inserter(rOutput), mConstitutiveLawVector.size(), r_prop.GetValue(rVariable)); } @@ -629,24 +604,22 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const V rOutput.resize(mConstitutiveLawVector.size()); for ( unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i ) { - rOutput[i] = mConstitutiveLawVector[i]->GetValue( rVariable, rOutput[i] ); + rOutput[i] = mConstitutiveLawVector[i]->GetValue(rVariable, rOutput[i]); } } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateOnIntegrationPoints( const Variable>& rVariable, - std::vector>& rOutput, - const ProcessInfo& rCurrentProcessInfo ) +template +void UPwSmallStrainElement::CalculateOnIntegrationPoints(const Variable>& rVariable, + std::vector>& rOutput, + const ProcessInfo& rCurrentProcessInfo) { KRATOS_TRY const GeometryType& rGeom = this->GetGeometry(); - const IndexType NumGPoints = rGeom.IntegrationPointsNumber( mThisIntegrationMethod ); + const IndexType NumGPoints = rGeom.IntegrationPointsNumber(mThisIntegrationMethod); if ( rOutput.size() != NumGPoints ) rOutput.resize(NumGPoints); @@ -659,29 +632,27 @@ void UPwSmallStrainElement:: array_1d GradPressureTerm; array_1d FluidFlux; - // create general parameters of retention law + // Create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); - //Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { - //Compute Np, GradNpT, B and StrainVector + // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - //Compute strain, needed for update permeability + // Compute strain, needed for update permeability this->CalculateStrain(Variables, GPoint); this->CalculatePermeabilityUpdateFactor(Variables); - GeoElementUtilities:: - InterpolateVariableWithComponents( Variables.BodyAcceleration, - Variables.NContainer, - Variables.VolumeAcceleration, - GPoint ); + GeoElementUtilities::InterpolateVariableWithComponents(Variables.BodyAcceleration, + Variables.NContainer, + Variables.VolumeAcceleration, + GPoint); Variables.FluidPressure = CalculateFluidPressure(Variables); RetentionParameters.SetFluidPressure(Variables.FluidPressure); Variables.RelativePermeability = mRetentionLawVector[GPoint]-> - CalculateRelativePermeability(RetentionParameters); + CalculateRelativePermeability(RetentionParameters); noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector); @@ -703,32 +674,29 @@ void UPwSmallStrainElement:: for ( unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i ) { noalias(rOutput[i]) = ZeroVector(3); - rOutput[i] = mConstitutiveLawVector[i]->GetValue( rVariable, rOutput[i] ); + rOutput[i] = mConstitutiveLawVector[i]->GetValue(rVariable, rOutput[i]); } } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector& rOutput, - const ProcessInfo& rCurrentProcessInfo ) +template +void UPwSmallStrainElement::CalculateOnIntegrationPoints(const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) { KRATOS_TRY //Defining necessary variables const GeometryType& rGeom = this->GetGeometry(); - const IndexType NumGPoints = rGeom.IntegrationPointsNumber( mThisIntegrationMethod ); + const IndexType NumGPoints = rGeom.IntegrationPointsNumber(mThisIntegrationMethod); const PropertiesType& rProp = this->GetProperties(); if ( rOutput.size() != NumGPoints ) rOutput.resize(NumGPoints); if (rVariable == CAUCHY_STRESS_VECTOR) { - //Loop over integration points for ( unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint ) { if ( rOutput[GPoint].size() != mStressVector[GPoint].size() ) rOutput[GPoint].resize( mStressVector[GPoint].size(), false ); @@ -736,18 +704,17 @@ void UPwSmallStrainElement:: rOutput[GPoint] = mStressVector[GPoint]; } } else if (rVariable == TOTAL_STRESS_VECTOR) { - //Defining necessary variables - + // Defining necessary variables ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom,rProp,rCurrentProcessInfo); ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR); ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); - //Element variables + // Element variables ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // create general parameters of retention law + // Create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); Vector VoigtVector(mStressVector[0].size()); @@ -759,18 +726,18 @@ void UPwSmallStrainElement:: Vector TotalStressVector(mStressVector[0].size()); - //Loop over integration points + // Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { - //Compute Np, GradNpT, B and StrainVector + // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - //Compute infinitesimal strain + // Compute infinitesimal strain this->CalculateStrain(Variables, GPoint); - //set Gauss points variables to constitutive law parameters + // Set Gauss points variables to constitutive law parameters this->SetConstitutiveParameters(Variables, ConstitutiveParameters); - //compute constitutive tensor and/or stresses + // Compute constitutive tensor and/or stresses noalias(Variables.StressVector) = mStressVector[GPoint]; ConstitutiveParameters.SetStressVector(Variables.StressVector); mConstitutiveLawVector[GPoint]->CalculateMaterialResponseCauchy(ConstitutiveParameters); @@ -786,18 +753,17 @@ void UPwSmallStrainElement:: * Variables.FluidPressure * VoigtVector; - if ( rOutput[GPoint].size() != TotalStressVector.size() ) - rOutput[GPoint].resize( TotalStressVector.size(), false ); + if (rOutput[GPoint].size() != TotalStressVector.size()) + rOutput[GPoint].resize(TotalStressVector.size(), false); rOutput[GPoint] = TotalStressVector; } } else if (rVariable == ENGINEERING_STRAIN_VECTOR) { - //Definition of variables + // Definition of variables ElementVariables Variables; this->InitializeElementVariables(Variables,rCurrentProcessInfo); - //Loop over integration points - for ( unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint ) { + for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { noalias(Variables.Np) = row(Variables.NContainer, GPoint); Matrix J0,InvJ0; @@ -810,66 +776,56 @@ void UPwSmallStrainElement:: // Calculating operator B this->CalculateBMatrix( Variables.B, Variables.GradNpTInitialConfiguration, Variables.Np); - //Compute infinitessimal strain - this->CalculateCauchyStrain( Variables ); + // Compute infinitesimal strain + this->CalculateCauchyStrain(Variables); - if ( rOutput[GPoint].size() != Variables.StrainVector.size() ) - rOutput[GPoint].resize( Variables.StrainVector.size(), false ); + if (rOutput[GPoint].size() != Variables.StrainVector.size()) + rOutput[GPoint].resize(Variables.StrainVector.size(), false); rOutput[GPoint] = Variables.StrainVector; } } else if (rVariable == GREEN_LAGRANGE_STRAIN_VECTOR) { - //Definition of variables ElementVariables Variables; this->InitializeElementVariables(Variables,rCurrentProcessInfo); - //Loop over integration points for ( unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint ) { - //compute element kinematics (Np, gradNpT, |J|, B, strains) + // Compute element kinematics (Np, gradNpT, |J|, B, strains) this->CalculateKinematics(Variables,GPoint); - //Compute strain this->CalculateStrain(Variables, GPoint); - if ( rOutput[GPoint].size() != Variables.StrainVector.size() ) - rOutput[GPoint].resize( Variables.StrainVector.size(), false ); + if (rOutput[GPoint].size() != Variables.StrainVector.size()) + rOutput[GPoint].resize(Variables.StrainVector.size(), false); rOutput[GPoint] = Variables.StrainVector; } - } - else if (rProp.Has(rVariable)) - { - // map initial material property to gauss points, as required for the output + } else if (rProp.Has(rVariable)){ + // Map initial material property to Gauss points, as required for the output rOutput.clear(); std::fill_n(std::back_inserter(rOutput), mConstitutiveLawVector.size(), rProp.GetValue(rVariable)); - } - else { - for ( unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i ) + } else { + for (unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i) rOutput[i] = mConstitutiveLawVector[i]->GetValue( rVariable , rOutput[i] ); } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateOnIntegrationPoints(const Variable& rVariable, - std::vector& rOutput, - const ProcessInfo& rCurrentProcessInfo ) +template +void UPwSmallStrainElement::CalculateOnIntegrationPoints(const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo) { KRATOS_TRY - //Defining necessary variables const GeometryType& rGeom = this->GetGeometry(); const IndexType NumGPoints = rGeom.IntegrationPointsNumber( mThisIntegrationMethod ); - if ( rOutput.size() != NumGPoints ) + if (rOutput.size() != NumGPoints) rOutput.resize(NumGPoints); if (rVariable == CAUCHY_STRESS_TENSOR) { - //Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { rOutput[GPoint].resize(StressTensorSize, StressTensorSize, false ); rOutput[GPoint] = MathUtils::StressVectorToTensor(mStressVector[GPoint]); @@ -879,7 +835,6 @@ void UPwSmallStrainElement:: this->CalculateOnIntegrationPoints(TOTAL_STRESS_VECTOR, StressVector, rCurrentProcessInfo); - //loop integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { if (rOutput[GPoint].size2() != TDim) rOutput[GPoint].resize(TDim, TDim, false); @@ -890,171 +845,154 @@ void UPwSmallStrainElement:: } else if (rVariable == ENGINEERING_STRAIN_TENSOR) { std::vector StrainVector; - CalculateOnIntegrationPoints( ENGINEERING_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo ); + CalculateOnIntegrationPoints(ENGINEERING_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo); - //loop integration points - for ( unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint ) { - if ( rOutput[GPoint].size2() != TDim ) - rOutput[GPoint].resize( TDim, TDim, false ); + for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { + if (rOutput[GPoint].size2() != TDim) + rOutput[GPoint].resize(TDim, TDim, false); rOutput[GPoint] = MathUtils::StrainVectorToTensor(StrainVector[GPoint]); } } else if (rVariable == GREEN_LAGRANGE_STRAIN_TENSOR) { std::vector StrainVector; - CalculateOnIntegrationPoints( GREEN_LAGRANGE_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo ); + CalculateOnIntegrationPoints(GREEN_LAGRANGE_STRAIN_VECTOR, StrainVector, rCurrentProcessInfo); - //loop integration points for ( unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint ) { - if ( rOutput[GPoint].size2() != TDim ) - rOutput[GPoint].resize( TDim, TDim, false ); + if (rOutput[GPoint].size2() != TDim) + rOutput[GPoint].resize(TDim, TDim, false); rOutput[GPoint] = MathUtils::StrainVectorToTensor(StrainVector[GPoint]); } } else if (rVariable == PERMEABILITY_MATRIX) { //If the permeability of the element is a given property - BoundedMatrix PermeabilityMatrix; - GeoElementUtilities::FillPermeabilityMatrix(PermeabilityMatrix,this->GetProperties()); + BoundedMatrix PermeabilityMatrix; + GeoElementUtilities::FillPermeabilityMatrix(PermeabilityMatrix, this->GetProperties()); - //Loop over integration points - for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) - { + for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { rOutput[GPoint].resize(TDim,TDim,false); noalias(rOutput[GPoint]) = PermeabilityMatrix; } } else { - if ( rOutput.size() != mConstitutiveLawVector.size() ) + if (rOutput.size() != mConstitutiveLawVector.size()) rOutput.resize(mConstitutiveLawVector.size()); - for ( unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i ) { - rOutput[i].resize(TDim,TDim,false); - noalias(rOutput[i]) = ZeroMatrix(TDim,TDim); - rOutput[i] = mConstitutiveLawVector[i]->GetValue( rVariable, rOutput[i] ); + for (unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i) { + rOutput[i].resize(TDim, TDim, false); + noalias(rOutput[i]) = ZeroMatrix(TDim, TDim); + rOutput[i] = mConstitutiveLawVector[i]->GetValue(rVariable, rOutput[i]); } } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateMaterialStiffnessMatrix( MatrixType& rStiffnessMatrix, - const ProcessInfo& rCurrentProcessInfo ) +template +void UPwSmallStrainElement::CalculateMaterialStiffnessMatrix(MatrixType& rStiffnessMatrix, + const ProcessInfo& rCurrentProcessInfo) { KRATOS_TRY const IndexType N_DOF = TNumNodes * (TDim + 1); - //Resizing stiffness matrix - if ( rStiffnessMatrix.size1() != N_DOF ) - rStiffnessMatrix.resize( N_DOF, N_DOF, false ); - noalias( rStiffnessMatrix ) = ZeroMatrix( N_DOF, N_DOF ); + if (rStiffnessMatrix.size1() != N_DOF) + rStiffnessMatrix.resize(N_DOF, N_DOF, false); + noalias(rStiffnessMatrix) = ZeroMatrix(N_DOF, N_DOF); - //Previous definitions + // Previous definitions const PropertiesType& rProp = this->GetProperties(); const GeometryType& rGeom = this->GetGeometry(); - const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints( mThisIntegrationMethod ); + const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(mThisIntegrationMethod); const IndexType NumGPoints = IntegrationPoints.size(); - //Constitutive Law parameters + // Constitutive Law parameters ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom,rProp,rCurrentProcessInfo); ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR); ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); - //Element variables ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - //Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - //Compute Np, GradNpT, B and StrainVector + // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - //Compute infinitesimal strain + // Compute infinitesimal strain this->CalculateStrain(Variables, GPoint); - //set Gauss points variables to constitutive law parameters + // Set Gauss points variables to constitutive law parameters this->SetConstitutiveParameters(Variables, ConstitutiveParameters); - //Compute constitutive tensor + // Compute constitutive tensor noalias(Variables.StressVector) = mStressVector[GPoint]; ConstitutiveParameters.SetStressVector(Variables.StressVector); mConstitutiveLawVector[GPoint]->CalculateMaterialResponseCauchy(ConstitutiveParameters); - //Compute weighting coefficient for integration - Variables.IntegrationCoefficient = - this->CalculateIntegrationCoefficient(IntegrationPoints, - GPoint, - Variables.detJ); + // Compute weighting coefficient for integration + Variables.IntegrationCoefficient = this->CalculateIntegrationCoefficient(IntegrationPoints, + GPoint, + Variables.detJ); - //Compute stiffness matrix + // Compute stiffness matrix this->CalculateAndAddStiffnessMatrix(rStiffnessMatrix, Variables); } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddGeometricStiffnessMatrix( MatrixType& rLeftHandSideMatrix, - ElementVariables& rVariables, - unsigned int GPoint ) +template +void UPwSmallStrainElement::CalculateAndAddGeometricStiffnessMatrix(MatrixType& rLeftHandSideMatrix, + ElementVariables& rVariables, + unsigned int GPoint) { KRATOS_TRY - Matrix StressTensor = MathUtils::StressVectorToTensor( mStressVector[GPoint] ); + Matrix StressTensor = MathUtils::StressVectorToTensor(mStressVector[GPoint]); Matrix ReducedKgMatrix = prod( rVariables.GradNpT, rVariables.IntegrationCoefficient * Matrix( prod( StressTensor, trans(rVariables.GradNpT) ) ) ); //to be optimized Matrix UMatrix(TNumNodes*TDim, TNumNodes*TDim); noalias(UMatrix) = ZeroMatrix(TNumNodes*TDim, TNumNodes*TDim); - MathUtils::ExpandAndAddReducedMatrix( UMatrix, ReducedKgMatrix, TDim ); + MathUtils::ExpandAndAddReducedMatrix(UMatrix, ReducedKgMatrix, TDim); - //Distribute stiffness block matrix into the elemental matrix + // Distribute stiffness block matrix into the elemental matrix GeoElementUtilities::AssembleUBlockMatrix(rLeftHandSideMatrix, UMatrix, TNumNodes, TDim); - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateMassMatrix( MatrixType& rMassMatrix, - const ProcessInfo& rCurrentProcessInfo ) +template +void UPwSmallStrainElement::CalculateMassMatrix(MatrixType& rMassMatrix, + const ProcessInfo& rCurrentProcessInfo) { KRATOS_TRY const IndexType N_DOF = this->GetNumberOfDOF(); - //Resizing mass matrix - if ( rMassMatrix.size1() != N_DOF ) - rMassMatrix.resize( N_DOF, N_DOF, false ); - noalias( rMassMatrix ) = ZeroMatrix( N_DOF, N_DOF ); + if (rMassMatrix.size1() != N_DOF) + rMassMatrix.resize(N_DOF, N_DOF, false); + noalias(rMassMatrix) = ZeroMatrix(N_DOF, N_DOF); const GeometryType& rGeom = this->GetGeometry(); const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(this->mThisIntegrationMethod); const IndexType NumGPoints = IntegrationPoints.size(); - //Element variables ElementVariables Variables; - this->InitializeElementVariables( Variables, - rCurrentProcessInfo ); + this->InitializeElementVariables(Variables, + rCurrentProcessInfo); - // create general parametes of retention law + // Create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); - //Defining shape functions at all integration points - //Defining necessary variables + // Defining shape functions at all integration points + // Defining necessary variables BoundedMatrix Nut = ZeroMatrix(TDim+1, TNumNodes*(TDim+1)); BoundedMatrix AuxDensityMatrix = ZeroMatrix(TDim+1, TNumNodes*(TDim+1)); BoundedMatrix DensityMatrix = ZeroMatrix(TDim+1, TDim+1); - //Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint ) { GeoElementUtilities::CalculateNuElementMatrix(Nut, Variables.NContainer, GPoint); @@ -1065,7 +1003,7 @@ void UPwSmallStrainElement:: DNu_DX0, GPoint); - //calculating weighting coefficient for integration + // Calculating weighting coefficient for integration Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, @@ -1075,110 +1013,99 @@ void UPwSmallStrainElement:: this->CalculateSoilDensity(Variables); - GeoElementUtilities:: - AssembleDensityMatrix< TDim >(DensityMatrix, Variables.Density); + GeoElementUtilities::AssembleDensityMatrix< TDim >(DensityMatrix, Variables.Density); noalias(AuxDensityMatrix) = prod(DensityMatrix, Nut); - //Adding contribution to Mass matrix - noalias(rMassMatrix) += prod(trans(Nut), AuxDensityMatrix) - * Variables.IntegrationCoefficientInitialConfiguration; + // Adding contribution to Mass matrix + noalias(rMassMatrix) += prod(trans(Nut), AuxDensityMatrix) + * Variables.IntegrationCoefficientInitialConfiguration; } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } - -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAll( MatrixType& rLeftHandSideMatrix, - VectorType& rRightHandSideVector, - const ProcessInfo& rCurrentProcessInfo, - const bool CalculateStiffnessMatrixFlag, - const bool CalculateResidualVectorFlag) +template +void UPwSmallStrainElement::CalculateAll(MatrixType& rLeftHandSideMatrix, + VectorType& rRightHandSideVector, + const ProcessInfo& rCurrentProcessInfo, + bool CalculateStiffnessMatrixFlag, + bool CalculateResidualVectorFlag) { KRATOS_TRY - //Previous definitions + // Previous definitions const PropertiesType& rProp = this->GetProperties(); const GeometryType& rGeom = this->GetGeometry(); const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(this->mThisIntegrationMethod); const IndexType NumGPoints = IntegrationPoints.size(); - //Constitutive Law parameters + // Constitutive Law parameters ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom, rProp, rCurrentProcessInfo); - // stiffness matrix is needed to calculate Biot coefficient + // Stiffness matrix is needed to calculate Biot coefficient ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR); if (CalculateResidualVectorFlag) ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_STRESS); ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); - //Element variables ElementVariables Variables; - this->InitializeElementVariables( Variables, - rCurrentProcessInfo ); + this->InitializeElementVariables(Variables, + rCurrentProcessInfo); - // create general parameters of retention law + // Create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT); - //Loop over integration points for ( unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - //Compute GradNpT, B and StrainVector + // Compute GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - //Compute infinitesimal strain + // Compute infinitesimal strain this->CalculateStrain(Variables, GPoint); - //set Gauss points variables to constitutive law parameters + // Set Gauss points variables to constitutive law parameters this->SetConstitutiveParameters(Variables, ConstitutiveParameters); - //Compute Nu and BodyAcceleration + // Compute Nu and BodyAcceleration GeoElementUtilities::CalculateNuMatrix(Variables.Nu, Variables.NContainer, GPoint); - GeoElementUtilities:: - InterpolateVariableWithComponents( Variables.BodyAcceleration, - Variables.NContainer, - Variables.VolumeAcceleration, - GPoint ); + GeoElementUtilities::InterpolateVariableWithComponents(Variables.BodyAcceleration, + Variables.NContainer, + Variables.VolumeAcceleration, + GPoint); - - //Compute constitutive tensor and stresses + // Compute constitutive tensor and stresses ConstitutiveParameters.SetStressVector(mStressVector[GPoint]); mConstitutiveLawVector[GPoint]->CalculateMaterialResponseCauchy(ConstitutiveParameters); CalculateRetentionResponse(Variables, RetentionParameters, GPoint); this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); - this->CalculatePermeabilityUpdateFactor( Variables ); + this->CalculatePermeabilityUpdateFactor(Variables); - //Compute weighting coefficient for integration - Variables.IntegrationCoefficient = - this->CalculateIntegrationCoefficient(IntegrationPoints, - GPoint, - Variables.detJ); + // Compute weighting coefficient for integration + Variables.IntegrationCoefficient = this->CalculateIntegrationCoefficient(IntegrationPoints, + GPoint, + Variables.detJ); Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, Variables.detJInitialConfiguration); - //Contributions to the left hand side + // Contributions to the left hand side if (CalculateStiffnessMatrixFlag) this->CalculateAndAddLHS(rLeftHandSideMatrix, Variables); - //Contributions to the right hand side + // Contributions to the right hand side if (CalculateResidualVectorFlag) this->CalculateAndAddRHS(rRightHandSideVector, Variables, GPoint); } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -double UPwSmallStrainElement:: - CalculateBiotCoefficient( const ElementVariables& rVariables, - const bool &hasBiotCoefficient) const +template +double UPwSmallStrainElement::CalculateBiotCoefficient(const ElementVariables& rVariables, + bool hasBiotCoefficient) const { KRATOS_TRY @@ -1188,25 +1115,23 @@ double UPwSmallStrainElement:: if (hasBiotCoefficient) { return rProp[BIOT_COEFFICIENT]; } else { - // calculate Bulk modulus from stiffness matrix + // Calculate Bulk modulus from stiffness matrix const double BulkModulus = CalculateBulkModulus(rVariables.ConstitutiveMatrix); return 1.0 - BulkModulus / rProp[BULK_MODULUS_SOLID]; } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - InitializeBiotCoefficients( ElementVariables& rVariables, - const bool &hasBiotCoefficient) +template +void UPwSmallStrainElement::InitializeBiotCoefficients(ElementVariables& rVariables, + bool hasBiotCoefficient) { KRATOS_TRY const PropertiesType& rProp = this->GetProperties(); - //Properties variables + // Properties variables rVariables.BiotCoefficient = CalculateBiotCoefficient(rVariables, hasBiotCoefficient); if (!rProp[IGNORE_UNDRAINED]) { @@ -1220,13 +1145,11 @@ void UPwSmallStrainElement:: rVariables.BiotModulusInverse *= rVariables.DegreeOfSaturation; rVariables.BiotModulusInverse -= rVariables.DerivativeOfSaturation*rProp[POROSITY]; - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculatePermeabilityUpdateFactor( ElementVariables &rVariables) +template +void UPwSmallStrainElement::CalculatePermeabilityUpdateFactor(ElementVariables& rVariables) { KRATOS_TRY @@ -1244,85 +1167,71 @@ void UPwSmallStrainElement:: rVariables.PermeabilityUpdateFactor = 1.0; } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -double UPwSmallStrainElement:: - CalculateBulkModulus(const Matrix &ConstitutiveMatrix) const +template +double UPwSmallStrainElement::CalculateBulkModulus(const Matrix& ConstitutiveMatrix) const { KRATOS_TRY const int IndexG = ConstitutiveMatrix.size1() - 1; + return ConstitutiveMatrix(0, 0) - (4.0/3.0)*ConstitutiveMatrix(IndexG, IndexG); - const double M = ConstitutiveMatrix(0, 0); - const double G = ConstitutiveMatrix(IndexG, IndexG); - - return M - (4.0/3.0)*G; - - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - InitializeElementVariables( ElementVariables& rVariables, - const ProcessInfo& rCurrentProcessInfo ) +template +void UPwSmallStrainElement::InitializeElementVariables(ElementVariables& rVariables, + const ProcessInfo& rCurrentProcessInfo) { KRATOS_TRY - //Properties variables - this->InitializeProperties( rVariables ); + // Properties variables + this->InitializeProperties(rVariables); - //ProcessInfo variables + // ProcessInfo variables rVariables.VelocityCoefficient = rCurrentProcessInfo[VELOCITY_COEFFICIENT]; rVariables.DtPressureCoefficient = rCurrentProcessInfo[DT_PRESSURE_COEFFICIENT]; - //Nodal Variables - this->InitializeNodalDisplacementVariables( rVariables ); - this->InitializeNodalPorePressureVariables( rVariables ); - this->InitializeNodalVolumeAccelerationVariables( rVariables ); + // Nodal Variables + this->InitializeNodalDisplacementVariables(rVariables); + this->InitializeNodalPorePressureVariables(rVariables); + this->InitializeNodalVolumeAccelerationVariables(rVariables); - //Variables computed at each GP - noalias(rVariables.Nu) = ZeroMatrix(TDim, TNumNodes*TDim); + // Variables computed at each GP + rVariables.Nu = ZeroMatrix(TDim, TNumNodes*TDim); rVariables.Np.resize(TNumNodes,false); rVariables.GradNpT.resize(TNumNodes,TDim,false); - rVariables.F.resize(TDim,TDim,false); - - noalias(rVariables.F) = identity_matrix(TDim); + rVariables.F = identity_matrix(TDim); rVariables.detF = 1.0; - //General Variables - rVariables.VoigtVector.resize(VoigtSize); - noalias(rVariables.VoigtVector) = ZeroVector(VoigtSize); - for (unsigned int i=0; i < StressTensorSize; ++i) rVariables.VoigtVector[i] = 1.0; + // General Variables + rVariables.VoigtVector = ZeroVector(VoigtSize); + std::fill_n(rVariables.VoigtVector.begin(), StressTensorSize, 1.0); - rVariables.B.resize(VoigtSize,TNumNodes*TDim,false); - noalias(rVariables.B) = ZeroMatrix(VoigtSize,TNumNodes*TDim); + rVariables.B = ZeroMatrix(VoigtSize,TNumNodes*TDim); const GeometryType& rGeom = this->GetGeometry(); - const IndexType NumGPoints = rGeom.IntegrationPointsNumber( mThisIntegrationMethod ); + const IndexType NumGPoints = rGeom.IntegrationPointsNumber(mThisIntegrationMethod); - // shape functions - (rVariables.NContainer).resize(NumGPoints, TNumNodes, false); - rVariables.NContainer = rGeom.ShapeFunctionsValues( mThisIntegrationMethod ); + // Shape functions + rVariables.NContainer = rGeom.ShapeFunctionsValues(mThisIntegrationMethod); - // gradient of shape functions and determinant of Jacobian - (rVariables.detJContainer).resize(NumGPoints,false); + // Gradient of shape functions and determinant of Jacobian + rVariables.detJContainer.resize(NumGPoints,false); - rGeom.ShapeFunctionsIntegrationPointsGradients( rVariables.DN_DXContainer, + rGeom.ShapeFunctionsIntegrationPointsGradients(rVariables.DN_DXContainer, rVariables.detJContainer, - mThisIntegrationMethod ); + mThisIntegrationMethod); - //Constitutive Law parameters - + // Constitutive Law parameters rVariables.StressVector.resize(VoigtSize, false); rVariables.StrainVector.resize(VoigtSize, false); rVariables.ConstitutiveMatrix.resize(VoigtSize, VoigtSize, false); - //Auxiliary variables + // Auxiliary variables rVariables.UVoigtMatrix.resize(TNumNodes*TDim, VoigtSize, false); // Retention law @@ -1332,22 +1241,19 @@ void UPwSmallStrainElement:: rVariables.RelativePermeability = 1.0; rVariables.BishopCoefficient = 1.0; - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateBMatrix(Matrix& rB, - const Matrix& GradNpT, - const Vector &Np) +template +void UPwSmallStrainElement::CalculateBMatrix(Matrix& rB, + const Matrix& GradNpT, + const Vector& Np) { KRATOS_TRY unsigned int index; - if constexpr (TDim > 2) { - for ( unsigned int i = 0; i < TNumNodes; ++i ) { + for (unsigned int i = 0; i < TNumNodes; ++i) { index = TDim * i; rB( INDEX_3D_XX, index + INDEX_X ) = GradNpT( i, INDEX_X ); @@ -1362,7 +1268,7 @@ void UPwSmallStrainElement:: } } else { // 2D plane strain - for ( unsigned int i = 0; i < TNumNodes; ++i ) { + for (unsigned int i = 0; i < TNumNodes; ++i) { index = TDim * i; rB( INDEX_2D_PLANE_STRAIN_XX, index + INDEX_X ) = GradNpT( i, INDEX_X ); @@ -1372,15 +1278,14 @@ void UPwSmallStrainElement:: } } - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddLHS(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables) +template +void UPwSmallStrainElement::CalculateAndAddLHS(MatrixType& rLeftHandSideMatrix, + ElementVariables& rVariables) { - KRATOS_TRY; + KRATOS_TRY this->CalculateAndAddStiffnessMatrix(rLeftHandSideMatrix,rVariables); this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix,rVariables); @@ -1390,33 +1295,31 @@ void UPwSmallStrainElement:: this->CalculateAndAddPermeabilityMatrix(rLeftHandSideMatrix,rVariables); } - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddStiffnessMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables) +template +void UPwSmallStrainElement::CalculateAndAddStiffnessMatrix(MatrixType& rLeftHandSideMatrix, + ElementVariables& rVariables) { - KRATOS_TRY; + KRATOS_TRY noalias(rVariables.UVoigtMatrix) = prod(trans(rVariables.B), rVariables.ConstitutiveMatrix); - noalias(rVariables.UMatrix) = prod(rVariables.UVoigtMatrix, rVariables.B)*rVariables.IntegrationCoefficient; + noalias(rVariables.UMatrix) = prod(rVariables.UVoigtMatrix, rVariables.B) * rVariables.IntegrationCoefficient; - //Distribute stiffness block matrix into the elemental matrix + // Distribute stiffness block matrix into the elemental matrix GeoElementUtilities::AssembleUBlockMatrix(rLeftHandSideMatrix, rVariables.UMatrix); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddCouplingMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables) +template +void UPwSmallStrainElement::CalculateAndAddCouplingMatrix(MatrixType& rLeftHandSideMatrix, + ElementVariables& rVariables) { - KRATOS_TRY; + KRATOS_TRY - noalias(rVariables.UVector) = prod(trans(rVariables.B),rVariables.VoigtVector); + noalias(rVariables.UVector) = prod(trans(rVariables.B), rVariables.VoigtVector); noalias(rVariables.UPMatrix) = PORE_PRESSURE_SIGN_FACTOR * rVariables.BiotCoefficient @@ -1424,7 +1327,7 @@ void UPwSmallStrainElement:: * outer_prod(rVariables.UVector,rVariables.Np) * rVariables.IntegrationCoefficient; - //Distribute coupling block matrix into the elemental matrix + // Distribute coupling block matrix into the elemental matrix GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, rVariables.UPMatrix); if (!rVariables.IgnoreUndrained) { @@ -1434,96 +1337,85 @@ void UPwSmallStrainElement:: * rVariables.VelocityCoefficient * trans(rVariables.UPMatrix); - //Distribute transposed coupling block matrix into the elemental matrix + // Distribute transposed coupling block matrix into the elemental matrix GeoElementUtilities::AssemblePUBlockMatrix(rLeftHandSideMatrix,rVariables.PUMatrix); } - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateCompressibilityMatrix(BoundedMatrix &PMatrix, - const ElementVariables &rVariables) const +template +void UPwSmallStrainElement::CalculateCompressibilityMatrix(BoundedMatrix& rPMatrix, + const ElementVariables& rVariables) const { - KRATOS_TRY; + KRATOS_TRY - noalias(PMatrix) = - PORE_PRESSURE_SIGN_FACTOR + noalias(rPMatrix) = - PORE_PRESSURE_SIGN_FACTOR * rVariables.DtPressureCoefficient * rVariables.BiotModulusInverse * outer_prod(rVariables.Np, rVariables.Np) * rVariables.IntegrationCoefficient; - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, - ElementVariables& rVariables) +template +void UPwSmallStrainElement::CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, + ElementVariables& rVariables) { - KRATOS_TRY; + KRATOS_TRY this->CalculateCompressibilityMatrix(rVariables.PMatrix, rVariables); - //Distribute compressibility block matrix into the elemental matrix - GeoElementUtilities:: - AssemblePBlockMatrix< TDim, TNumNodes >(rLeftHandSideMatrix, - rVariables.PMatrix); + // Distribute compressibility block matrix into the elemental matrix + GeoElementUtilities::AssemblePBlockMatrix(rLeftHandSideMatrix, + rVariables.PMatrix); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculatePermeabilityMatrix(BoundedMatrix &PDimMatrix, - BoundedMatrix &PMatrix, - const ElementVariables &rVariables) const +template +void UPwSmallStrainElement::CalculatePermeabilityMatrix(BoundedMatrix& rPDimMatrix, + BoundedMatrix& rPMatrix, + const ElementVariables& rVariables) const { - KRATOS_TRY; + KRATOS_TRY - noalias(PDimMatrix) = - PORE_PRESSURE_SIGN_FACTOR + noalias(rPDimMatrix) = - PORE_PRESSURE_SIGN_FACTOR * prod(rVariables.GradNpT, rVariables.PermeabilityMatrix); - noalias(PMatrix) = rVariables.DynamicViscosityInverse + noalias(rPMatrix) = rVariables.DynamicViscosityInverse * rVariables.RelativePermeability * rVariables.PermeabilityUpdateFactor - * prod(PDimMatrix, trans(rVariables.GradNpT)) + * prod(rPDimMatrix, trans(rVariables.GradNpT)) * rVariables.IntegrationCoefficient; - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddPermeabilityMatrix(MatrixType &rLeftHandSideMatrix, - ElementVariables &rVariables) +template +void UPwSmallStrainElement::CalculateAndAddPermeabilityMatrix(MatrixType &rLeftHandSideMatrix, + ElementVariables &rVariables) { - KRATOS_TRY; + KRATOS_TRY this->CalculatePermeabilityMatrix(rVariables.PDimMatrix, rVariables.PMatrix, rVariables); - //Distribute permeability block matrix into the elemental matrix - GeoElementUtilities:: - AssemblePBlockMatrix< TDim, TNumNodes >(rLeftHandSideMatrix, rVariables.PMatrix); + // Distribute permeability block matrix into the elemental matrix + GeoElementUtilities::AssemblePBlockMatrix(rLeftHandSideMatrix, + rVariables.PMatrix); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddRHS(VectorType& rRightHandSideVector, - ElementVariables& rVariables, - unsigned int GPoint) +template +void UPwSmallStrainElement::CalculateAndAddRHS(VectorType& rRightHandSideVector, + ElementVariables& rVariables, + unsigned int GPoint) { - KRATOS_TRY; + KRATOS_TRY this->CalculateAndAddStiffnessForce(rRightHandSideVector, rVariables, GPoint); @@ -1539,83 +1431,72 @@ void UPwSmallStrainElement:: this->CalculateAndAddFluidBodyFlow(rRightHandSideVector, rVariables); } - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddStiffnessForce( VectorType& rRightHandSideVector, - ElementVariables& rVariables, - unsigned int GPoint ) +template +void UPwSmallStrainElement::CalculateAndAddStiffnessForce(VectorType& rRightHandSideVector, + ElementVariables& rVariables, + unsigned int GPoint) { - KRATOS_TRY; + KRATOS_TRY noalias(rVariables.UVector) = -1.0 * prod(trans(rVariables.B), mStressVector[GPoint]) * rVariables.IntegrationCoefficient; - //Distribute stiffness block vector into elemental vector + // Distribute stiffness block vector into elemental vector GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, rVariables.UVector); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddMixBodyForce(VectorType& rRightHandSideVector, ElementVariables& rVariables) +template +void UPwSmallStrainElement::CalculateAndAddMixBodyForce(VectorType& rRightHandSideVector, + ElementVariables& rVariables) { - KRATOS_TRY; + KRATOS_TRY this->CalculateSoilGamma(rVariables); noalias(rVariables.UVector) = prod(trans(rVariables.Nu), rVariables.SoilGamma) * rVariables.IntegrationCoefficientInitialConfiguration; - //Distribute body force block vector into elemental vector - GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector,rVariables.UVector); + // Distribute body force block vector into elemental vector + GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, rVariables.UVector); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateSoilDensity(ElementVariables &rVariables) +template +void UPwSmallStrainElement::CalculateSoilDensity(ElementVariables &rVariables) { - KRATOS_TRY; - - rVariables.Density = ( rVariables.DegreeOfSaturation - * rVariables.Porosity - * rVariables.FluidDensity ) - + (1.0 - rVariables.Porosity)*rVariables.SolidDensity; + KRATOS_TRY - KRATOS_CATCH(""); + rVariables.Density = rVariables.DegreeOfSaturation * rVariables.Porosity * rVariables.FluidDensity + + (1.0 - rVariables.Porosity) * rVariables.SolidDensity; + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateSoilGamma(ElementVariables &rVariables) +template +void UPwSmallStrainElement::CalculateSoilGamma(ElementVariables &rVariables) { - KRATOS_TRY; + KRATOS_TRY this->CalculateSoilDensity(rVariables); noalias(rVariables.SoilGamma) = rVariables.Density * rVariables.BodyAcceleration; - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddCouplingTerms(VectorType& rRightHandSideVector, ElementVariables& rVariables) +template +void UPwSmallStrainElement::CalculateAndAddCouplingTerms(VectorType& rRightHandSideVector, + ElementVariables& rVariables) { - KRATOS_TRY; + KRATOS_TRY - noalias(rVariables.UVector) = prod(trans(rVariables.B),rVariables.VoigtVector); + noalias(rVariables.UVector) = prod(trans(rVariables.B), rVariables.VoigtVector); noalias(rVariables.UPMatrix) = - PORE_PRESSURE_SIGN_FACTOR * rVariables.BiotCoefficient @@ -1625,7 +1506,7 @@ void UPwSmallStrainElement:: noalias(rVariables.UVector) = prod(rVariables.UPMatrix, rVariables.PressureVector); - //Distribute coupling block vector 1 into elemental vector + // Distribute coupling block vector 1 into elemental vector GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, rVariables.UVector); if (!rVariables.IgnoreUndrained) { @@ -1634,136 +1515,123 @@ void UPwSmallStrainElement:: * SaturationCoefficient * prod(trans(rVariables.UPMatrix), rVariables.VelocityVector); - //Distribute coupling block vector 2 into elemental vector - GeoElementUtilities::AssemblePBlockVector< TDim, TNumNodes >(rRightHandSideVector, rVariables.PVector); + // Distribute coupling block vector 2 into elemental vector + GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, rVariables.PVector); } - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateCompressibilityFlow(BoundedMatrix &PMatrix, - array_1d &PVector, - const ElementVariables &rVariables) const +template +void UPwSmallStrainElement::CalculateCompressibilityFlow(BoundedMatrix& rPMatrix, + array_1d& rPVector, + const ElementVariables& rVariables) const { - KRATOS_TRY; + KRATOS_TRY - noalias(PMatrix) = - PORE_PRESSURE_SIGN_FACTOR + noalias(rPMatrix) = - PORE_PRESSURE_SIGN_FACTOR * rVariables.BiotModulusInverse * outer_prod(rVariables.Np,rVariables.Np) * rVariables.IntegrationCoefficient; - noalias(PVector) = - prod(PMatrix, rVariables.DtPressureVector); + noalias(rPVector) = - prod(rPMatrix, rVariables.DtPressureVector); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddCompressibilityFlow( VectorType& rRightHandSideVector, - ElementVariables& rVariables ) +template +void UPwSmallStrainElement::CalculateAndAddCompressibilityFlow(VectorType& rRightHandSideVector, + ElementVariables& rVariables) { - KRATOS_TRY; + KRATOS_TRY this->CalculateCompressibilityFlow(rVariables.PMatrix, rVariables.PVector, rVariables); - //Distribute compressibility block vector into elemental vector + // Distribute compressibility block vector into elemental vector GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, rVariables.PVector); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculatePermeabilityFlow(BoundedMatrix &PDimMatrix, - BoundedMatrix &PMatrix, - array_1d &PVector, - const ElementVariables &rVariables) const +template +void UPwSmallStrainElement::CalculatePermeabilityFlow(BoundedMatrix& rPDimMatrix, + BoundedMatrix& rPMatrix, + array_1d& rPVector, + const ElementVariables &rVariables) const { - KRATOS_TRY; + KRATOS_TRY - noalias(PDimMatrix) = prod(rVariables.GradNpT, rVariables.PermeabilityMatrix); + noalias(rPDimMatrix) = prod(rVariables.GradNpT, rVariables.PermeabilityMatrix); - noalias(PMatrix) = - PORE_PRESSURE_SIGN_FACTOR + noalias(rPMatrix) = - PORE_PRESSURE_SIGN_FACTOR * rVariables.DynamicViscosityInverse * rVariables.RelativePermeability * rVariables.PermeabilityUpdateFactor - * prod(PDimMatrix,trans(rVariables.GradNpT)) + * prod(rPDimMatrix,trans(rVariables.GradNpT)) * rVariables.IntegrationCoefficient; - noalias(PVector) = - prod(PMatrix, rVariables.PressureVector); + noalias(rPVector) = - prod(rPMatrix, rVariables.PressureVector); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddPermeabilityFlow( VectorType &rRightHandSideVector, - ElementVariables &rVariables ) +template +void UPwSmallStrainElement::CalculateAndAddPermeabilityFlow(VectorType &rRightHandSideVector, + ElementVariables &rVariables) { - KRATOS_TRY; + KRATOS_TRY this->CalculatePermeabilityFlow(rVariables.PDimMatrix, rVariables.PMatrix, rVariables.PVector, rVariables); - //Distribute permeability block vector into elemental vector + // Distribute permeability block vector into elemental vector GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, rVariables.PVector); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateFluidBodyFlow(BoundedMatrix &PDimMatrix, - array_1d &PVector, - const ElementVariables &rVariables) const +template +void UPwSmallStrainElement::CalculateFluidBodyFlow(BoundedMatrix& rPDimMatrix, + array_1d& rPVector, + const ElementVariables& rVariables) const { - KRATOS_TRY; + KRATOS_TRY - noalias(PDimMatrix) = prod(rVariables.GradNpT,rVariables.PermeabilityMatrix) + noalias(rPDimMatrix) = prod(rVariables.GradNpT, rVariables.PermeabilityMatrix) * rVariables.IntegrationCoefficient; - noalias(PVector) = rVariables.DynamicViscosityInverse + noalias(rPVector) = rVariables.DynamicViscosityInverse * rVariables.FluidDensity * rVariables.RelativePermeability * rVariables.PermeabilityUpdateFactor - * prod(PDimMatrix,rVariables.BodyAcceleration); + * prod(rPDimMatrix, rVariables.BodyAcceleration); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateAndAddFluidBodyFlow(VectorType& rRightHandSideVector, - ElementVariables& rVariables) +template +void UPwSmallStrainElement::CalculateAndAddFluidBodyFlow(VectorType& rRightHandSideVector, + ElementVariables& rVariables) { - KRATOS_TRY; + KRATOS_TRY this->CalculateFluidBodyFlow(rVariables.PDimMatrix, rVariables.PVector, rVariables); - //Distribute fluid body flow block vector into elemental vector + // Distribute fluid body flow block vector into elemental vector GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector,rVariables.PVector); - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateStrain( ElementVariables& rVariables, unsigned int GPoint ) +template +void UPwSmallStrainElement::CalculateStrain(ElementVariables& rVariables, + unsigned int GPoint) { if (rVariables.UseHenckyStrain) { this->CalculateDeformationGradient(rVariables, GPoint); @@ -1773,22 +1641,18 @@ void UPwSmallStrainElement:: } } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateCauchyStrain( ElementVariables& rVariables ) +template +void UPwSmallStrainElement::CalculateCauchyStrain(ElementVariables& rVariables) { noalias(rVariables.StrainVector) = prod(rVariables.B, rVariables.DisplacementVector); } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateCauchyGreenStrain( ElementVariables& rVariables ) +template +void UPwSmallStrainElement::CalculateCauchyGreenStrain(ElementVariables& rVariables) { - KRATOS_TRY; + KRATOS_TRY - //-Compute total deformation gradient + // Compute total deformation gradient const Matrix& F = rVariables.F; Matrix ETensor; @@ -1809,16 +1673,13 @@ void UPwSmallStrainElement:: noalias(rVariables.StrainVector) = MathUtils::StrainTensorToVector(ETensor); } - KRATOS_CATCH(""); + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateCauchyAlmansiStrain(ElementVariables& rVariables) +template +void UPwSmallStrainElement::CalculateCauchyAlmansiStrain(ElementVariables& rVariables) { - - //-Compute total deformation gradient + // Compute total deformation gradient const Matrix& F = rVariables.F; Matrix LeftCauchyGreen; @@ -1826,7 +1687,7 @@ void UPwSmallStrainElement:: Matrix ETensor; double det; - MathUtils::InvertMatrix(LeftCauchyGreen, ETensor, det ); + MathUtils::InvertMatrix(LeftCauchyGreen, ETensor, det); for (unsigned int i=0; i:: } -//---------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateDeformationGradient( ElementVariables& rVariables, - unsigned int GPoint) +template +void UPwSmallStrainElement::CalculateDeformationGradient(ElementVariables& rVariables, + unsigned int GPoint) { KRATOS_TRY - // calculation of derivative of shape function with respect to reference configuration + // Calculation of derivative of shape function with respect to reference configuration // derivative of shape function (displacement) Matrix J0, InvJ0, DNu_DX0; double detJ0; @@ -1864,7 +1723,7 @@ void UPwSmallStrainElement:: DNu_DX0, GPoint); - //Calculating current Jacobian in order to find deformation gradient + // Calculating current Jacobian in order to find deformation gradient Matrix J, InvJ, DNu_DX; double detJ; this->CalculateJacobianOnCurrentConfiguration(detJ, @@ -1884,62 +1743,56 @@ void UPwSmallStrainElement:: noalias(rVariables.F) = prod( J, InvJ0 ); rVariables.detF = MathUtils::Det(rVariables.F); - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - InitializeNodalPorePressureVariables( ElementVariables& rVariables ) +template +void UPwSmallStrainElement::InitializeNodalPorePressureVariables(ElementVariables& rVariables) { KRATOS_TRY const GeometryType& rGeom = this->GetGeometry(); - //Nodal Variables + // Nodal variables for (unsigned int i=0; i -void UPwSmallStrainElement:: - InitializeNodalDisplacementVariables( ElementVariables& rVariables ) +template +void UPwSmallStrainElement::InitializeNodalDisplacementVariables(ElementVariables& rVariables) { KRATOS_TRY const GeometryType& rGeom = this->GetGeometry(); - //Nodal Variables + // Nodal variables GeoElementUtilities::GetNodalVariableVector(rVariables.DisplacementVector, rGeom, DISPLACEMENT); GeoElementUtilities::GetNodalVariableVector(rVariables.VelocityVector, rGeom, VELOCITY); - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - InitializeNodalVolumeAccelerationVariables( ElementVariables& rVariables ) +template +void UPwSmallStrainElement::InitializeNodalVolumeAccelerationVariables(ElementVariables& rVariables) { KRATOS_TRY const GeometryType& rGeom = this->GetGeometry(); - //Nodal Variables - GeoElementUtilities::GetNodalVariableVector(rVariables.VolumeAcceleration, rGeom, VOLUME_ACCELERATION); + // Nodal variables + GeoElementUtilities::GetNodalVariableVector(rVariables.VolumeAcceleration, + rGeom, + VOLUME_ACCELERATION); - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - InitializeProperties( ElementVariables& rVariables ) +template +void UPwSmallStrainElement::InitializeProperties(ElementVariables& rVariables) { KRATOS_TRY @@ -1954,7 +1807,7 @@ void UPwSmallStrainElement:: if (rProp.Has(CONSIDER_GEOMETRIC_STIFFNESS)) rVariables.ConsiderGeometricStiffness = rProp[CONSIDER_GEOMETRIC_STIFFNESS]; - rVariables.DynamicViscosityInverse = 1.0/rProp[DYNAMIC_VISCOSITY]; + rVariables.DynamicViscosityInverse = 1.0 / rProp[DYNAMIC_VISCOSITY]; rVariables.FluidDensity = rProp[DENSITY_WATER]; rVariables.SolidDensity = rProp[DENSITY_SOLID]; rVariables.Porosity = rProp[POROSITY]; @@ -1962,24 +1815,22 @@ void UPwSmallStrainElement:: rVariables.PermeabilityUpdateFactor = 1.0; - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateKinematics( ElementVariables& rVariables, - unsigned int GPoint ) +template +void UPwSmallStrainElement::CalculateKinematics(ElementVariables& rVariables, + unsigned int GPoint) { KRATOS_TRY - //Setting the vector of shape functions and the matrix of the shape functions global gradients + // Setting the vector of shape functions and the matrix of the shape functions global gradients noalias(rVariables.Np) = row(rVariables.NContainer, GPoint); noalias(rVariables.GradNpT) = rVariables.DN_DXContainer[GPoint]; rVariables.detJ = rVariables.detJContainer[GPoint]; - //Compute the deformation matrix B + // Compute the deformation matrix B this->CalculateBMatrix(rVariables.B, rVariables.GradNpT, rVariables.Np); Matrix J0,InvJ0; @@ -1989,14 +1840,12 @@ void UPwSmallStrainElement:: rVariables.GradNpTInitialConfiguration, GPoint); - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - SetConstitutiveParameters(ElementVariables& rVariables, - ConstitutiveLaw::Parameters& rConstitutiveParameters) +template +void UPwSmallStrainElement::SetConstitutiveParameters(ElementVariables& rVariables, + ConstitutiveLaw::Parameters& rConstitutiveParameters) { KRATOS_TRY @@ -2007,40 +1856,34 @@ void UPwSmallStrainElement:: rConstitutiveParameters.SetDeformationGradientF(rVariables.F); rConstitutiveParameters.SetDeterminantF(rVariables.detF); - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - SetRetentionParameters(const ElementVariables& rVariables, - RetentionLaw::Parameters& rRetentionParameters) +template +void UPwSmallStrainElement::SetRetentionParameters(const ElementVariables& rVariables, + RetentionLaw::Parameters& rRetentionParameters) { KRATOS_TRY rRetentionParameters.SetFluidPressure(rVariables.FluidPressure); - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -double UPwSmallStrainElement:: - CalculateFluidPressure(const ElementVariables &rVariables) +template +double UPwSmallStrainElement::CalculateFluidPressure(const ElementVariables& rVariables) { KRATOS_TRY return inner_prod(rVariables.Np, rVariables.PressureVector); - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateRetentionResponse( ElementVariables& rVariables, - RetentionLaw::Parameters& rRetentionParameters, - unsigned int GPoint ) +template +void UPwSmallStrainElement::CalculateRetentionResponse(ElementVariables& rVariables, + RetentionLaw::Parameters& rRetentionParameters, + unsigned int GPoint) { KRATOS_TRY @@ -2053,12 +1896,11 @@ void UPwSmallStrainElement:: rVariables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(rRetentionParameters); rVariables.EffectiveSaturation = mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(rRetentionParameters); - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< unsigned int TDim, unsigned int TNumNodes > -void UPwSmallStrainElement:: - CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix) + +template +void UPwSmallStrainElement::CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix) { KRATOS_TRY @@ -2068,62 +1910,50 @@ void UPwSmallStrainElement:: << this->Id() << std::endl; - KRATOS_CATCH( "" ) + KRATOS_CATCH("") } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< > -void UPwSmallStrainElement< 2, 3 >:: - CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix) +template < > +void UPwSmallStrainElement<2, 3>::CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix) { - /// The matrix contains the shape functions at each GP evaluated at each node. - /// Rows: nodes - /// Columns: GP - - //Triangle_2d_3 - //GI_GAUSS_2 + // The matrix contains the shape functions at each GP evaluated at each node. + // Rows: nodes + // Columns: GP + // Triangle_2d_3 + // GI_GAUSS_2 rExtrapolationMatrix(0,0) = 1.6666666666666666666; rExtrapolationMatrix(0,1) = -0.33333333333333333333; rExtrapolationMatrix(0,2) = -0.33333333333333333333; rExtrapolationMatrix(1,0) = -0.33333333333333333333; rExtrapolationMatrix(1,1) = 1.6666666666666666666; rExtrapolationMatrix(1,2) = -0.33333333333333333333; rExtrapolationMatrix(2,0) = -0.33333333333333333333; rExtrapolationMatrix(2,1) = -0.33333333333333333333; rExtrapolationMatrix(2,2) = 1.6666666666666666666; } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< > -void UPwSmallStrainElement< 2, 4 >:: - CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix) +template < > +void UPwSmallStrainElement<2, 4>::CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix) { - //Quadrilateral_2d_4 - //GI_GAUSS_2 - + // Quadrilateral_2d_4 + // GI_GAUSS_2 rExtrapolationMatrix(0,0) = 1.8660254037844386; rExtrapolationMatrix(0,1) = -0.5; rExtrapolationMatrix(0,2) = 0.13397459621556132; rExtrapolationMatrix(0,3) = -0.5; rExtrapolationMatrix(1,0) = -0.5; rExtrapolationMatrix(1,1) = 1.8660254037844386; rExtrapolationMatrix(1,2) = -0.5; rExtrapolationMatrix(1,3) = 0.13397459621556132; rExtrapolationMatrix(2,0) = 0.13397459621556132; rExtrapolationMatrix(2,1) = -0.5; rExtrapolationMatrix(2,2) = 1.8660254037844386; rExtrapolationMatrix(2,3) = -0.5; rExtrapolationMatrix(3,0) = -0.5; rExtrapolationMatrix(3,1) = 0.13397459621556132; rExtrapolationMatrix(3,2) = -0.5; rExtrapolationMatrix(3,3) = 1.8660254037844386; } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< > -void UPwSmallStrainElement< 3, 4 >:: - CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix) +template < > +void UPwSmallStrainElement<3, 4>::CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix) { - //Tetrahedra_3d_4 - //GI_GAUSS_2 - + // Tetrahedra_3d_4 + // GI_GAUSS_2 rExtrapolationMatrix(0,0) = -0.309016988749894905; rExtrapolationMatrix(0,1) = -0.3090169887498949046; rExtrapolationMatrix(0,2) = -0.309016988749894905; rExtrapolationMatrix(0,3) = 1.9270509662496847144; rExtrapolationMatrix(1,0) = 1.9270509662496847144; rExtrapolationMatrix(1,1) = -0.30901698874989490481; rExtrapolationMatrix(1,2) = -0.3090169887498949049; rExtrapolationMatrix(1,3) = -0.30901698874989490481; rExtrapolationMatrix(2,0) = -0.30901698874989490473; rExtrapolationMatrix(2,1) = 1.9270509662496847143; rExtrapolationMatrix(2,2) = -0.3090169887498949049; rExtrapolationMatrix(2,3) = -0.30901698874989490481; rExtrapolationMatrix(3,0) = -0.3090169887498949048; rExtrapolationMatrix(3,1) = -0.30901698874989490471; rExtrapolationMatrix(3,2) = 1.9270509662496847143; rExtrapolationMatrix(3,3) = -0.30901698874989490481; } -//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -template< > -void UPwSmallStrainElement< 3, 8 >:: - CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix) +template < > +void UPwSmallStrainElement<3, 8>::CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix) { - //Hexahedra_3d_8 - //GI_GAUSS_2 - + // Hexahedra_3d_8 + // GI_GAUSS_2 rExtrapolationMatrix(0,0) = 2.549038105676658; rExtrapolationMatrix(0,1) = -0.6830127018922192; rExtrapolationMatrix(0,2) = 0.18301270189221927; rExtrapolationMatrix(0,3) = -0.6830127018922192; rExtrapolationMatrix(0,4) = -0.6830127018922192; rExtrapolationMatrix(0,5) = 0.18301270189221927; rExtrapolationMatrix(0,6) = -0.04903810567665795; rExtrapolationMatrix(0,7) = 0.18301270189221927; @@ -2149,8 +1979,6 @@ void UPwSmallStrainElement< 3, 8 >:: rExtrapolationMatrix(7,4) = -0.6830127018922192; rExtrapolationMatrix(7,5) = 0.18301270189221927; rExtrapolationMatrix(7,6) = -0.6830127018922192; rExtrapolationMatrix(7,7) = 2.549038105676658; } -//---------------------------------------------------------------------------------------------------- - template class UPwSmallStrainElement<2,3>; template class UPwSmallStrainElement<2,4>; template class UPwSmallStrainElement<3,4>; diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp index 7626c10366d9..2bcba0a36b4f 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -25,9 +25,9 @@ namespace Kratos { -template< unsigned int TDim, unsigned int TNumNodes > +template class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : - public UPwBaseElement + public UPwBaseElement { public: @@ -43,31 +43,28 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : using MatrixType = Matrix; /// The definition of the sizetype using SizeType = std::size_t; - using UPwBaseElement::mConstitutiveLawVector; - using UPwBaseElement::mRetentionLawVector; - using UPwBaseElement::mStressVector; - using UPwBaseElement::mStateVariablesFinalized; - using UPwBaseElement::mIsInitialised; - using UPwBaseElement::CalculateDerivativesOnInitialConfiguration; - using UPwBaseElement::mThisIntegrationMethod; + using UPwBaseElement::mConstitutiveLawVector; + using UPwBaseElement::mRetentionLawVector; + using UPwBaseElement::mStressVector; + using UPwBaseElement::mStateVariablesFinalized; + using UPwBaseElement::mIsInitialised; + using UPwBaseElement::CalculateDerivativesOnInitialConfiguration; + using UPwBaseElement::mThisIntegrationMethod; -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - - /// Default Constructor - explicit UPwSmallStrainElement(IndexType NewId = 0) : UPwBaseElement( NewId ) {} + explicit UPwSmallStrainElement(IndexType NewId = 0) : UPwBaseElement(NewId) {} /// Constructor using an array of nodes UPwSmallStrainElement(IndexType NewId, - const NodesArrayType& ThisNodes) : UPwBaseElement(NewId, ThisNodes) {} + const NodesArrayType& ThisNodes) : UPwBaseElement(NewId, ThisNodes) {} /// Constructor using Geometry UPwSmallStrainElement(IndexType NewId, - GeometryType::Pointer pGeometry) : UPwBaseElement(NewId, pGeometry) {} + GeometryType::Pointer pGeometry) : UPwBaseElement(NewId, pGeometry) {} /// Constructor using Properties UPwSmallStrainElement(IndexType NewId, GeometryType::Pointer pGeometry, - PropertiesType::Pointer pProperties) : UPwBaseElement( NewId, pGeometry, pProperties ) {} + PropertiesType::Pointer pProperties) : UPwBaseElement(NewId, pGeometry, pProperties) {} ~UPwSmallStrainElement() override = default; UPwSmallStrainElement(const UPwSmallStrainElement&) = delete; @@ -75,20 +72,17 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : UPwSmallStrainElement(UPwSmallStrainElement&&) = delete; UPwSmallStrainElement& operator=(UPwSmallStrainElement&&) = delete; -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - Element::Pointer Create( IndexType NewId, - NodesArrayType const& ThisNodes, - PropertiesType::Pointer pProperties ) const override; + Element::Pointer Create(IndexType NewId, + NodesArrayType const& ThisNodes, + PropertiesType::Pointer pProperties) const override; - Element::Pointer Create( IndexType NewId, - GeometryType::Pointer pGeom, - PropertiesType::Pointer pProperties ) const override; + Element::Pointer Create(IndexType NewId, + GeometryType::Pointer pGeom, + PropertiesType::Pointer pProperties) const override; int Check(const ProcessInfo& rCurrentProcessInfo) const override; -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - void InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override; void InitializeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) override; @@ -97,7 +91,11 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : void FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override; -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + void SetValuesOnIntegrationPoints(const Variable& rVariable, + const std::vector& rValues, + const ProcessInfo& rCurrentProcessInfo) override; + + using UPwBaseElement::SetValuesOnIntegrationPoints; void CalculateOnIntegrationPoints(const Variable& rVariable, std::vector& rOutput, @@ -115,22 +113,16 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : std::vector& rOutput, const ProcessInfo& rCurrentProcessInfo) override; - // Turn back information as a string. std::string Info() const override { - std::stringstream buffer; - buffer << "U-Pw small strain Element #" << this->Id() << "\nConstitutive law: " << mConstitutiveLawVector[0]->Info(); - return buffer.str(); + return "U-Pw small strain Element #" + std::to_string(this->Id()) + "\nConstitutive law: " + mConstitutiveLawVector[0]->Info(); } - // Print information about this object. void PrintInfo(std::ostream& rOStream) const override { - rOStream << "U-Pw small strain Element #" << this->Id() << "\nConstitutive law: " << mConstitutiveLawVector[0]->Info(); + rOStream << Info(); } -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - protected: static constexpr SizeType VoigtSize = ( TDim == N_DIM_3D ? VOIGT_SIZE_3D : VOIGT_SIZE_2D_PLANE_STRAIN); @@ -151,18 +143,18 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : double BiotCoefficient; double BiotModulusInverse; - BoundedMatrix PermeabilityMatrix; + BoundedMatrix PermeabilityMatrix; ///ProcessInfo variables double VelocityCoefficient; double DtPressureCoefficient; ///Nodal variables - array_1d PressureVector; - array_1d DtPressureVector; - array_1d DisplacementVector; - array_1d VelocityVector; - array_1d VolumeAcceleration; + array_1d PressureVector; + array_1d DtPressureVector; + array_1d DisplacementVector; + array_1d VelocityVector; + array_1d VolumeAcceleration; ///General elemental variables Vector VoigtVector; @@ -202,146 +194,141 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : double IntegrationCoefficientInitialConfiguration; //Auxiliary Variables - BoundedMatrix UMatrix; - BoundedMatrix UPMatrix; - BoundedMatrix PUMatrix; - BoundedMatrix PMatrix; + BoundedMatrix UMatrix; + BoundedMatrix UPMatrix; + BoundedMatrix PUMatrix; + BoundedMatrix PMatrix; Matrix UVoigtMatrix; - BoundedMatrix PDimMatrix; - array_1d UVector; - array_1d PVector; + BoundedMatrix PDimMatrix; + array_1d UVector; + array_1d PVector; }; - void SaveGPStress( Matrix &rStressContainer, - const Vector &StressVector, - unsigned int GPoint ); + void SaveGPStress(Matrix &rStressContainer, + const Vector& rStressVector, + unsigned int GPoint); - void ExtrapolateGPValues( const Matrix &StressContainer); + void ExtrapolateGPValues(const Matrix& rStressContainer); - void CalculateMaterialStiffnessMatrix( MatrixType& rStiffnessMatrix, - const ProcessInfo& CurrentProcessInfo ) override; + void CalculateMaterialStiffnessMatrix(MatrixType& rStiffnessMatrix, + const ProcessInfo& CurrentProcessInfo) override; - void CalculateMassMatrix( MatrixType& rMassMatrix, - const ProcessInfo& rCurrentProcessInfo ) override; + void CalculateMassMatrix(MatrixType& rMassMatrix, + const ProcessInfo& rCurrentProcessInfo) override; void CalculateAll( MatrixType &rLeftHandSideMatrix, VectorType &rRightHandSideVector, const ProcessInfo& CurrentProcessInfo, - const bool CalculateStiffnessMatrixFlag, - const bool CalculateResidualVectorFlag ) override; + bool CalculateStiffnessMatrixFlag, + bool CalculateResidualVectorFlag) override; - virtual void InitializeElementVariables( ElementVariables &rVariables, - const ProcessInfo &CurrentProcessInfo ); + virtual void InitializeElementVariables(ElementVariables& rVariables, + const ProcessInfo& CurrentProcessInfo); - void SetConstitutiveParameters(ElementVariables &rVariables, - ConstitutiveLaw::Parameters &rConstitutiveParameters); + void SetConstitutiveParameters(ElementVariables& rVariables, + ConstitutiveLaw::Parameters& rConstitutiveParameters); - void SetRetentionParameters(const ElementVariables &rVariables, - RetentionLaw::Parameters &rRetentionParameters); + void SetRetentionParameters(const ElementVariables& rVariables, + RetentionLaw::Parameters& rRetentionParameters); virtual void CalculateKinematics( ElementVariables &rVariables, unsigned int PointNumber ); - void InitializeBiotCoefficients( ElementVariables &rVariables, - const bool &hasBiotCoefficient=false ); + void InitializeBiotCoefficients(ElementVariables& rVariables, + bool hasBiotCoefficient=false ); - void CalculatePermeabilityUpdateFactor( ElementVariables &rVariables); + void CalculatePermeabilityUpdateFactor(ElementVariables &rVariables); virtual void CalculateBMatrix( Matrix &rB, const Matrix &GradNpT, const Vector &Np ); - virtual void CalculateAndAddLHS(MatrixType &rLeftHandSideMatrix, ElementVariables &rVariables); + virtual void CalculateAndAddLHS(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables); - void CalculateAndAddStiffnessMatrix(MatrixType &rLeftHandSideMatrix, ElementVariables &rVariables); + void CalculateAndAddStiffnessMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables); - void CalculateAndAddCouplingMatrix(MatrixType &rLeftHandSideMatrix, ElementVariables &rVariables); + void CalculateAndAddCouplingMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables); - virtual void CalculateAndAddCompressibilityMatrix(MatrixType &rLeftHandSideMatrix, + virtual void CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables); - virtual void CalculateCompressibilityMatrix(BoundedMatrix &PMatrix, - const ElementVariables &rVariables) const; + virtual void CalculateCompressibilityMatrix(BoundedMatrix& rPMatrix, + const ElementVariables& rVariables) const; - virtual void CalculateAndAddPermeabilityMatrix(MatrixType &rLeftHandSideMatrix, - ElementVariables &rVariables); + virtual void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, + ElementVariables& rVariables); - virtual void CalculatePermeabilityMatrix(BoundedMatrix &PDimMatrix, - BoundedMatrix &PMatrix, - const ElementVariables &rVariables) const; + virtual void CalculatePermeabilityMatrix(BoundedMatrix& rPDimMatrix, + BoundedMatrix& rPMatrix, + const ElementVariables& rVariables) const; - virtual void CalculateAndAddRHS(VectorType &rRightHandSideVector, - ElementVariables &rVariables, + virtual void CalculateAndAddRHS(VectorType& rRightHandSideVector, + ElementVariables& rVariables, unsigned int GPoint); - void CalculateAndAddStiffnessForce(VectorType &rRightHandSideVector, + void CalculateAndAddStiffnessForce(VectorType& rRightHandSideVector, ElementVariables& rVariables, unsigned int GPoint); - void CalculateAndAddMixBodyForce(VectorType &rRightHandSideVector, ElementVariables &rVariables); + void CalculateAndAddMixBodyForce(VectorType& rRightHandSideVector, ElementVariables& rVariables); - void CalculateAndAddCouplingTerms(VectorType& rRightHandSideVector, ElementVariables &rVariables); + void CalculateAndAddCouplingTerms(VectorType& rRightHandSideVector, ElementVariables& rVariables); - virtual void CalculateAndAddCompressibilityFlow(VectorType &rRightHandSideVector, - ElementVariables &rVariables); + virtual void CalculateAndAddCompressibilityFlow(VectorType& rRightHandSideVector, + ElementVariables& rVariables); - virtual void CalculateCompressibilityFlow(BoundedMatrix &PMatrix, - array_1d &PVector, - const ElementVariables &rVariables) const; + virtual void CalculateCompressibilityFlow(BoundedMatrix& rPMatrix, + array_1d& rPVector, + const ElementVariables& rVariables) const; - virtual void CalculateAndAddPermeabilityFlow(VectorType &rRightHandSideVector, - ElementVariables &rVariables); + virtual void CalculateAndAddPermeabilityFlow(VectorType& rRightHandSideVector, + ElementVariables& rVariables); - virtual void CalculatePermeabilityFlow(BoundedMatrix &PDimMatrix, - BoundedMatrix &PMatrix, - array_1d &PVector, - const ElementVariables &rVariables) const; + virtual void CalculatePermeabilityFlow(BoundedMatrix& rPDimMatrix, + BoundedMatrix& rPMatrix, + array_1d& rPVector, + const ElementVariables& rVariables) const; - virtual void CalculateAndAddFluidBodyFlow(VectorType &rRightHandSideVector, - ElementVariables &rVariables); - virtual void CalculateFluidBodyFlow(BoundedMatrix &PDimMatrix, - array_1d &PVector, - const ElementVariables &rVariables) const; + virtual void CalculateAndAddFluidBodyFlow(VectorType& rRightHandSideVector, + ElementVariables& rVariables); + virtual void CalculateFluidBodyFlow(BoundedMatrix& rPDimMatrix, + array_1d& rPVector, + const ElementVariables& rVariables) const; double CalculateBulkModulus(const Matrix &ConstitutiveMatrix) const; - double CalculateBiotCoefficient( const ElementVariables &rVariables, - const bool &hasBiotCoefficient) const; + double CalculateBiotCoefficient(const ElementVariables& rVariables, + bool hasBiotCoefficient) const; - virtual void CalculateCauchyAlmansiStrain( ElementVariables &rVariables ); - virtual void CalculateCauchyGreenStrain( ElementVariables &rVariables ); - virtual void CalculateCauchyStrain( ElementVariables &rVariables ); - virtual void CalculateStrain( ElementVariables &rVariables, unsigned int GPoint ); - virtual void CalculateDeformationGradient( ElementVariables& rVariables, - unsigned int GPoint ); + virtual void CalculateCauchyAlmansiStrain(ElementVariables& rVariables); + virtual void CalculateCauchyGreenStrain(ElementVariables& rVariables); + virtual void CalculateCauchyStrain(ElementVariables& rVariables); + virtual void CalculateStrain(ElementVariables &rVariables, unsigned int GPoint); + virtual void CalculateDeformationGradient(ElementVariables& rVariables, unsigned int GPoint); - void InitializeNodalDisplacementVariables( ElementVariables &rVariables ); - void InitializeNodalPorePressureVariables( ElementVariables &rVariables ); - void InitializeNodalVolumeAccelerationVariables( ElementVariables &rVariables ); + void InitializeNodalDisplacementVariables(ElementVariables& rVariables); + void InitializeNodalPorePressureVariables(ElementVariables& rVariables); + void InitializeNodalVolumeAccelerationVariables(ElementVariables& rVariables); - void InitializeProperties( ElementVariables &rVariables ); - double CalculateFluidPressure( const ElementVariables &rVariables); + void InitializeProperties(ElementVariables& rVariables); + double CalculateFluidPressure(const ElementVariables& rVariables); - void CalculateRetentionResponse( ElementVariables &rVariables, - RetentionLaw::Parameters &rRetentionParameters, - unsigned int GPoint ); + void CalculateRetentionResponse(ElementVariables& rVariables, + RetentionLaw::Parameters& rRetentionParameters, + unsigned int GPoint); - void CalculateExtrapolationMatrix(BoundedMatrix &rExtrapolationMatrix); + void CalculateExtrapolationMatrix(BoundedMatrix& rExtrapolationMatrix); void ResetHydraulicDischarge(); void CalculateHydraulicDischarge(const ProcessInfo& rCurrentProcessInfo); - void CalculateSoilGamma(ElementVariables &rVariables); - virtual void CalculateSoilDensity(ElementVariables &rVariables); - - virtual void CalculateAndAddGeometricStiffnessMatrix( MatrixType& rLeftHandSideMatrix, - ElementVariables& rVariables, - unsigned int GPoint ); + void CalculateSoilGamma(ElementVariables& rVariables); + virtual void CalculateSoilDensity(ElementVariables& rVariables); -///---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + virtual void CalculateAndAddGeometricStiffnessMatrix(MatrixType& rLeftHandSideMatrix, + ElementVariables& rVariables, + unsigned int GPoint); private: - /// Serialization - friend class Serializer; void save(Serializer& rSerializer) const override @@ -354,10 +341,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : KRATOS_SERIALIZE_LOAD_BASE_CLASS( rSerializer, Element ) } - // Private Operations - - template < class TValueType > - inline void ThreadSafeNodeWrite(NodeType& rNode, const Variable &Var, const TValueType Value) + template + inline void ThreadSafeNodeWrite(NodeType& rNode, const Variable& Var, const TValueType Value) { rNode.SetLock(); rNode.FastGetSolutionStepValue(Var) = Value; diff --git a/applications/GeoMechanicsApplication/tests/test_k0_procedure_process/test_k0_procedure_k0_nc_layers/mesh_stage1.mdpa b/applications/GeoMechanicsApplication/tests/test_k0_procedure_process/test_k0_procedure_k0_nc_layers/mesh_stage1.mdpa index a681c15644a7..37e06baa5746 100644 --- a/applications/GeoMechanicsApplication/tests/test_k0_procedure_process/test_k0_procedure_k0_nc_layers/mesh_stage1.mdpa +++ b/applications/GeoMechanicsApplication/tests/test_k0_procedure_process/test_k0_procedure_k0_nc_layers/mesh_stage1.mdpa @@ -59,7 +59,7 @@ Begin Nodes 53 10.0000000000 -100.0000000000 0.0000000000 End Nodes -Begin Elements SmallStrainUPwDiffOrderElement2D8N +Begin Elements UPwSmallStrainElement2D8N 1 3 51 53 48 46 52 50 47 49 2 3 46 48 43 41 47 45 42 44 3 3 41 43 38 36 42 40 37 39