[GeoMechanicsApplication] Move all non-template code to .cpp files and make sure any .h file does not contain any implementations (elements)

applications/GeoMechanicsApplication/custom_elements/Pw_element.hpp

@@ -0,0 +1,235 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Mohamed Nabi
+//                   John van Esch
+//                   Richard Faasse
+//                   Gennady Markelov
+//
+
+#pragma once
+
+#include "calculation_contribution.h"
+#include "compressibility_calculator.hpp"
+#include "custom_retention/retention_law_factory.h"
+#include "custom_utilities/check_utilities.h"
+#include "custom_utilities/constitutive_law_utilities.h"
+#include "custom_utilities/dof_utilities.h"
+#include "custom_utilities/element_utilities.hpp"
+#include "custom_utilities/hydraulic_discharge.h"
+#include "custom_utilities/transport_equation_utilities.hpp"
+#include "custom_utilities/variables_utilities.hpp"
+#include "filter_compressibility_calculator.hpp"
+#include "fluid_body_flow_calculator.hpp"
+#include "geo_mechanics_application_variables.h"
+#include "includes/cfd_variables.h"
+#include "includes/element.h"
+#include "includes/serializer.h"
+#include "integration_coefficients_calculator.hpp"
+#include "permeability_calculator.hpp"
+
+#include <optional>
+
+namespace Kratos
+{
+
+template <unsigned int TDim, unsigned int TNumNodes>
+class KRATOS_API(GEO_MECHANICS_APPLICATION) PwElement : public Element
+{
+public:
+    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(PwElement);
+
+    explicit PwElement(IndexType NewId = 0);
+
+    PwElement(IndexType                                       NewId,
+              const GeometryType::Pointer&                    pGeometry,
+              const std::vector<CalculationContribution>&     rContributions,
+              std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier);
+
+    PwElement(IndexType                                       NewId,
+              const GeometryType::Pointer&                    pGeometry,
+              const PropertiesType::Pointer&                  pProperties,
+              const std::vector<CalculationContribution>&     rContributions,
+              std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier);
+
+    ~PwElement() override                      = default;
+    PwElement(const PwElement&)                = delete;
+    PwElement& operator=(const PwElement&)     = delete;
+    PwElement(PwElement&&) noexcept            = delete;
+    PwElement& operator=(PwElement&&) noexcept = delete;
+
+    Element::Pointer Create(IndexType               NewId,
+                            const NodesArrayType&   rThisNodes,
+                            PropertiesType::Pointer pProperties) const override;
+
+    Element::Pointer Create(IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const override;
+
+    void GetDofList(DofsVectorType& rElementalDofList, const ProcessInfo&) const override;
+
+    void EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo&) const override;
+
+    void Initialize(const ProcessInfo&) override;
+
+    void InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
+
+    void FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
+
+    using Element::CalculateOnIntegrationPoints;
+
+    void CalculateOnIntegrationPoints(const Variable<double>& rVariable,
+                                      std::vector<double>&    rOutput,
+                                      const ProcessInfo&      rCurrentProcessInfo) override;
+
+    void CalculateOnIntegrationPoints(const Variable<array_1d<double, 3>>& rVariable,
+                                      std::vector<array_1d<double, 3>>&    rOutput,
+                                      const ProcessInfo& rCurrentProcessInfo) override;
+
+    void CalculateLocalSystem(MatrixType&        rLeftHandSideMatrix,
+                              VectorType&        rRightHandSideVector,
+                              const ProcessInfo& rCurrentProcessInfo) override;
+
+    void CalculateRightHandSide(VectorType& rRightHandSideVector, const ProcessInfo& rCurrentProcessInfo) override;
+
+    void CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, const ProcessInfo& rCurrentProcessInfo) override;
+
+    GeometryData::IntegrationMethod GetIntegrationMethod() const override;
+
+    int Check(const ProcessInfo& rCurrentProcessInfo) const override;
+
+    void CalculateOnIntegrationPoints(const Variable<Matrix>& rVariable,
+                                      std::vector<Matrix>&    rOutput,
+                                      const ProcessInfo&      rCurrentProcessInfo) override;
+
+private:
+    std::vector<CalculationContribution> mContributions;
+    IntegrationCoefficientsCalculator    mIntegrationCoefficientsCalculator;
+    std::vector<RetentionLaw::Pointer>   mRetentionLawVector;
+    Vector                               mIntegrationCoefficients;
+    Matrix                               mNContainer;
+    Vector                               mDetJCcontainer;
+    std::vector<double>                  mFluidPressures;
+
+    std::vector<double> CalculateIntegrationCoefficients(const Vector& rDetJs) const;
+
+    std::vector<Vector> CalculateProjectedGravityAtIntegrationPoints(const Matrix& rNContainer) const;
+
+    std::unique_ptr<IntegrationCoefficientModifier> CloneIntegrationCoefficientModifier() const;
+
+    std::unique_ptr<ContributionCalculator<TNumNodes>> CreateCalculator(const CalculationContribution& rContribution,
+                                                                        const ProcessInfo& rCurrentProcessInfo);
+
+    void CachingDataForCalculator();
+
+    typename CompressibilityCalculator<TNumNodes>::InputProvider CreateCompressibilityInputProvider(
+        const ProcessInfo& rCurrentProcessInfo);
+
+    typename FilterCompressibilityCalculator<TNumNodes>::InputProvider CreateFilterCompressibilityInputProvider(
+        const ProcessInfo& rCurrentProcessInfo);
+
+    typename PermeabilityCalculator<TNumNodes>::InputProvider CreatePermeabilityInputProvider();
+
+    typename FluidBodyFlowCalculator<TNumNodes>::InputProvider CreateFluidBodyFlowInputProvider();
+
+    auto MakePropertiesGetter()
+    {
+        return [this]() -> const Properties& { return GetProperties(); };
+    }
+
+    auto MakeRetentionLawsGetter()
+    {
+        return [this]() -> const std::vector<RetentionLaw::Pointer>& { return mRetentionLawVector; };
+    }
+
+    auto GetNContainer()
+    {
+        return [this]() -> const Matrix& { return mNContainer; };
+    }
+
+    Matrix CalculateNContainer()
+    {
+        return GetGeometry().ShapeFunctionsValues(GetIntegrationMethod());
+    }
+
+    auto GetIntegrationCoefficients()
+    {
+        return [this]() -> const Vector& { return mIntegrationCoefficients; };
+    }
+
+    Vector CalculateIntegrationCoefficients()
+    {
+        GetGeometry().DeterminantOfJacobian(mDetJCcontainer, this->GetIntegrationMethod());
+        return mIntegrationCoefficientsCalculator.Run<Vector>(
+            GetGeometry().IntegrationPoints(GetIntegrationMethod()), mDetJCcontainer, this);
+    }
+
+    auto GetFluidPressures()
+    {
+        return [this]() -> const std::vector<double>& { return mFluidPressures; };
+    }
+
+    std::vector<double> CalculateFluidPressure()
+    {
+        return GeoTransportEquationUtilities::CalculateFluidPressures(
+            mNContainer, VariablesUtilities::GetNodalValuesOf<TNumNodes>(WATER_PRESSURE, this->GetGeometry()));
+    }
+
+    auto MakeProjectedGravityForIntegrationPointsGetter() const
+    {
+        return [this]() -> std::vector<Vector> {
+            return CalculateProjectedGravityAtIntegrationPoints(mNContainer);
+        };
+    }
+
+    static auto MakeMatrixScalarFactorGetter(const ProcessInfo& rCurrentProcessInfo)
+    {
+        return [&rCurrentProcessInfo]() { return rCurrentProcessInfo[DT_PRESSURE_COEFFICIENT]; };
+    }
+
+    auto MakeNodalVariableGetter() const
+    {
+        return [this](const Variable<double>& rVariable) -> Vector {
+            return VariablesUtilities::GetNodalValuesOf<TNumNodes>(rVariable, this->GetGeometry());
+        };
+    }
+
+    auto MakeShapeFunctionLocalGradientsGetter()
+    {
+        return [this]() {
+            GeometryType::ShapeFunctionsGradientsType dN_dX_container;
+            if (GetGeometry().LocalSpaceDimension() == 1) {
+                dN_dX_container = GetGeometry().ShapeFunctionsLocalGradients(this->GetIntegrationMethod());
+                std::transform(dN_dX_container.begin(), dN_dX_container.end(),
+                               mDetJCcontainer.begin(), dN_dX_container.begin(), std::divides<>());
+            } else {
+                GetGeometry().ShapeFunctionsIntegrationPointsGradients(
+                    dN_dX_container, mDetJCcontainer, this->GetIntegrationMethod());
+            }
+
+            return dN_dX_container;
+        };
+    }
+
+    auto MakeLocalSpaceDimensionGetter() const
+    {
+        return [this]() -> std::size_t { return this->GetGeometry().LocalSpaceDimension(); };
+    }
+
+    [[nodiscard]] DofsVectorType GetDofs() const
+    {
+        return Geo::DofUtilities::ExtractDofsFromNodes(GetGeometry(), WATER_PRESSURE);
+    }
+
+    friend class Serializer;
+
+    void save(Serializer& rSerializer) const override;
+
+    void load(Serializer& rSerializer) override;
+};
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp

