[OptApp] Mass Response with shape sens with FD

applications/OptimizationApplication/custom_python/add_custom_response_utilities_to_python.cpp

@@ -35,7 +35,7 @@ void  AddCustomResponseUtilitiesToPython(pybind11::module& m)
     m.def_submodule("MassResponseUtils")
         .def("Check", &MassResponseUtils::Check)
         .def("CalculateValue", &MassResponseUtils::CalculateValue)
-        .def("CalculateGradient", &MassResponseUtils::CalculateGradient, py::arg("list_of_gradient_variables"), py::arg("list_of_gradient_required_model_parts"), py::arg("list_of_gradient_computed_model_parts"), py::arg("list_of_container_expressions"))
+        .def("CalculateGradient", &MassResponseUtils::CalculateGradient, py::arg("list_of_gradient_variables"), py::arg("list_of_gradient_required_model_parts"), py::arg("list_of_gradient_computed_model_parts"), py::arg("list_of_container_expressions"), py::arg("perturbation_size"))
         ;
 
     m.def_submodule("LinearStrainEnergyResponseUtils")

applications/OptimizationApplication/custom_utilities/response/mass_response_utils.cpp

@@ -120,7 +120,8 @@ void MassResponseUtils::CalculateGradient(
     const PhysicalFieldVariableTypes& rPhysicalVariable,
     ModelPart& rGradientRequiredModelPart,
     ModelPart& rGradientComputedModelPart,
-    std::vector<ContainerExpressionType>& rListOfContainerExpressions)
+    std::vector<ContainerExpressionType>& rListOfContainerExpressions,
+    const double PerturbationSize)
 {
     KRATOS_TRY
 
@@ -148,7 +149,7 @@ void MassResponseUtils::CalculateGradient(
             VariableUtils().SetNonHistoricalVariableToZero(SHAPE_SENSITIVITY, rGradientRequiredModelPart.Nodes());
 
             // computes density sensitivty and store it within each elements' properties
-            CalculateMassShapeGradient(rGradientComputedModelPart, SHAPE_SENSITIVITY);
+            CalculateMassShapeGradient(rGradientComputedModelPart, SHAPE_SENSITIVITY, PerturbationSize);
         } else {
             KRATOS_ERROR
                 << "Unsupported sensitivity w.r.t. " << pVariable->Name()
@@ -198,13 +199,12 @@ void MassResponseUtils::CalculateGradient(
 
 void MassResponseUtils::CalculateMassShapeGradient(
     ModelPart& rModelPart,
-    const Variable<array_1d<double, 3>>& rOutputGradientVariable)
+    const Variable<array_1d<double, 3>>& rOutputGradientVariable,
+    const double PerturbationSize)
 {
     KRATOS_TRY
 
-    if (rModelPart.NumberOfElements() == 0) {
-        return;
-    }
+    KRATOS_ERROR_IF(rModelPart.NumberOfElements() == 0) << rModelPart.FullName() << " does not contain any elements.\n";
 
     KRATOS_ERROR_IF_NOT(HasVariableInProperties(rModelPart, DENSITY))
         << "DENSITY is not found in element properties of " << rModelPart.FullName() << ".\n";
@@ -227,12 +227,13 @@ void MassResponseUtils::CalculateMassShapeGradient(
         get_cross_area = [](const ModelPart::ElementType& rElement) -> double { return 1.0; };
     }
 
-    using VolumeDerivativeMethodType = std::function<double(IndexType, IndexType, const GeometryType&)>;
+    using VolumeDerivativeMethodType = std::function<double(IndexType, IndexType, GeometryType&)>;
 
     VolumeDerivativeMethodType volume_derivative_method;
+    bool is_parallel_computation_allowed = true;
     switch (rModelPart.Elements().begin()->GetGeometry().GetGeometryType()) {
         case GeometryData::KratosGeometryType::Kratos_Line2D2:
-            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  const GeometryType& rGeometry) {
+            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  GeometryType& rGeometry) {
                 const double lx = rGeometry[0].X() - rGeometry[1].X();
                 const double lx_derivative = ((NodeIndex == 0) - (NodeIndex == 1)) * (DirectionIndex == 0);
 
@@ -246,7 +247,7 @@ void MassResponseUtils::CalculateMassShapeGradient(
             };
             break;
         case GeometryData::KratosGeometryType::Kratos_Line3D2:
-            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  const GeometryType& rGeometry) {
+            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  GeometryType& rGeometry) {
                 const double lx = rGeometry[0].X() - rGeometry[1].X();
                 const double lx_derivative = ((NodeIndex == 0) - (NodeIndex == 1)) * (DirectionIndex == 0);
 
@@ -263,52 +264,76 @@ void MassResponseUtils::CalculateMassShapeGradient(
             };
             break;
         case GeometryData::KratosGeometryType::Kratos_Triangle2D3:
-            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  const GeometryType& rGeometry) {
+            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  GeometryType& rGeometry) {
                 return 2.0 * ElementSizeCalculator<2, 3>::AverageElementSize(rGeometry) * ElementSizeCalculator<2, 3>::AverageElementSizeDerivative(NodeIndex, DirectionIndex, rGeometry);
             };
             break;
         case GeometryData::KratosGeometryType::Kratos_Quadrilateral2D4:
-            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  const GeometryType& rGeometry) {
+            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  GeometryType& rGeometry) {
                 return 2.0 * ElementSizeCalculator<2, 4>::AverageElementSize(rGeometry) * ElementSizeCalculator<2, 4>::AverageElementSizeDerivative(NodeIndex, DirectionIndex, rGeometry);
             };
             break;
         case GeometryData::KratosGeometryType::Kratos_Tetrahedra3D4:
-            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  const GeometryType& rGeometry) {
+            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  GeometryType& rGeometry) {
                 return 3.0 * std::pow(ElementSizeCalculator<3, 4>::AverageElementSize(rGeometry), 2) * ElementSizeCalculator<3, 4>::AverageElementSizeDerivative(NodeIndex, DirectionIndex, rGeometry);
             };
             break;
         case GeometryData::KratosGeometryType::Kratos_Prism3D6:
-            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  const GeometryType& rGeometry) {
+            volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  GeometryType& rGeometry) {
                 return 3.0 * std::pow(ElementSizeCalculator<3, 6>::AverageElementSize(rGeometry), 2) * ElementSizeCalculator<3, 6>::AverageElementSizeDerivative(NodeIndex, DirectionIndex, rGeometry);
             };
             break;
