[FastPR][ConvDiff] Enhance element substitution

applications/ConvectionDiffusionApplication/convection_diffusion_application.cpp

@@ -34,9 +34,9 @@ KratosConvectionDiffusionApplication::KratosConvectionDiffusionApplication()
     : KratosApplication("ConvectionDiffusionApplication"),
       mAxisymmetricEulerianConvectionDiffusion2D3N(0, Element::GeometryType::Pointer(new Triangle2D3<Node >(Element::GeometryType::PointsArrayType(3)))),
       mAxisymmetricEulerianConvectionDiffusion2D4N(0, Element::GeometryType::Pointer(new Quadrilateral2D4<Node >(Element::GeometryType::PointsArrayType(4)))),
-      mEulerianConvDiff2D(0, Element::GeometryType::Pointer(new Triangle2D3<Node >(Element::GeometryType::PointsArrayType(3)))),
+      mEulerianConvDiff2D3N(0, Element::GeometryType::Pointer(new Triangle2D3<Node >(Element::GeometryType::PointsArrayType(3)))),
       mEulerianConvDiff2D4N(0, Element::GeometryType::Pointer(new Quadrilateral2D4<Node >(Element::GeometryType::PointsArrayType(4)))),
-      mEulerianConvDiff3D(0, Element::GeometryType::Pointer(new Tetrahedra3D4<Node >(Element::GeometryType::PointsArrayType(4)))),
+      mEulerianConvDiff3D4N(0, Element::GeometryType::Pointer(new Tetrahedra3D4<Node >(Element::GeometryType::PointsArrayType(4)))),
       mEulerianConvDiff3D8N(0, Element::GeometryType::Pointer(new Hexahedra3D8<Node >(Element::GeometryType::PointsArrayType(8)))),
       mEulerianDiffusion2D3N(0, Element::GeometryType::Pointer(new Triangle2D3<Node >(Element::GeometryType::PointsArrayType(3)))),
       mEulerianDiffusion3D4N(0, Element::GeometryType::Pointer(new Tetrahedra3D4<Node >(Element::GeometryType::PointsArrayType(4)))),
@@ -95,9 +95,11 @@ void KratosConvectionDiffusionApplication::Register() {
     // Registering elements and conditions here
     KRATOS_REGISTER_ELEMENT("AxisymmetricEulerianConvectionDiffusion2D3N", mAxisymmetricEulerianConvectionDiffusion2D3N);
     KRATOS_REGISTER_ELEMENT("AxisymmetricEulerianConvectionDiffusion2D4N", mAxisymmetricEulerianConvectionDiffusion2D4N);
-    KRATOS_REGISTER_ELEMENT("EulerianConvDiff2D", mEulerianConvDiff2D);
+    KRATOS_REGISTER_ELEMENT("EulerianConvDiff2D", mEulerianConvDiff2D3N); //TODO: To be removed as it does not follow the naming convention
+    KRATOS_REGISTER_ELEMENT("EulerianConvDiff2D3N", mEulerianConvDiff2D3N);
     KRATOS_REGISTER_ELEMENT("EulerianConvDiff2D4N", mEulerianConvDiff2D4N);
-    KRATOS_REGISTER_ELEMENT("EulerianConvDiff3D", mEulerianConvDiff3D);
+    KRATOS_REGISTER_ELEMENT("EulerianConvDiff3D", mEulerianConvDiff3D4N); //TODO: To be removed as it does not follow the naming convention
+    KRATOS_REGISTER_ELEMENT("EulerianConvDiff3D4N", mEulerianConvDiff3D4N);
     KRATOS_REGISTER_ELEMENT("EulerianConvDiff3D8N", mEulerianConvDiff3D8N);
     KRATOS_REGISTER_ELEMENT("EulerianDiffusion2D3N", mEulerianDiffusion2D3N);
     KRATOS_REGISTER_ELEMENT("EulerianDiffusion3D4N", mEulerianDiffusion3D4N);

applications/ConvectionDiffusionApplication/convection_diffusion_application.h

@@ -233,9 +233,9 @@ class KRATOS_API(CONVECTION_DIFFUSION_APPLICATION) KratosConvectionDiffusionAppl
     const AxisymmetricEulerianConvectionDiffusionElement<2,3>  mAxisymmetricEulerianConvectionDiffusion2D3N;
     const AxisymmetricEulerianConvectionDiffusionElement<2,4>  mAxisymmetricEulerianConvectionDiffusion2D4N;
 
-    const EulerianConvectionDiffusionElement<2,3>  mEulerianConvDiff2D;
+    const EulerianConvectionDiffusionElement<2,3>  mEulerianConvDiff2D3N;
     const EulerianConvectionDiffusionElement<2,4>  mEulerianConvDiff2D4N;
-    const EulerianConvectionDiffusionElement<3,4>  mEulerianConvDiff3D;
+    const EulerianConvectionDiffusionElement<3,4>  mEulerianConvDiff3D4N;
     const EulerianConvectionDiffusionElement<3,8>  mEulerianConvDiff3D8N;
     const EulerianDiffusionElement<2,3>  mEulerianDiffusion2D3N;
     const EulerianDiffusionElement<3,4>  mEulerianDiffusion3D4N;

applications/ConvectionDiffusionApplication/python_scripts/convection_diffusion_solver.py

@@ -533,55 +533,49 @@ def _set_and_fill_buffer(self):
             self.main_model_part.CloneTimeStep(time)
 
     def _get_element_condition_replace_settings(self):
+        ## Get and check domain size
         domain_size = self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]
