Reorder DOFs of non-diff order U-Pw elements and conditions

This commit reorders the degrees of freedom (DOFs) of the non-diff order elements and conditions, such that they are consistent with the ordering applied to the diff order elements and conditions. In essence, this boils down to first enumerating all displacement DOFs, followed by enumerating the pore water pressure DOFs.

Changes in this commit include:
- Added utility functions that extract all degrees of freedom of U-Pw elements given the element's nodes and its dimension. These new functions are now used by the diff order elements (and conditions) as well as the non-diff order ones.
- Removed several duplicated or unused (utility) functions.
- Renamed a few utility functions for assembling matrices to make clear which part of the partitioned matrices they consider.
- Re-implemented the utility functions that assemble matrices and vectors (to avoid lots of raw loops).
- Replaced the use of matrix/vector assemble functions where a dimension of 0 was supplied to force the position of addition. Instead, we now use the matrix/vector addition operators, which is much clearer and cleaner.
- Renamed a few function parameters and template type names to be more accurate.
- Added a header file for common (type) aliases.
- Fixed several (minor) code smells, including the removal of commented out code.
- Extended the utility functions for setting up test models. Also added a few helper functions for carrying out common checks in tests.
    Added documentation which describes how the degrees of freedom of the U-Pw elements have been ordered.

applications/GeoMechanicsApplication/custom_conditions/Pw_normal_flux_condition.cpp

@@ -83,9 +83,7 @@ void PwNormalFluxCondition<TDim,TNumNodes>::
 {
     noalias(rVariables.PVector) = - rVariables.NormalFlux * rVariables.Np * rVariables.IntegrationCoefficient;
 
-    GeoElementUtilities::
-        AssemblePBlockVector< 0, TNumNodes >(rRightHandSideVector,
-                                                rVariables.PVector);
+    rRightHandSideVector += rVariables.PVector;
 }
 
 

applications/GeoMechanicsApplication/custom_conditions/T_microclimate_flux_condition.cpp

@@ -151,7 +151,7 @@ void GeoTMicroClimateFluxCondition<TDim, TNumNodes>::CalculateAndAddLHS(
     const auto flux_matrix = BoundedMatrix<double, TNumNodes, TNumNodes>{
         outer_prod(rN, element_prod(rN, rLeftHandSideFluxes)) * IntegrationCoefficient};
 
-    GeoElementUtilities::AssemblePBlockMatrix<0, TNumNodes>(rLeftHandSideMatrix, flux_matrix);
+    rLeftHandSideMatrix += flux_matrix;
 
     KRATOS_CATCH("")
 }
@@ -167,16 +167,11 @@ void GeoTMicroClimateFluxCondition<TDim, TNumNodes>::CalculateAndAddRHS(
 {
     auto temporary_matrix = BoundedMatrix<double, TNumNodes, TNumNodes>{
         outer_prod(rN, rN) * IntegrationCoefficient};
-    auto temporary_vector =
-        array_1d<double, TNumNodes>{prod(temporary_matrix, rRightHandSideFluxes)};
-    GeoElementUtilities::AssemblePBlockVector<0, TNumNodes>(
-        rRightHandSideVector, temporary_vector);
+    rRightHandSideVector += prod(temporary_matrix, rRightHandSideFluxes);
 
     temporary_matrix = BoundedMatrix<double, TNumNodes, TNumNodes>{
         outer_prod(rN, element_prod(rN, rLeftHandSideFluxes)) * IntegrationCoefficient};
-    temporary_vector = -prod(temporary_matrix, rNodalTemperatures);
-    GeoElementUtilities::AssemblePBlockVector<0, TNumNodes>(
-        rRightHandSideVector, temporary_vector);
+    rRightHandSideVector -= prod(temporary_matrix, rNodalTemperatures);
 }
 
 template <unsigned int TDim, unsigned int TNumNodes>

applications/GeoMechanicsApplication/custom_conditions/T_normal_flux_condition.cpp

@@ -80,9 +80,7 @@ void GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateRHS(Vector& rRightHandSi
                 r_integration_points[integration_point].Weight());
 
         // Contributions to the right hand side
-        auto normal_flux_on_DOF = normal_flux_on_integration_point * N * weighting_integration_coefficient;
-        GeoElementUtilities::AssemblePBlockVector<0, TNumNodes>(rRightHandSideVector,
-                                                                normal_flux_on_DOF);
+        rRightHandSideVector += (normal_flux_on_integration_point * N * weighting_integration_coefficient);
     }
 }
 