@@ -12,7 +12,7 @@
 //
 
 // Application includes
-#include "custom_elements/U_Pw_base_element.hpp"
+#include "custom_elements/U_Pw_base_element.h"
 #include "custom_utilities/check_utilities.h"
 #include "custom_utilities/dof_utilities.h"
 #include "custom_utilities/equation_of_motion_utilities.h"
@@ -21,6 +21,38 @@
 
 namespace Kratos
 {
+UPwBaseElement::UPwBaseElement(IndexType                                       NewId,
+                               const NodesArrayType&                           ThisNodes,
+                               std::unique_ptr<StressStatePolicy>              pStressStatePolicy,
+                               std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier)
+    : Element(NewId, ThisNodes),
+      mpStressStatePolicy{std::move(pStressStatePolicy)},
+      mIntegrationCoefficientsCalculator{std::move(pCoefficientModifier)}
+{
+}
+
+UPwBaseElement::UPwBaseElement(IndexType                                       NewId,
+                               GeometryType::Pointer                           pGeometry,
+                               std::unique_ptr<StressStatePolicy>              pStressStatePolicy,
+                               std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier)
+    : Element(NewId, pGeometry),
+      mpStressStatePolicy{std::move(pStressStatePolicy)},
+      mIntegrationCoefficientsCalculator{std::move(pCoefficientModifier)}
+{
+}
+
+UPwBaseElement::UPwBaseElement(IndexType                                       NewId,
+                               GeometryType::Pointer                           pGeometry,
+                               PropertiesType::Pointer                         pProperties,
+                               std::unique_ptr<StressStatePolicy>              pStressStatePolicy,
+                               std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier)
+    : Element(NewId, pGeometry, pProperties),
+      mpStressStatePolicy{std::move(pStressStatePolicy)},
+      mIntegrationCoefficientsCalculator{std::move(pCoefficientModifier)}
+{
+    // this is needed for interface elements
+    mThisIntegrationMethod = this->GetIntegrationMethod();
+}
 
 int UPwBaseElement::Check(const ProcessInfo& rCurrentProcessInfo) const
 {
@@ -438,4 +470,17 @@ std::unique_ptr<IntegrationCoefficientModifier> UPwBaseElement::CloneIntegration
     return mIntegrationCoefficientsCalculator.CloneModifier();
 }
 
+std::string UPwBaseElement::Info() const
+{
+    const std::string constitutive_info =
+        !mConstitutiveLawVector.empty() ? mConstitutiveLawVector[0]->Info() : "not defined";
+
+    std::ostringstream oss;
+    oss << "U-Pw Base class Element #" << Id() << "\nConstitutive law: " << constitutive_info;
+
+    return oss.str();
+}
+
+void UPwBaseElement::PrintInfo(std::ostream& rOStream) const { rOStream << Info(); }
+
 } // Namespace Kratos

applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.hpp → applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.h

@@ -25,7 +25,7 @@
 
 // Application includes
 #include "custom_utilities/element_utilities.hpp"
-#include "integration_coefficients_calculator.h"
+#include "integration_coefficients_calculator.hpp"
 #include "stress_state_policy.h"
 
 namespace Kratos
@@ -44,37 +44,18 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element
     UPwBaseElement(IndexType                                       NewId,
                    const NodesArrayType&                           ThisNodes,
                    std::unique_ptr<StressStatePolicy>              pStressStatePolicy,
-                   std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier = nullptr)
-        : Element(NewId, ThisNodes),
-          mpStressStatePolicy{std::move(pStressStatePolicy)},
-          mIntegrationCoefficientsCalculator{std::move(pCoefficientModifier)}
-    {
-    }
-
-    /// Constructor using Geometry
+                   std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier = nullptr);
+
     UPwBaseElement(IndexType                                       NewId,
                    GeometryType::Pointer                           pGeometry,
                    std::unique_ptr<StressStatePolicy>              pStressStatePolicy,
-                   std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier = nullptr)
-        : Element(NewId, pGeometry),
-          mpStressStatePolicy{std::move(pStressStatePolicy)},
-          mIntegrationCoefficientsCalculator{std::move(pCoefficientModifier)}
-    {
-    }
-
-    /// Constructor using Properties
+                   std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier = nullptr);
+
     UPwBaseElement(IndexType                                       NewId,
                    GeometryType::Pointer                           pGeometry,
                    PropertiesType::Pointer                         pProperties,
                    std::unique_ptr<StressStatePolicy>              pStressStatePolicy,
-                   std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier = nullptr)
-        : Element(NewId, pGeometry, pProperties),
-          mpStressStatePolicy{std::move(pStressStatePolicy)},
-          mIntegrationCoefficientsCalculator{std::move(pCoefficientModifier)}
-    {
-        // this is needed for interface elements
-        mThisIntegrationMethod = this->GetIntegrationMethod();
-    }
+                   std::unique_ptr<IntegrationCoefficientModifier> pCoefficientModifier = nullptr);
 
     ~UPwBaseElement() override                           = default;
     UPwBaseElement(const UPwBaseElement&)                = delete;
@@ -148,14 +129,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element
 
     using Element::SetValuesOnIntegrationPoints;
 
-    std::string Info() const override
-    {
-        const std::string constitutive_info =
-            !mConstitutiveLawVector.empty() ? mConstitutiveLawVector[0]->Info() : "not defined";
-        return "U-Pw Base class Element #" + std::to_string(Id()) + "\nConstitutive law: " + constitutive_info;
-    }
+    std::string Info() const override;
 
-    void PrintInfo(std::ostream& rOStream) const override { rOStream << Info(); }
+    void PrintInfo(std::ostream& rOStream) const override;
 
 protected:
     /// Member Variables