+        if domain_size not in [2,3]:
+            raise Exception("DOMAIN_SIZE is not set in ProcessInfo container.")
 
-        ## Elements
-        num_nodes_elements = 0
-        if (len(self.main_model_part.Elements) > 0):
-            for elem in self.main_model_part.Elements:
-                num_nodes_elements = len(elem.GetNodes())
-                break
-
-        num_nodes_elements = self.main_model_part.GetCommunicator().GetDataCommunicator().MaxAll(num_nodes_elements)
+        ## Validate the replace settings
+        default_replace_settings = self.GetDefaultParameters()["element_replace_settings"]
+        self.settings["element_replace_settings"].ValidateAndAssignDefaults(default_replace_settings)
 
+        ## Elements
+        ## Note that we check for the elements that require substitution to allow for custom elements
         element_name = self.settings["element_replace_settings"]["element_name"].GetString()
+        element_list = ["EulerianConvDiff","LaplacianElement","MixedLaplacianElement","AdjointHeatDiffusionElement","QSConvectionDiffusionExplicit","DConvectionDiffusionExplicit","AxisymmetricEulerianConvectionDiffusion"]
+        if element_name in element_list:
+            num_nodes_elements = 0
+            if (len(self.main_model_part.Elements) > 0):
+                for elem in self.main_model_part.Elements:
+                    num_nodes_elements = len(elem.GetNodes())
+                    break
 
-        if domain_size not in (2,3):
-            raise Exception("DOMAIN_SIZE not set")
+            num_nodes_elements = self.main_model_part.GetCommunicator().GetDataCommunicator().MaxAll(num_nodes_elements)
+            if not num_nodes_elements:
+                num_nodes_elements = domain_size + 1
 
-        if element_name == "EulerianConvDiff":
-            if domain_size == 2:
-                if num_nodes_elements == 3:
-                    self.settings["element_replace_settings"]["element_name"].SetString("EulerianConvDiff2D")
-                else:
-                    self.settings["element_replace_settings"]["element_name"].SetString("EulerianConvDiff2D4N")
-            else:
-                if num_nodes_elements == 4:
-                    self.settings["element_replace_settings"]["element_name"].SetString("EulerianConvDiff3D")
-                else:
-                    self.settings["element_replace_settings"]["element_name"].SetString("EulerianConvDiff3D8N")
-        elif element_name in ("LaplacianElement","MixedLaplacianElement","AdjointHeatDiffusionElement","QSConvectionDiffusionExplicit","DConvectionDiffusionExplicit"):
-            name_string = "{0}{1}D{2}N".format(element_name,domain_size, num_nodes_elements)
+            name_string = f"{element_name}{domain_size}D{num_nodes_elements}N"
             self.settings["element_replace_settings"]["element_name"].SetString(name_string)
 
         ## Conditions