applications/GeoMechanicsApplication/custom_conditions/U_Pw_condition.cpp

@@ -127,14 +127,7 @@ void UPwCondition<TDim,TNumNodes>::
 template <unsigned int TDim, unsigned int TNumNodes>
 Condition::DofsVectorType UPwCondition<TDim, TNumNodes>::GetDofs() const
 {
-    auto result = Condition::DofsVectorType{};
-    for (const auto& r_node : GetGeometry()) {
-        result.push_back(r_node.pGetDof(DISPLACEMENT_X));
-        result.push_back(r_node.pGetDof(DISPLACEMENT_Y));
-        if constexpr (TDim == 3) result.push_back(r_node.pGetDof(DISPLACEMENT_Z));
-        result.push_back(r_node.pGetDof(WATER_PRESSURE));
-    }
-    return result;
+    return Geo::DofUtilities::ExtractUPwDofsFromNodes(GetGeometry(), TDim);
 }
 
 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

applications/GeoMechanicsApplication/custom_conditions/U_Pw_face_load_condition.cpp

@@ -15,6 +15,7 @@
 
 // Application includes
 #include "custom_conditions/U_Pw_face_load_condition.hpp"
+#include "custom_utilities/element_utilities.hpp"
 
 namespace Kratos
 {
@@ -71,7 +72,7 @@ void UPwFaceLoadCondition<TDim,TNumNodes>::
 
         //Contributions to the right hand side
         noalias(UVector) = prod(trans(Nu),TractionVector) * integration_coefficient;
-        ConditionUtilities::AssembleUBlockVector<TDim, TNumNodes>(rRightHandSideVector, UVector);
+        GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, UVector);
     }
 }
 

applications/GeoMechanicsApplication/custom_conditions/U_Pw_face_load_interface_condition.cpp

@@ -14,6 +14,7 @@
 
 // Application includes
 #include "custom_conditions/U_Pw_face_load_interface_condition.hpp"
+#include "custom_utilities/element_utilities.hpp"
 
 namespace Kratos
 {
@@ -119,7 +120,7 @@ void UPwFaceLoadInterfaceCondition<TDim,TNumNodes>::
 
         //Contributions to the right hand side
         noalias(UVector) = prod(trans(Nu),TractionVector) * integration_coefficient;
-        ConditionUtilities::AssembleUBlockVector<TDim, TNumNodes>(rRightHandSideVector,UVector);
+        GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector,UVector);
     }
 }
 

applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_face_load_condition.cpp

@@ -15,6 +15,7 @@
 
 // Application includes
 #include "custom_conditions/U_Pw_normal_face_load_condition.hpp"
+#include "custom_utilities/element_utilities.hpp"
 
 namespace Kratos
 {
@@ -68,7 +69,7 @@ void UPwNormalFaceLoadCondition<TDim,TNumNodes>::
 
         //Contributions to the right hand side
         noalias(UVector) = prod(trans(Nu),TractionVector) * IntegrationCoefficient;
-        ConditionUtilities::AssembleUBlockVector<TDim, TNumNodes>(rRightHandSideVector, UVector);
+        GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, UVector);
     }
 }
 

applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_flux_FIC_condition.cpp

@@ -199,11 +199,11 @@ template< unsigned int TDim, unsigned int TNumNodes >
 void UPwNormalFluxFICCondition<TDim,TNumNodes>::CalculateAndAddBoundaryMassMatrix(MatrixType& rLeftHandSideMatrix, NormalFluxVariables& rVariables,
                                                                                     NormalFluxFICVariables& rFICVariables)
 {
-    noalias(rFICVariables.PMatrix) = -rFICVariables.DtPressureCoefficient*rFICVariables.ElementLength*rFICVariables.BiotModulusInverse/6.0*
+    noalias(rFICVariables.PPMatrix) = -rFICVariables.DtPressureCoefficient*rFICVariables.ElementLength*rFICVariables.BiotModulusInverse/6.0*
                                         outer_prod(rVariables.Np,rVariables.Np)*rVariables.IntegrationCoefficient;
 
     //Distribute boundary mass matrix into the elemental matrix
-    GeoElementUtilities::AssemblePBlockMatrix< TDim, TNumNodes >(rLeftHandSideMatrix,rFICVariables.PMatrix);
+    GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix,rFICVariables.PPMatrix);
 }
 
 //----------------------------------------------------------------------------------------
