Merge pull request #11450 from KratosMultiphysics/optapp/geometric_re…

…sponses

[OptApp] Adding geometric centroid deviation response function

applications/OptimizationApplication/python_scripts/responses/geometric_centroid_deviation_response_function.py

@@ -0,0 +1,105 @@
+from typing import Optional
+
+import KratosMultiphysics as Kratos
+import KratosMultiphysics.OptimizationApplication as KratosOA
+from KratosMultiphysics.OptimizationApplication.responses.response_function import ResponseFunction
+from KratosMultiphysics.OptimizationApplication.utilities.union_utilities import SupportedSensitivityFieldVariableTypes
+from KratosMultiphysics.OptimizationApplication.utilities.model_part_utilities import ModelPartOperation
+from KratosMultiphysics.OptimizationApplication.utilities.model_part_utilities import ModelPartUtilities
+
+def Factory(model: Kratos.Model, parameters: Kratos.Parameters, _) -> ResponseFunction:
+    if not parameters.Has("name"):
+        raise RuntimeError(f"GeometricCentroidDeviationResponseFunction instantiation requires a \"name\" in parameters [ parameters = {parameters}].")
+    if not parameters.Has("settings"):
+        raise RuntimeError(f"GeometricCentroidDeviationResponseFunction instantiation requires a \"settings\" in parameters [ parameters = {parameters}].")
+
+    return GeometricCentroidDeviationResponseFunction(parameters["name"].GetString(), model, parameters["settings"])
+
+class GeometricCentroidDeviationResponseFunction(ResponseFunction):
+    def __init__(self, name: str, model: Kratos.Model, parameters: Kratos.Parameters):
+        super().__init__(name)
+
+        default_settings = Kratos.Parameters("""{
+            "evaluated_model_part_names"     : [
+                "PLEASE_PROVIDE_A_MODEL_PART_NAME"
+            ]
+        }""")
+        parameters.ValidateAndAssignDefaults(default_settings)
+
+        self.model = model
+
+        evaluated_model_part_names = parameters["evaluated_model_part_names"].GetStringArray()
+        if len(evaluated_model_part_names) == 0:
+            raise RuntimeError(f"No model parts were provided for GeometricCentroidDeviationResponseFunction. [ response name = \"{self.GetName()}\"]")
+
+        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
+
+    def GetImplementedPhysicalKratosVariables(self) -> list[SupportedSensitivityFieldVariableTypes]:
+        return [KratosOA.SHAPE]
+
+    def GetEvaluatedModelPart(self) -> Kratos.ModelPart:
+        if self.model_part is None:
+            raise RuntimeError("Please call GeometricCentroidDeviationResponseFunction::Initialize first.")
+        return self.model_part
+
+    def GetAnalysisModelPart(self) -> None:
+        return None
+
+    def Initialize(self) -> None:
+        self.model_part = self.model_part_operation.GetModelPart()
+
+        self.model_part_center = Kratos.Array3(0.0)
+        number_of_nodes = self.model_part.GetCommunicator().GlobalNumberOfNodes()
+
+        if number_of_nodes == 0:
+            raise RuntimeError(f"The provided model part {self.model_part.FullName()} does not contain any nodes.")
+
+        for node in self.model_part.Nodes:
+            self.model_part_center[0] += node.X
+            self.model_part_center[1] += node.Y
+            self.model_part_center[2] += node.Z
+
+        self.model_part_center = self.model_part.GetCommunicator().GetDataCommunicator().SumAll(self.model_part_center)
+        self.model_part_center /= number_of_nodes
+
+    def Check(self) -> None:
+        pass
+
+    def Finalize(self) -> None:
+        pass
+
+    def CalculateValue(self) -> float:
+        average_location = Kratos.Array3(0.0)
+        for node in self.model_part.Nodes:
+            average_location[0] += node.X
+            average_location[1] += node.Y
+            average_location[2] += node.Z
+
+        average_location = self.model_part.GetCommunicator().GetDataCommunicator().SumAll(average_location)
+
+        number_of_nodes = self.model_part.GetCommunicator().GlobalNumberOfNodes()
+        self.value_array = (average_location / number_of_nodes  - self.model_part_center)
+        return self.value_array[0] ** 2 + self.value_array[1] ** 2 + self.value_array[2] ** 2
+
+    def CalculateGradient(self, physical_variable_collective_expressions: dict[SupportedSensitivityFieldVariableTypes, KratosOA.CollectiveExpression]) -> None:
+        # first merge all the model parts
+        merged_model_part_map = ModelPartUtilities.GetMergedMap(physical_variable_collective_expressions, False)
+
+        # now get the intersected model parts
+        intersected_model_part_map = ModelPartUtilities.GetIntersectedMap(self.model_part, merged_model_part_map, False)
+
+        number_of_nodes = self.model_part.GetCommunicator().GlobalNumberOfNodes()
+
+        # calculate the gradients
+        for physical_variable, merged_model_part in merged_model_part_map.items():
+            if physical_variable == KratosOA.SHAPE:
+                Kratos.VariableUtils().SetNonHistoricalVariableToZero(Kratos.SHAPE_SENSITIVITY, merged_model_part.Nodes)
+                Kratos.VariableUtils().SetNonHistoricalVariable(Kratos.SHAPE_SENSITIVITY, 2.0 * self.value_array / number_of_nodes, intersected_model_part_map[physical_variable].Nodes)
+                for container_expression in physical_variable_collective_expressions[physical_variable].GetContainerExpressions():
+                    if isinstance(container_expression, Kratos.Expression.NodalExpression):
+                        Kratos.Expression.VariableExpressionIO.Read(container_expression, Kratos.SHAPE_SENSITIVITY, False)
+                    else:
+                        raise RuntimeError(f"Requesting sensitivity w.r.t. SHAPE for a Container expression which is not a NodalExpression. [ Requested container expression = {container_expression} ].")
+            else:
+                raise RuntimeError(f"Unsupported sensitivity w.r.t. {physical_variable.Name()} requested. Followings are supported sensitivity variables:\n\tSHAPE")

