diff --git a/applications/OptimizationApplication/python_scripts/responses/geometric_centroid_deviation_response_function.py b/applications/OptimizationApplication/python_scripts/responses/geometric_centroid_deviation_response_function.py new file mode 100644 index 000000000000..a174ba94bab9 --- /dev/null +++ b/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") diff --git a/applications/OptimizationApplication/tests/responses_tests/test_geometric_centroid_deviation_response_function.py b/applications/OptimizationApplication/tests/responses_tests/test_geometric_centroid_deviation_response_function.py new file mode 100644 index 000000000000..68ec4d2392c9 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/applications/OptimizationApplication/tests/test_OptimizationApplication.py b/applications/OptimizationApplication/tests/test_OptimizationApplication.py index bbd74ad9d259..b0f94feaab19 100644 --- a/applications/OptimizationApplication/tests/test_OptimizationApplication.py +++ b/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]))