[GeoMechanicsApplication] Remove swapping faces in condition process

applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_face_load_condition.cpp

@@ -98,7 +98,7 @@ void UPwNormalFaceLoadCondition<TDim, TNumNodes>::CalculateTractionVector(
     const unsigned int& GPoint)
 {
     Vector NormalVector = ZeroVector(TDim);
-    const double NormalStress = MathUtils<>::Dot(row(NContainer, GPoint), Variables.NormalStressVector);
+    double NormalStress = MathUtils<>::Dot(row(NContainer, GPoint), Variables.NormalStressVector);
 
     if constexpr (TDim == 2) {
         const double TangentialStress = MathUtils<>::Dot(row(NContainer, GPoint), Variables.TangentialStressVector);
@@ -107,6 +107,10 @@ void UPwNormalFaceLoadCondition<TDim, TNumNodes>::CalculateTractionVector(
         rTractionVector[1] = NormalStress * NormalVector[0] + TangentialStress * NormalVector[1];
     }
     else if constexpr (TDim == 3) {
+        // Since the normal vector is pointing outwards for the 3D conditions, the normal stress
+        // should switch sign, such that positive normal contact stress is defined inwards.
+        NormalStress *= -1;
+
         MathUtils<double>::CrossProduct(NormalVector, column(Jacobian, 0), column(Jacobian, 1));
         rTractionVector = NormalStress * NormalVector;
     }

applications/GeoMechanicsApplication/custom_conditions/surface_normal_load_3D_diff_order_condition.cpp

@@ -15,6 +15,7 @@
 
 // Project includes
 #include "custom_conditions/surface_normal_load_3D_diff_order_condition.hpp"
+#include "custom_utilities/math_utilities.hpp"
 
 namespace Kratos
 {
@@ -36,37 +37,28 @@ Condition::Pointer SurfaceNormalLoad3DDiffOrderCondition::Create(IndexType NewId
     return Condition::Pointer(new SurfaceNormalLoad3DDiffOrderCondition(NewId, GetGeometry().Create(ThisNodes), pProperties));
 }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-void SurfaceNormalLoad3DDiffOrderCondition::
-CalculateConditionVector(ConditionVariables& rVariables, unsigned int PointNumber)
+void SurfaceNormalLoad3DDiffOrderCondition::CalculateConditionVector(ConditionVariables& rVariables,
+                                                                     unsigned int PointNumber)
 {
     KRATOS_TRY
 
-    double NormalVector[3];
-
-    NormalVector[0] = rVariables.JContainer[PointNumber](1,0) * rVariables.JContainer[PointNumber](2,1) -
-                      rVariables.JContainer[PointNumber](2,0) * rVariables.JContainer[PointNumber](1,1);
-
-    NormalVector[1] = rVariables.JContainer[PointNumber](2,0) * rVariables.JContainer[PointNumber](0,1) -
-                      rVariables.JContainer[PointNumber](0,0) * rVariables.JContainer[PointNumber](2,1);
+    Vector normal_vector(3);
+    MathUtils<double>::CrossProduct(normal_vector, column(rVariables.JContainer[PointNumber], 0),
+                                    column(rVariables.JContainer[PointNumber], 1));
 
-    NormalVector[2] = rVariables.JContainer[PointNumber](0,0) * rVariables.JContainer[PointNumber](1,1) -
-                      rVariables.JContainer[PointNumber](1,0) * rVariables.JContainer[PointNumber](0,1);
-
-    const GeometryType& rGeom = GetGeometry();
-    const SizeType NumUNodes = rGeom.PointsNumber();
-    double NormalStress = 0.0;
-    rVariables.ConditionVector.resize(3,false);
+    const auto& r_geometry                   = GetGeometry();
 
-    for ( SizeType i = 0; i < NumUNodes; ++i ) {
-        NormalStress += rVariables.Nu[i]*rGeom[i].FastGetSolutionStepValue(NORMAL_CONTACT_STRESS);
-    }
+    Vector normal_stresses(r_geometry.PointsNumber());
+    std::transform(r_geometry.begin(), r_geometry.end(), normal_stresses.begin(), [](const auto& r_node) {
+        return r_node.FastGetSolutionStepValue(NORMAL_CONTACT_STRESS);
+    });
 
-    rVariables.ConditionVector[0] = NormalStress * NormalVector[0];
-    rVariables.ConditionVector[1] = NormalStress * NormalVector[1];
-    rVariables.ConditionVector[2] = NormalStress * NormalVector[2];
+    // Since the normal vector is pointing outwards for the 3D conditions, the normal stress
+    // should switch sign, such that positive normal contact stress is defined inwards.
+    const double normal_stress = -1 * MathUtils<>::Dot(rVariables.Nu, normal_stresses);
+    rVariables.ConditionVector = normal_stress * normal_vector;
 
-    KRATOS_CATCH( "" )
+    KRATOS_CATCH("")
 }
 
 //----------------------------------------------------------------------------------------
@@ -92,7 +84,7 @@ void SurfaceNormalLoad3DDiffOrderCondition::
 
     for ( SizeType i = 0; i < NumUNodes; ++i ) {
         Index = i * 3;
-        
+
         rRightHandSideVector[Index]   += rVariables.Nu[i] * rVariables.ConditionVector[0] * rVariables.IntegrationCoefficient;
         rRightHandSideVector[Index+1] += rVariables.Nu[i] * rVariables.ConditionVector[1] * rVariables.IntegrationCoefficient;
         rRightHandSideVector[Index+2] += rVariables.Nu[i] * rVariables.ConditionVector[2] * rVariables.IntegrationCoefficient;

applications/GeoMechanicsApplication/custom_processes/find_neighbour_elements_of_conditions_process.cpp

@@ -36,30 +36,6 @@ void FindNeighbourElementsOfConditionsProcess::Execute()
             Ids[i] = rGeometry[i].Id();
         }
 
-        if (rGeometry.GetGeometryType() == GeometryData::KratosGeometryType::Kratos_Triangle3D3) {
-            // reverse ordering to be consistent with face ordering
-            std::swap(Ids[1], Ids[2]);
-        }
-
-        if (rGeometry.GetGeometryType() == GeometryData::KratosGeometryType::Kratos_Triangle3D6) {
-            // reverse ordering to be consistent with face ordering
-            std::swap(Ids[1], Ids[2]);
-            std::swap(Ids[3], Ids[5]);
-        }
-
-        if (rGeometry.GetGeometryType() == GeometryData::KratosGeometryType::Kratos_Quadrilateral3D4) {
-            // reverse ordering to be consistent with face ordering
-            std::swap(Ids[0], Ids[3]);
-            std::swap(Ids[1], Ids[2]);
-        }
-
-        if (rGeometry.GetGeometryType() == GeometryData::KratosGeometryType::Kratos_Quadrilateral3D8) {
-            // reverse ordering to be consistent with face ordering
-            std::swap(Ids[0], Ids[3]);
-            std::swap(Ids[1], Ids[2]);
-            std::swap(Ids[4], Ids[6]);
-        }
-
         // adds to the map
         FacesMap.insert( hashmap::value_type(Ids, std::vector<Condition::Pointer>({*itCond.base()})) );
 

applications/GeoMechanicsApplication/tests/hexa_8n_normal_load/MaterialParameters.json

@@ -0,0 +1,29 @@
+{
+   "properties": [{
+      "model_part_name":         "PorousDomain.Soil",
+      "properties_id":           1,
+      "Material": {
+          "constitutive_law": {
+              "name"             :  "LinearElastic3DLaw" 
+          },
+          "Variables": {
+              "IGNORE_UNDRAINED"         :  true,
+              "YOUNG_MODULUS"            :  3.0e7,
+              "POISSON_RATIO"            :  0,
+              "DENSITY_SOLID"            :  0,
+              "DENSITY_WATER"            :  1.0e3,
+              "POROSITY"                 :  0.3,
+              "BULK_MODULUS_SOLID"       :  1.0e16,
+              "BULK_MODULUS_FLUID"       :  2.0e-30,
+              "PERMEABILITY_XX"          :  4.5e-30,
+              "PERMEABILITY_YY"          :  4.5e-30,
+              "PERMEABILITY_ZZ"          :  4.5e-30,
+              "PERMEABILITY_XY"          :  0.0,
+              "PERMEABILITY_YZ"          :  0.0,
+              "PERMEABILITY_ZX"          :  0.0,
+              "DYNAMIC_VISCOSITY"        :  1.0e-3
+         },
+         "Tables": {}
+      }
+   }]
+}