-        // Following is not consistent with the DomainSize calculation, hence commented out for the time being.
-        // case GeometryData::KratosGeometryType::Kratos_Hexahedra3D8:
-        //     volume_derivative_method = [](IndexType NodeIndex, IndexType DirectionIndex,  const GeometryType& rGeometry) {
-        //         return 3.0 * std::pow(ElementSizeCalculator<3, 8>::AverageElementSize(rGeometry), 2) * ElementSizeCalculator<3, 8>::AverageElementSizeDerivative(NodeIndex, DirectionIndex, rGeometry);
-        //     };
-        //     break;
         default:
-            KRATOS_ERROR << "Non supported geometry type for mass shape sensitivity calculation (CalculateMassShapeGradient())." << std::endl;
+            is_parallel_computation_allowed = false;
+            volume_derivative_method = [PerturbationSize](IndexType NodeIndex, IndexType DirectionIndex, GeometryType& rGeometry) {
+                auto& coordinates = rGeometry[NodeIndex].Coordinates();
+                coordinates[DirectionIndex] += PerturbationSize;
+                const double perturbed_domain_size = rGeometry.DomainSize();
+                coordinates[DirectionIndex] -= PerturbationSize;
+                return perturbed_domain_size;
+            };
+            break;
     }
 
-    block_for_each(rModelPart.Elements(), [&](auto& rElement){
-        auto& r_geometry = rElement.GetGeometry();
-        const IndexType dimension = r_geometry.WorkingSpaceDimension();
+    if (is_parallel_computation_allowed) {
+        block_for_each(rModelPart.Elements(), [&](auto& rElement){
+            auto& r_geometry = rElement.GetGeometry();
+            const IndexType dimension = r_geometry.WorkingSpaceDimension();
 
-        const double density = rElement.GetProperties()[DENSITY];
-        const double thickness = get_thickness(rElement);
-        const double cross_area = get_cross_area(rElement);
+            const double density = rElement.GetProperties()[DENSITY];
+            const double thickness = get_thickness(rElement);
+            const double cross_area = get_cross_area(rElement);
 
-        for (IndexType c = 0; c < r_geometry.PointsNumber(); ++c) {
-            auto& r_derivative_value = r_geometry[c].GetValue(rOutputGradientVariable);
+            for (IndexType c = 0; c < r_geometry.PointsNumber(); ++c) {
+                auto& r_derivative_value = r_geometry[c].GetValue(rOutputGradientVariable);
 
-            for (IndexType k = 0; k < dimension; ++k) {
-                const double derivative_value = volume_derivative_method(c, k, r_geometry) * thickness * density * cross_area;
-                AtomicAdd(r_derivative_value[k], derivative_value);
+                for (IndexType k = 0; k < dimension; ++k) {
+                    const double derivative_value = volume_derivative_method(c, k, r_geometry) * thickness * density * cross_area;
+                    AtomicAdd(r_derivative_value[k], derivative_value);
+                }
             }
-        }
-    });
+        });
+    } else {
+        for (auto& r_element : rModelPart.Elements()) {
+            auto& r_geometry = r_element.GetGeometry();
+            const IndexType dimension = r_geometry.WorkingSpaceDimension();
+
+            const double density = r_element.GetProperties()[DENSITY];
+            const double thickness = get_thickness(r_element);
+            const double cross_area = get_cross_area(r_element);
+            const double initial_domain_size = r_geometry.DomainSize();
+
+            for (IndexType c = 0; c < r_geometry.PointsNumber(); ++c) {
+                auto& r_derivative_value = r_geometry[c].GetValue(rOutputGradientVariable);
+
+                for (IndexType k = 0; k < dimension; ++k) {
+                    const double perturbed_domain_size = volume_derivative_method(c, k, r_geometry);
+                    const double domain_size_derivative = (perturbed_domain_size - initial_domain_size) * thickness * density * cross_area / PerturbationSize;
+                    r_derivative_value[k] += domain_size_derivative;
+                }
+            }
+        };
+    }
 
     rModelPart.GetCommunicator().AssembleNonHistoricalData(rOutputGradientVariable);
 

applications/OptimizationApplication/custom_utilities/response/mass_response_utils.h

@@ -56,7 +56,8 @@ class KRATOS_API(OPTIMIZATION_APPLICATION) MassResponseUtils
         const PhysicalFieldVariableTypes& rPhysicalVariable,
         ModelPart& rGradientRequiredModelPart,
         ModelPart& rGradientComputedModelPart,
-        std::vector<ContainerExpressionType>& rListOfContainerExpressions);
+        std::vector<ContainerExpressionType>& rListOfContainerExpressions,
+        const double PerturbationSize);
 
     ///@}
 private:
