[Core] ComputeNormalsOnSimplex - fix when there are other geometries (#…

…11876) * skipping non triangle conditions * adding test * fix for both dimensions * only run test in serial --------- Co-authored-by: Daniel Diez <[email protected]>
KratosMultiphysics · Dec 6, 2023 · bba1966 · bba1966
1 parent 45ff863
commit bba1966
Show file tree

Hide file tree

Showing 2 changed files with 47 additions and 8 deletions.
diff --git a/kratos/tests/test_normal_utils.py b/kratos/tests/test_normal_utils.py
@@ -130,6 +130,19 @@ def test_ComputeSimplexNormalModelPart(self):
 
             self.assertLess(CalculateNorm(normal - solution_normal), 0.15)
 
+    @KratosUnittest.skipIf(KratosMultiphysics.IsDistributedRun(), "This test is designed for serial runs only.")
+    def test_ComputeSimplexNormalModelPartWithLineCondition(self):
+        #Adding one line, to make sure it is getting ignored
+        self.model_part.CreateNewCondition("LineCondition3D2N", 1000, [1,2], self.model_part.GetProperties()[1])
+        KratosMultiphysics.NormalCalculationUtils().CalculateOnSimplex(self.model_part)
+        for node in self.model_part.GetSubModelPart("Skin_Part").Nodes:
+            normal = CalculateAnalyticalNormal(node)
+
+            solution_normal = node.GetSolutionStepValue(KratosMultiphysics.NORMAL)
+            solution_normal /= CalculateNorm(solution_normal)
+
+            self.assertLess(CalculateNorm(normal - solution_normal), 0.15)
+
     def test_ComputeSimplexNormalModelPartWithCustomVariable(self):
         ## Calculate the normals using NODAL_VAUX as custom variable
         KratosMultiphysics.NormalCalculationUtils().CalculateOnSimplex(

diff --git a/kratos/utilities/normal_calculation_utils.cpp b/kratos/utilities/normal_calculation_utils.cpp
@@ -29,6 +29,8 @@ void NormalCalculationUtils::CalculateNormalsInContainer(
     const Variable<array_1d<double,3>>& rNormalVariable
     )
 {
+    KRATOS_TRY;
+
     // Declare auxiliar coordinates
     Point::CoordinatesArrayType aux_coords;
 
@@ -50,6 +52,8 @@ void NormalCalculationUtils::CalculateNormalsInContainer(
         r_geometry.PointLocalCoordinates(aux_coords, r_geometry.Center());
         it_entity->SetValue(rNormalVariable, r_geometry.UnitNormal(aux_coords));
     }
+
+    KRATOS_CATCH("Error in CalculateNormalsInContainer");
 }
 
 /***********************************************************************************/
@@ -61,6 +65,8 @@ void NormalCalculationUtils::InitializeNormals(
     const Variable<array_1d<double,3>>& rNormalVariable
     )
 {
+    KRATOS_TRY
+
     // Resetting the normals
     const array_1d<double, 3> zero = ZeroVector(3);
 
@@ -94,6 +100,8 @@ void NormalCalculationUtils::InitializeNormals(
                 }
         });
     }
+
+    KRATOS_CATCH("Error in InitializeNormals");
 }
 
 /***********************************************************************************/
@@ -107,13 +115,17 @@ KRATOS_API(KRATOS_CORE) void NormalCalculationUtils::CalculateNormals<ModelPart:
     const Variable<array_1d<double,3>>& rNormalVariable
     )
 {
+    KRATOS_TRY
+
     if (CheckUseSimplex(rModelPart, EnforceGenericGeometryAlgorithm)) {
         const auto& r_process_info = rModelPart.GetProcessInfo();
         const SizeType dimension = r_process_info.GetValue(DOMAIN_SIZE);
         CalculateOnSimplex(rModelPart, dimension, rNormalVariable);
     } else {
         CalculateNormalsUsingGenericAlgorithm<ModelPart::ConditionsContainerType, true>(rModelPart, ConsiderUnitNormal, rNormalVariable);
     }
+
+    KRATOS_CATCH("Error in CalculateNormals");
 }
 
 /***********************************************************************************/
@@ -297,6 +309,8 @@ void NormalCalculationUtils::CalculateOnSimplex(
     const Variable<array_1d<double,3>>& rNormalVariable
     )
 {
+    KRATOS_TRY
+
     // Initialize the normals
     InitializeNormals<ModelPart::ConditionsContainerType,true>(rModelPart, rNormalVariable);
 
@@ -305,6 +319,8 @@ void NormalCalculationUtils::CalculateOnSimplex(
 
     // Synchronize the normal
     rModelPart.GetCommunicator().AssembleCurrentData(rNormalVariable);
+
+    KRATOS_CATCH("Error in CalculateOnSimplex");
 }
 
 /***********************************************************************************/
@@ -713,6 +729,8 @@ void NormalCalculationUtils::AuxiliaryCalculateOnSimplex(
     const NormalVariableType& rNormalVariable
     )
 {
+    KRATOS_TRY
+
     // Calculating the normals and storing on the conditions
     if (Dimension == 2) {
         block_for_each(rConditions, array_1d<double,3>(), [&rNormalVariable](Condition& rCondition, array_1d<double,3>& rAn){
@@ -736,18 +754,26 @@ void NormalCalculationUtils::AuxiliaryCalculateOnSimplex(
         });
     }
 
+    const auto condition_simplex_type = Dimension == 2
+            ? GeometryData::KratosGeometryType::Kratos_Line2D2
+            : GeometryData::KratosGeometryType::Kratos_Triangle3D3;
+
     // Adding the normals to the nodes
-    block_for_each(rConditions, [&rNormalVariable, this](Condition& rCondition) {
+    block_for_each(rConditions, [&rNormalVariable, &condition_simplex_type, this](Condition& rCondition) {
         auto& r_geometry = rCondition.GetGeometry();
-        const double coeff = 1.0 / r_geometry.size();
-        const auto& r_normal = rCondition.GetValue(rNormalVariable);
-        for (unsigned int i = 0; i < r_geometry.size(); ++i) {
-            auto& r_node = r_geometry[i];
-            r_node.SetLock();
-            noalias(GetNormalValue<TIsHistorical>(r_node, rNormalVariable)) += coeff * r_normal;
-            r_node.UnSetLock();
+        if (rCondition.GetGeometry().GetGeometryType() == condition_simplex_type) {
+            const double coeff = 1.0 / r_geometry.size();
+            const auto& r_normal = rCondition.GetValue(rNormalVariable);
+            for (unsigned int i = 0; i < r_geometry.size(); ++i) {
+                auto& r_node = r_geometry[i];
+                r_node.SetLock();
+                noalias(GetNormalValue<TIsHistorical>(r_node, rNormalVariable)) += coeff * r_normal;
+                r_node.UnSetLock();
+            }
         }
     });
+
+    KRATOS_CATCH("AuxiliaryCalculateOnSimplex");
 }
 
 // template instantiations