-        num_conditions = self.main_model_part.GetCommunicator().GetDataCommunicator().SumAll(len(self.main_model_part.Conditions))
-
-        if num_conditions > 0:
+        condition_name = self.settings["element_replace_settings"]["condition_name"].GetString()
+        condition_list = ["FluxCondition","ThermalFace","AxisymmetricThermalFace","LineCondition","SurfaceCondition"]
+        if condition_name in condition_list:
             num_nodes_conditions = 0
             if (len(self.main_model_part.Conditions) > 0):
                 for cond in self.main_model_part.Conditions:
                     num_nodes_conditions = len(cond.GetNodes())
                     break
 
             num_nodes_conditions = self.main_model_part.GetCommunicator().GetDataCommunicator().MaxAll(num_nodes_conditions)
+            if not num_nodes_conditions:
+                num_nodes_conditions = domain_size
 
-            condition_name = self.settings["element_replace_settings"]["condition_name"].GetString()
-            if condition_name in ("FluxCondition","ThermalFace","Condition"):
-                name_string = "{0}{1}D{2}N".format(condition_name,domain_size, num_nodes_conditions)
-                self.settings["element_replace_settings"]["condition_name"].SetString(name_string)
-        else:
-            self.settings["element_replace_settings"]["condition_name"].SetString("")
+            name_string = f"{condition_name}{domain_size}D{num_nodes_conditions}N"
+            self.settings["element_replace_settings"]["condition_name"].SetString(name_string)
 
         return self.settings["element_replace_settings"]
 

applications/ConvectionDiffusionApplication/python_scripts/coupled_fluid_thermal_solver.py

@@ -1,11 +1,9 @@
-import sys
-
 # Importing the Kratos Library
 import KratosMultiphysics
 
-# Import applications
-import KratosMultiphysics.FluidDynamicsApplication as KratosCFD
-import KratosMultiphysics.ConvectionDiffusionApplication as ConvDiff
+# Import applications modules
+from KratosMultiphysics.FluidDynamicsApplication import python_solvers_wrapper_fluid
+from KratosMultiphysics.ConvectionDiffusionApplication import python_solvers_wrapper_convection_diffusion
 
 # Importing the base class
 from KratosMultiphysics.python_solver import PythonSolver
@@ -31,13 +29,13 @@ def GetDefaultParameters(cls):
                 }
             },
             "thermal_solver_settings": {
-                "solver_type": "Transient",
+                "solver_type": "transient",
                 "analysis_type": "linear",
                 "model_import_settings": {
                     "input_type": "use_input_model_part"
                 },
                 "material_import_settings": {
-                        "materials_filename": "ThermalMaterials.json"
+                    "materials_filename": "ThermalMaterials.json"
                 }
             }
         }
@@ -47,16 +45,14 @@ def GetDefaultParameters(cls):
         return default_settings
 
     def __init__(self, model, custom_settings):
-
-        super(CoupledFluidThermalSolver, self).__init__(model, custom_settings)
+        ## Cal base class constructor
+        super().__init__(model, custom_settings)
 
         ## Get domain size
         self.domain_size = self.settings["domain_size"].GetInt()
 
-        from KratosMultiphysics.FluidDynamicsApplication import python_solvers_wrapper_fluid
+        ## Create subdomain solvers
         self.fluid_solver = python_solvers_wrapper_fluid.CreateSolverByParameters(self.model, self.settings["fluid_solver_settings"],"OpenMP")