@@ -222,14 +222,14 @@ template< unsigned int TDim, unsigned int TNumNodes >
 void UPwNormalFluxFICCondition<TDim,TNumNodes>::CalculateAndAddBoundaryMassFlow(VectorType& rRightHandSideVector, NormalFluxVariables& rVariables,
                                                                                     NormalFluxFICVariables& rFICVariables)
 {
-    noalias(rFICVariables.PMatrix) = rFICVariables.ElementLength*rFICVariables.BiotModulusInverse/6.0*
+    noalias(rFICVariables.PPMatrix) = rFICVariables.ElementLength*rFICVariables.BiotModulusInverse/6.0*
                                         outer_prod(rVariables.Np,rVariables.Np)*rVariables.IntegrationCoefficient;
 
 
-    noalias(rVariables.PVector) = prod(rFICVariables.PMatrix,rFICVariables.DtPressureVector);
+    noalias(rVariables.PVector) = prod(rFICVariables.PPMatrix,rFICVariables.DtPressureVector);
 
     //Distribute boundary mass flow vector into elemental vector
-    GeoElementUtilities::AssemblePBlockVector< TDim, TNumNodes >(rRightHandSideVector,rVariables.PVector);
+    GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector,rVariables.PVector);
 }
 
 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_flux_FIC_condition.hpp

@@ -79,7 +79,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION)
         double BiotModulusInverse;
 
         array_1d<double,TNumNodes> DtPressureVector;
-        BoundedMatrix<double,TNumNodes,TNumNodes> PMatrix;
+        BoundedMatrix<double,TNumNodes,TNumNodes> PPMatrix;
     };
 
     // Member Variables

applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_flux_condition.cpp

@@ -83,9 +83,7 @@ void UPwNormalFluxCondition<TDim,TNumNodes>::
 {
     noalias(rVariables.PVector) = - rVariables.NormalFlux * rVariables.Np * rVariables.IntegrationCoefficient;
 
-    GeoElementUtilities::
-        AssemblePBlockVector< TDim, TNumNodes >(rRightHandSideVector,
-                                                rVariables.PVector);
+    GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, rVariables.PVector);
 }
 
 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_flux_interface_condition.cpp

@@ -90,7 +90,7 @@ void UPwNormalFluxInterfaceCondition<TDim,TNumNodes>::
 
         //Contributions to the right hand side
         noalias(PVector) = -NormalFlux * Np * IntegrationCoefficient;
-        GeoElementUtilities::AssemblePBlockVector< TDim, TNumNodes >(rRightHandSideVector,PVector);
+        GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector,PVector);
     }
 }
 

applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.cpp

@@ -104,7 +104,7 @@ void UPwLysmerAbsorbingCondition<TDim, TNumNodes>::CalculateRightHandSide(Vector
     this->CalculateConditionStiffnessMatrix(stiffness_matrix, rCurrentProcessInfo);
 
     MatrixType global_stiffness_matrix = ZeroMatrix(CONDITION_SIZE, CONDITION_SIZE);
-    GeoElementUtilities::AssembleUBlockMatrix< TDim, TNumNodes >(global_stiffness_matrix, stiffness_matrix);
+    GeoElementUtilities::AssembleUUBlockMatrix(global_stiffness_matrix, stiffness_matrix);
 
     this->CalculateAndAddRHS(rRightHandSideVector, global_stiffness_matrix);
 }