@@ -75,7 +76,8 @@ class KRATOS_API(OPTIMIZATION_APPLICATION) MassResponseUtils
 
     static void CalculateMassShapeGradient(
         ModelPart& rModelPart,
-        const Variable<array_1d<double, 3>>& rOutputGradientVariable);
+        const Variable<array_1d<double, 3>>& rOutputGradientVariable,
+        const double PerturbationSize);
 
     static void CalculateMassDensityGradient(
         ModelPart& rModelPart,

applications/OptimizationApplication/python_scripts/responses/mass_response_function.py

@@ -23,7 +23,8 @@ def __init__(self, name: str, model: Kratos.Model, parameters: Kratos.Parameters
         default_settings = Kratos.Parameters("""{
             "evaluated_model_part_names"     : [
                 "PLEASE_PROVIDE_A_MODEL_PART_NAME"
-            ]
+            ],
+            "perturbation_size": 1e-6
         }""")
         parameters.ValidateAndAssignDefaults(default_settings)
 
@@ -35,6 +36,7 @@ def __init__(self, name: str, model: Kratos.Model, parameters: Kratos.Parameters
 
         self.model_part_operation = ModelPartOperation(self.model, ModelPartOperation.OperationType.UNION, f"response_{self.GetName()}", evaluated_model_part_names, False)
         self.model_part: Optional[Kratos.ModelPart] = None
+        self.perturbation_size = parameters["perturbation_size"].GetDouble()
 
     def GetImplementedPhysicalKratosVariables(self) -> 'list[SupportedSensitivityFieldVariableTypes]':
         return [KratosOA.SHAPE, Kratos.DENSITY, Kratos.THICKNESS, KratosOA.CROSS_AREA]
@@ -74,7 +76,8 @@ def CalculateGradient(self, physical_variable_collective_expressions: 'dict[Supp
                 physical_variable,
                 merged_model_part,
                 intersected_model_part_map[physical_variable],
-                physical_variable_collective_expressions[physical_variable].GetContainerExpressions())
+                physical_variable_collective_expressions[physical_variable].GetContainerExpressions(),
+                self.perturbation_size)
 
     def __str__(self) -> str:
         return f"Response [type = {self.__class__.__name__}, name = {self.GetName()}, model part name = {self.model_part.FullName()}]"

applications/OptimizationApplication/tests/responses_tests/test_mass_response_function.py

@@ -6,6 +6,7 @@
 
 import KratosMultiphysics.KratosUnittest as kratos_unittest
 from KratosMultiphysics.OptimizationApplication.responses.mass_response_function import MassResponseFunction
+import KratosMultiphysics.StructuralMechanicsApplication
 
 class TestMassResponseFunctionBase(kratos_unittest.TestCase, ABC):
     @classmethod
@@ -283,6 +284,77 @@ def test_CalculateDensitySensitivity(self):
             1e-6,
             6)
 
+class TestMassResponseFunctionQuads(TestMassResponseFunctionBase):
+    @classmethod
+    def GetParameters(cls):
+        return Kratos.Parameters("""{
+            "evaluated_model_part_names": ["test"]
+        }""")
+
+    @classmethod
+    def CreateElements(cls):
+        cls.model_part.CreateNewNode(1, 0.0, 0.0, 0.0)
+        cls.model_part.CreateNewNode(2, 2.0, 0.0, 0.0)
+        cls.model_part.CreateNewNode(3, 2.0, 2.0, 0.0)
+        cls.model_part.CreateNewNode(4, 0.0, 2.0, 0.0)
+
+        properties = cls.model_part.CreateNewProperties(1)
+        properties[Kratos.DENSITY] = 2.0
+        properties[Kratos.THICKNESS] = 0.01
+        cls.model_part.CreateNewElement("ShellThickElement3D4N", 1, [1, 2, 3, 4], properties)
+
+    def test_CalculateValue(self):
+        v = 0.0
+        for element in self.model_part.Elements:
+            v += element.GetGeometry().DomainSize() * element.Properties[Kratos.DENSITY] * element.Properties[Kratos.THICKNESS]
+        self.assertAlmostEqual(self.ref_value, v, 12)
+
+    def test_CalculateShapeSensitivity(self):
+        sensitivity = KratosOA.CollectiveExpression([Kratos.Expression.NodalExpression(self.model_part)])
+        self.response_function.CalculateGradient({KratosOA.SHAPE: sensitivity})
+
+        # calculate nodal shape sensitivities
+        self._CheckSensitivity(
+            self.response_function,
+            self.model_part.Nodes,
+            lambda x: x.GetValue(Kratos.SHAPE_SENSITIVITY_X),
+            lambda x, y: self._UpdateNodalPositions(0, x, y),
+            sensitivity.GetContainerExpressions()[0].Evaluate()[:, 0],
+            1e-6,
+            4)
+
+        self._CheckSensitivity(
+            self.response_function,
+            self.model_part.Nodes,
+            lambda x: x.GetValue(Kratos.SHAPE_SENSITIVITY_Y),
+            lambda x, y: self._UpdateNodalPositions(1, x, y),
+            sensitivity.GetContainerExpressions()[0].Evaluate()[:, 1],
+            1e-6,
+            4)
+
+        self._CheckSensitivity(
+            self.response_function,
+            self.model_part.Nodes,
+            lambda x: x.GetValue(Kratos.SHAPE_SENSITIVITY_Z),
+            lambda x, y: self._UpdateNodalPositions(2, x, y),
+            sensitivity.GetContainerExpressions()[0].Evaluate()[:, 2],
+            1e-6,
+            4)
+
+    def test_CalculateDensitySensitivity(self):
+        sensitivity = KratosOA.CollectiveExpression([Kratos.Expression.ElementExpression(self.model_part)])
+        self.response_function.CalculateGradient({Kratos.DENSITY: sensitivity})
+
+        # calculate element density sensitivity
+        self._CheckSensitivity(
+            self.response_function,
+            self.model_part.Elements,
+            lambda x: x.Properties[KratosOA.DENSITY_SENSITIVITY],
+            lambda x, y: self._UpdateProperties(Kratos.DENSITY, x, y),
+            sensitivity.GetContainerExpressions()[0].Evaluate(),
+            1e-6,
+            6)
+
 
 if __name__ == "__main__":
     Kratos.Tester.SetVerbosity(Kratos.Tester.Verbosity.PROGRESS)  # TESTS_OUTPUTS