-
-        from KratosMultiphysics.ConvectionDiffusionApplication import python_solvers_wrapper_convection_diffusion
         self.thermal_solver = python_solvers_wrapper_convection_diffusion.CreateSolverByParameters(self.model,self.settings["thermal_solver_settings"],"OpenMP")
 
     def AddVariables(self):
@@ -73,17 +69,13 @@ def ImportModelPart(self):
         convection_diffusion_settings = self.thermal_solver.main_model_part.ProcessInfo.GetValue(KratosMultiphysics.CONVECTION_DIFFUSION_SETTINGS)
 
         # Here the fluid model part is cloned to be thermal model part so that the nodes are shared
+        element_name, condition_name = self.__GetElementAndConditionNames()
         modeler = KratosMultiphysics.ConnectivityPreserveModeler()
-        if self.domain_size == 2:
-            modeler.GenerateModelPart(self.fluid_solver.main_model_part,
-                                      self.thermal_solver.main_model_part,
-                                      "EulerianConvDiff2D",
-                                      "ThermalFace2D2N")
-        else:
-            modeler.GenerateModelPart(self.fluid_solver.main_model_part,
-                                      self.thermal_solver.main_model_part,
-                                      "EulerianConvDiff3D",
-                                      "ThermalFace3D3N")
+        modeler.GenerateModelPart(
+            self.fluid_solver.main_model_part,
+            self.thermal_solver.main_model_part,
+            element_name,
+            condition_name)
 
         # Set the saved convection diffusion settings to the new thermal model part
         self.thermal_solver.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.CONVECTION_DIFFUSION_SETTINGS, convection_diffusion_settings)
@@ -152,8 +144,46 @@ def FinalizeSolutionStep(self):
         self.fluid_solver.FinalizeSolutionStep()
         self.thermal_solver.FinalizeSolutionStep()
 
-    def Solve(self):
-        self.InitializeSolutionStep()
-        self.Predict()
-        self.SolveSolutionStep()
-        self.FinalizeSolutionStep()
+    def __GetElementAndConditionNames(self):
+        ''' Auxiliary function to get the element and condition names for the connectivity preserve modeler call
+
+        This function returns the element and condition names from the domain size and number of nodes.
+        Note that throughout all the substitution process a unique element type and condition is assumed.
+        Also note that the connectivity preserve modeler call will create standard base elements as these are to
+        be substituted by the corresponding ones in the PrepareModelPart call of the thermal solver.
+        '''
+
+        ## Get and check domain size
+        fluid_model_part = self.fluid_solver.main_model_part
+        domain_size = fluid_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]
+        if domain_size not in [2,3]:
+            raise Exception("DOMAIN_SIZE is not set in ProcessInfo container.")
+
+        ## Elements
+        ## Get the number of nodes from the fluid mesh elements (if there are no elements simplicial are assumed)
+        num_nodes_elements = 0
+        if (len(fluid_model_part.Elements) > 0):
+            for elem in fluid_model_part.Elements:
+                num_nodes_elements = len(elem.GetNodes())
+                break
+        num_nodes_elements = fluid_model_part.GetCommunicator().GetDataCommunicator().MaxAll(num_nodes_elements)
+        if not num_nodes_elements:
+            num_nodes_elements = domain_size + 1
+
+        element_name = f"Element{domain_size}D{num_nodes_elements}N"
+
+        ## Conditions
+        ## Get the number of nodes from the fluid mesh conditions (if there are no elements simplicial are assumed)
+        num_nodes_conditions = 0
+        if (len(fluid_model_part.Conditions) > 0):
+            for cond in fluid_model_part.Conditions:
+                num_nodes_conditions = len(cond.GetNodes())
+                break
+        num_nodes_conditions = fluid_model_part.GetCommunicator().GetDataCommunicator().MaxAll(num_nodes_conditions)
+        if not num_nodes_conditions:
+            num_nodes_conditions = domain_size
+
+        aux_condition_name = "LineCondition" if domain_size == 2 else "SurfaceCondition"
+        condition_name = f"{aux_condition_name}{domain_size}D{num_nodes_conditions}N"
+
+        return element_name, condition_name