@@ -420,15 +420,13 @@ CalculateAndAddRHS(VectorType& rRightHandSideVector, const MatrixType& rStiffnes
 
 template< unsigned int TDim, unsigned int TNumNodes >
 void UPwLysmerAbsorbingCondition<TDim, TNumNodes>::
-AddLHS(MatrixType& rLeftHandSideMatrix, const ElementMatrixType& rUMatrix)
+AddLHS(MatrixType& rLeftHandSideMatrix, const ElementMatrixType& rUUMatrix)
 {
 	// assemble left hand side vector
     rLeftHandSideMatrix = ZeroMatrix(CONDITION_SIZE, CONDITION_SIZE);
 
     //Adding contribution to left hand side
-    GeoElementUtilities::
-        AssembleUBlockMatrix< TDim, TNumNodes >(rLeftHandSideMatrix,
-            rUMatrix);
+    GeoElementUtilities::AssembleUUBlockMatrix(rLeftHandSideMatrix, rUUMatrix);
 }
 
 template< unsigned int TDim, unsigned int TNumNodes >

applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.hpp

@@ -131,9 +131,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwLysmerAbsorbingCondition : public
     /**
     * @brief Adds the condition matrix to the global left hand side
     * @param rLeftHandSideMatrix Global Left hand side
-    * @param rUMatrix LHS displacement matrix of the current condition
+    * @param rUUMatrix LHS displacement matrix of the current condition
     */
-    void AddLHS(MatrixType& rLeftHandSideMatrix, const ElementMatrixType& rUMatrix);
+    void AddLHS(MatrixType& rLeftHandSideMatrix, const ElementMatrixType& rUUMatrix);
 
     /**
     * @brief Calculates and adds terms to the RHS

applications/GeoMechanicsApplication/custom_elements/README.md

@@ -1,3 +1,15 @@
+# U-Pw Continuum Elements
+
+For simulating soil behavior, the geomechanics application offers continuum elements that take into account the displacement field ($u$) as well as the pore water pressure field ($p_w$). In fact, the two fields can be coupled.
+
+Depending on your needs, you may opt for either so-called "diff order" elements or "non-diff order" elements. The former implies that the order of the interpolation polynomials that will be used for the pore water pressure field is one less than the order of the interpolation polynomials of the displacement field. The latter, on the other hand, assumes one and the same order for the interpolation polynomials of both fields.
+
+Regardless of the above distinction, the application orders the degrees of freedom (DOFs) in the same way. The first part relates to the displacement DOFs (grouped per node), followed by a second part that contains the pore water pressure DOFs. For instance, assuming a 6-noded triangular diff order element, the DOFs are ordered as follows:
+
+$$ u_{x,1}, u_{y,1}, u_{x,2}, u_{y,2}, \dots, u_{x,6}, u_{y,6}, p_1, p_2, p_3 $$
+
+where the directions of the displacements are denoted by a suffix (here either $x$ or $y$), and the node number is also added as a suffix. Note that in this case, there are only three pore water pressure DOFs (due to the diff order nature of the element). For a non-diff order quadratic element the above list of DOFs is appended by $p_4$, $p_5$, and $p_6$.
+
 # Transient Thermal Element
 
 ## Introduction

applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp

@@ -618,14 +618,8 @@ unsigned int UPwBaseElement<TDim, TNumNodes>::GetNumberOfDOF() const
 template <unsigned int TDim, unsigned int TNumNodes>
 Element::DofsVectorType UPwBaseElement<TDim, TNumNodes>::GetDofs() const
 {
-    auto result = Element::DofsVectorType{};
-    for (const auto& r_node : this->GetGeometry()) {
-        result.push_back(r_node.pGetDof(DISPLACEMENT_X));
-        result.push_back(r_node.pGetDof(DISPLACEMENT_Y));
-        if constexpr (TDim == 3) result.push_back(r_node.pGetDof(DISPLACEMENT_Z));
-        result.push_back(r_node.pGetDof(WATER_PRESSURE));
-    }
-    return result;
+    return Geo::DofUtilities::ExtractUPwDofsFromNodes(this->GetGeometry(),
+                                                      this->GetGeometry().WorkingSpaceDimension());
 }
 
 template <unsigned int TDim, unsigned int TNumNodes>

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp

@@ -917,7 +917,7 @@ void UPwSmallStrainFICElement<2, 4>::CalculateAndAddStrainGradientMatrix(MatrixT
         prod(rVariables.GradNpT, rFICVariables.StrainGradients) * rVariables.IntegrationCoefficient;
 
     // Distribute strain gradient matrix into the elemental matrix
-    GeoElementUtilities::AssemblePUBlockMatrix<2, 4>(rLeftHandSideMatrix, rVariables.PUMatrix);
+    GeoElementUtilities::AssemblePUBlockMatrix(rLeftHandSideMatrix, rVariables.PUMatrix);
 
     KRATOS_CATCH("")
 }
@@ -946,7 +946,7 @@ void UPwSmallStrainFICElement<3, 8>::CalculateAndAddStrainGradientMatrix(MatrixT
         prod(rVariables.GradNpT, rFICVariables.StrainGradients) * rVariables.IntegrationCoefficient;
 
     // Distribute strain gradient matrix into the elemental matrix
-    GeoElementUtilities::AssemblePUBlockMatrix<3, 8>(rLeftHandSideMatrix, rVariables.PUMatrix);
+    GeoElementUtilities::AssemblePUBlockMatrix(rLeftHandSideMatrix, rVariables.PUMatrix);
 
     KRATOS_CATCH("")
 }
@@ -970,7 +970,7 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAndAddDtStressGradientM
                                    rVariables.IntegrationCoefficient;
 
     // Distribute DtStressGradient Matrix into the elemental matrix
-    GeoElementUtilities::AssemblePUBlockMatrix<TDim, TNumNodes>(rLeftHandSideMatrix, rVariables.PUMatrix);
+    GeoElementUtilities::AssemblePUBlockMatrix(rLeftHandSideMatrix, rVariables.PUMatrix);
 
     KRATOS_CATCH("")
 }
@@ -1270,14 +1270,14 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAndAddPressureGradientM
     const double StabilizationParameter = rFICVariables.ElementLength * rFICVariables.ElementLength *
                                           SignBiotCoefficient / (8.0 * rFICVariables.ShearModulus);
 
-    noalias(rVariables.PMatrix) =
+    noalias(rVariables.PPMatrix) =
         rVariables.DtPressureCoefficient * StabilizationParameter *
         (SignBiotCoefficient - 2.0 * rFICVariables.ShearModulus * rVariables.BiotModulusInverse /
                                    (3.0 * SignBiotCoefficient)) *
         prod(rVariables.GradNpT, trans(rVariables.GradNpT)) * rVariables.IntegrationCoefficient;
 
     // Distribute pressure gradient block matrix into the elemental matrix
-    GeoElementUtilities::AssemblePBlockMatrix<TDim, TNumNodes>(rLeftHandSideMatrix, rVariables.PMatrix);
+    GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, rVariables.PPMatrix);
 
     KRATOS_CATCH("")
 }
@@ -1324,7 +1324,7 @@ void UPwSmallStrainFICElement<2, 4>::CalculateAndAddStrainGradientFlow(VectorTyp
     noalias(rVariables.PVector) = prod(rVariables.PUMatrix, rVariables.VelocityVector);
 
     // Distribute Strain Gradient vector into elemental vector
-    GeoElementUtilities::AssemblePBlockVector<2, 4>(rRightHandSideVector, rVariables.PVector);
+    GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, rVariables.PVector);
 
     KRATOS_CATCH("")
 }
@@ -1354,7 +1354,7 @@ void UPwSmallStrainFICElement<3, 8>::CalculateAndAddStrainGradientFlow(VectorTyp
     noalias(rVariables.PVector) = prod(rVariables.PUMatrix, rVariables.VelocityVector);
 
     // Distribute Strain Gradient vector into elemental vector
-    GeoElementUtilities::AssemblePBlockVector<3, 8>(rRightHandSideVector, rVariables.PVector);
+    GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, rVariables.PVector);
 
     KRATOS_CATCH("")
 }
@@ -1377,7 +1377,7 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAndAddDtStressGradientF
                                   rVariables.IntegrationCoefficient;
 
     // Distribute DtStressGradient block vector into elemental vector
-    GeoElementUtilities::AssemblePBlockVector<TDim, TNumNodes>(rRightHandSideVector, rVariables.PVector);
+    GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, rVariables.PVector);
 
     KRATOS_CATCH("")
 }
@@ -1430,16 +1430,16 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAndAddPressureGradientF
     double StabilizationParameter = rFICVariables.ElementLength * rFICVariables.ElementLength *
                                     SignBiotCoefficient / (8.0 * rFICVariables.ShearModulus);
 
-    noalias(rVariables.PMatrix) =
+    noalias(rVariables.PPMatrix) =
         StabilizationParameter *
         (SignBiotCoefficient - 2.0 * rFICVariables.ShearModulus * rVariables.BiotModulusInverse /
                                    (3.0 * SignBiotCoefficient)) *
         prod(rVariables.GradNpT, trans(rVariables.GradNpT)) * rVariables.IntegrationCoefficient;
 
-    noalias(rVariables.PVector) = -1.0 * prod(rVariables.PMatrix, rVariables.DtPressureVector);
+    noalias(rVariables.PVector) = -1.0 * prod(rVariables.PPMatrix, rVariables.DtPressureVector);
 
     // Distribute PressureGradient block vector into elemental vector
-    GeoElementUtilities::AssemblePBlockVector<TDim, TNumNodes>(rRightHandSideVector, rVariables.PVector);
+    GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, rVariables.PVector);
 
     KRATOS_CATCH("")
 }