applications/OptimizationApplication/tests/responses_tests/test_geometric_centroid_deviation_response_function.py

@@ -0,0 +1,61 @@
+import KratosMultiphysics as Kratos
+import KratosMultiphysics.OptimizationApplication as KratosOA
+
+import KratosMultiphysics.KratosUnittest as kratos_unittest
+from KratosMultiphysics.testing.utilities import ReadModelPart
+from KratosMultiphysics.OptimizationApplication.responses.geometric_centroid_deviation_response_function import GeometricCentroidDeviationResponseFunction
+
+class TestGeometricCentroidDeviationResponseFunction(kratos_unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        cls.model = Kratos.Model()
+        cls.model_part = cls.model.CreateModelPart("test")
+        with kratos_unittest.WorkFolderScope(".", __file__, True):
+            ReadModelPart("../model_part_utils_test/quads", cls.model_part)
+
+        cls.response_function = GeometricCentroidDeviationResponseFunction("geo_centroid", cls.model, Kratos.Parameters("""{"evaluated_model_part_names": ["test"]}"""))
+        cls.response_function.Initialize()
+        cls.response_function.Check()
+        cls.ref_value = cls.response_function.CalculateValue()
+
+    def _CheckSensitivity(self, update_method, expression_sensitivity_retrieval_method, delta, precision):
+        ref_value = self.response_function.CalculateValue()
+        for node in self.model_part.Nodes:
+            update_method(node, delta)
+            value = self.response_function.CalculateValue()
+            sensitivity = (value - ref_value)/delta
+            update_method(node, -delta)
+            self.assertAlmostEqual(sensitivity, expression_sensitivity_retrieval_method(node), precision)
+
+    def _UpdateNodalPositions(self, direction, entity, delta):
+        if direction == 0:
+            entity.X += delta
+        if direction == 1:
+            entity.Y += delta
+        if direction == 2:
+            entity.Z += delta
+
+    def test_CalculateValue(self):
+        self.assertAlmostEqual(self.ref_value, 0.0, 12)
+
+    def test_CalculateShapeSensitivity(self):
+        sensitivity = KratosOA.CollectiveExpression([Kratos.Expression.NodalExpression(self.model_part)])
+        self.response_function.CalculateGradient({KratosOA.SHAPE: sensitivity})
+        Kratos.Expression.VariableExpressionIO.Write(sensitivity.GetContainerExpressions()[0], KratosOA.SHAPE, False)
+
+        # calculate nodal shape sensitivities
+        self._CheckSensitivity(
+            lambda x, y: self._UpdateNodalPositions(0, x, y),
+            lambda x: x.GetValue(KratosOA.SHAPE_X),
+            1e-6,
+            4)
+
+        self._CheckSensitivity(
+            lambda x, y: self._UpdateNodalPositions(1, x, y),
+            lambda x: x.GetValue(KratosOA.SHAPE_Y),
+            1e-6,
+            4)
+
+if __name__ == "__main__":
+    Kratos.Tester.SetVerbosity(Kratos.Tester.Verbosity.PROGRESS)  # TESTS_OUTPUTS
+    kratos_unittest.main()

applications/OptimizationApplication/tests/test_OptimizationApplication.py

@@ -22,6 +22,7 @@
 import responses_tests.test_mass_response_function
 import responses_tests.test_linear_strain_energy_response_function
 import responses_tests.test_standardized_responses
+import responses_tests.test_geometric_centroid_deviation_response_function
 import test_model_part_utils
 import test_model_part_controllers
 import test_container_expression_utils
@@ -83,6 +84,7 @@ def AssembleTestSuites():
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([responses_tests.test_mass_response_function.TestMassResponseFunctionSolids]))
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([responses_tests.test_linear_strain_energy_response_function.TestLinearStrainEnergyResponseFunction]))
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([responses_tests.test_overhang_response_function.TestOverHangResponseFunction]))
+    smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([responses_tests.test_geometric_centroid_deviation_response_function.TestGeometricCentroidDeviationResponseFunction]))
 
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_container_expression.TestConditionPropertiesExpression]))
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_container_expression.TestElementPropertiesExpression]))