[SMApp] Reimplementing the thick shell element of 4 nodes

applications/ConstitutiveLawsApplication/custom_constitutive/structural_elements_constitutive_laws/thickness_integrated_isotropic_constitutive_law.cpp

@@ -533,7 +533,7 @@ void ThicknessIntegratedIsotropicConstitutiveLaw::CalculateMaterialResponseCauch
 
             if (flag_compute_constitutive_tensor) {
                 // membrane part
-                noalias(project(generalized_constitutive_matrix, range(0,3), range(0,3))) += weight * r_constitutive_matrix;
+                noalias(project(generalized_constitutive_matrix, range(0, 3), range(0, 3))) += weight * r_constitutive_matrix;
 
                 // bending part
                 noalias(project(generalized_constitutive_matrix, range(3, 6), range(3, 6))) += aux_weight2 * r_constitutive_matrix;

applications/StructuralMechanicsApplication/custom_constitutive/reissner_mindlin_shell_elastic_constitutive_law.h

@@ -149,17 +149,35 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) ReissnerMindlinShellElasticCo
     {
         double max_length = 0.0;
 
-        const auto& r_coord_1 = rGeometry[0].GetInitialPosition();
-        const auto& r_coord_2 = rGeometry[1].GetInitialPosition();
-        const auto& r_coord_3 = rGeometry[2].GetInitialPosition();
-
-        const double length_12 = norm_2(r_coord_2 - r_coord_1);
-        const double length_23 = norm_2(r_coord_3 - r_coord_2);
-        const double length_31 = norm_2(r_coord_1 - r_coord_3);
-
-        max_length = std::max(length_12, length_23);
-        max_length = std::max(max_length, length_31);
-
+        if (rGeometry.PointsNumber() == 3) {
+            const auto& r_coord_1 = rGeometry[0].GetInitialPosition();
+            const auto& r_coord_2 = rGeometry[1].GetInitialPosition();
+            const auto& r_coord_3 = rGeometry[2].GetInitialPosition();
+
+            const double length_12 = norm_2(r_coord_2 - r_coord_1);
+            const double length_23 = norm_2(r_coord_3 - r_coord_2);
+            const double length_31 = norm_2(r_coord_1 - r_coord_3);
+
+            max_length = std::max(length_12, length_23);
+            max_length = std::max(max_length, length_31);
+
+        } else if (rGeometry.PointsNumber() == 4) {
+            const auto& r_coord_1 = rGeometry[0].GetInitialPosition();
+            const auto& r_coord_2 = rGeometry[1].GetInitialPosition();
+            const auto& r_coord_3 = rGeometry[2].GetInitialPosition();
+            const auto& r_coord_4 = rGeometry[3].GetInitialPosition();
+
+            const double length_12 = norm_2(r_coord_2 - r_coord_1);
+            const double length_23 = norm_2(r_coord_3 - r_coord_2);
+            const double length_34 = norm_2(r_coord_4 - r_coord_3);
+            const double length_41 = norm_2(r_coord_1 - r_coord_4);
+
+            max_length = std::max(length_12, length_23);
+            max_length = std::max(max_length, length_34);
+            max_length = std::max(max_length, length_41);
+        } else {
+            KRATOS_ERROR << "GetMaxReferenceEdgeLength not implemented for geometries with " << rGeometry.PointsNumber() << " nodes." << std::endl;
+        }
         return max_length;
     }