Merge pull request #11492 from KratosMultiphysics/Kratos_M30_Deliverable

Kratos m30 deliverable eFlows4HPC (New Node-Neighbour Element Retrieval Function)

Author: SADPR

applications/RomApplication/custom_python/add_custom_utilities_to_python.cpp

@@ -50,6 +50,7 @@ void AddCustomUtilitiesToPython(pybind11::module& m)
         .def_static("SetHRomComputingModelPart", &RomAuxiliaryUtilities::SetHRomComputingModelPart)
         .def_static("SetHRomVolumetricVisualizationModelPart", &RomAuxiliaryUtilities::SetHRomVolumetricVisualizationModelPart)
         .def_static("GetHRomConditionParentsIds", &RomAuxiliaryUtilities::GetHRomConditionParentsIds)
+        .def_static("GetNodalNeighbouringElementIdsNotInHRom", &RomAuxiliaryUtilities::GetNodalNeighbouringElementIdsNotInHRom)
         .def_static("GetConditionIdsNotInHRomModelPart", &RomAuxiliaryUtilities::GetConditionIdsNotInHRomModelPart)
         .def_static("GetElementIdsNotInHRomModelPart", &RomAuxiliaryUtilities::GetElementIdsNotInHRomModelPart)
         .def_static("GetHRomMinimumConditionsIds", &RomAuxiliaryUtilities::GetHRomMinimumConditionsIds)

applications/RomApplication/custom_strategies/global_rom_builder_and_solver.h

@@ -640,11 +640,10 @@ class GlobalROMBuilderAndSolver : public ResidualBasedBlockBuilderAndSolver<TSpa
 
         Build(pScheme, rModelPart, rA, rb);
 
-        if (!mHromSimulation)
+        if (mMonotonicityPreservingFlag) {
             BaseType::ApplyDirichletConditions(pScheme, rModelPart, rA, rDx, rb);
-
-        if (mMonotonicityPreservingFlag)
             MonotonicityPreserving(rA, rb);
+        }
 
         ProjectROM(rModelPart, rA, rb);
 

applications/RomApplication/custom_utilities/rom_auxiliary_utilities.cpp

@@ -335,6 +335,34 @@ std::vector<IndexType> RomAuxiliaryUtilities::GetHRomConditionParentsIds(
     return parent_ids;
 }
 
+std::vector<IndexType> RomAuxiliaryUtilities::GetNodalNeighbouringElementIdsNotInHRom(
+    ModelPart& rModelPart,
+    ModelPart& rGivenModelPart,
+    const std::map<std::string, std::map<IndexType, double>>& rHRomWeights)
+{
+    std::vector<IndexType> new_element_ids;
+    const auto& r_elem_weights = rHRomWeights.at("Elements");
+    
+    FindGlobalNodalEntityNeighboursProcess<ModelPart::ElementsContainerType> find_nodal_elements_neighbours_process(rModelPart);
+    find_nodal_elements_neighbours_process.Execute();
+
+    for (const auto& r_node : rGivenModelPart.Nodes()) {
+        const auto& r_neigh = r_node.GetValue(NEIGHBOUR_ELEMENTS);
+
+        // Add the neighbour elements to the HROM weights
+        for (size_t i = 0; i < r_neigh.size(); ++i) {
+            const auto& r_elem = r_neigh[i];
+
+            // Note that we check if the element has been already added by the HROM element selection strategy
+            if (r_elem_weights.find(r_elem.Id() - 1) == r_elem_weights.end()) { //FIXME: FIX THE + 1 --> WE NEED TO WRITE REAL IDS IN THE WEIGHTS!!
+                new_element_ids.push_back(r_elem.Id() - 1); //FIXME: FIX THE + 1 --> WE NEED TO WRITE REAL IDS IN THE WEIGHTS!!
+            }
+        }
+    }
+
+    return new_element_ids;
+}
+
 std::vector<IndexType> RomAuxiliaryUtilities::GetElementIdsNotInHRomModelPart(
     const ModelPart& rModelPart,
     const ModelPart& rModelPartWithElementsToInclude,

applications/RomApplication/custom_utilities/rom_auxiliary_utilities.h

@@ -26,6 +26,7 @@
 #include "includes/model_part.h"
 #include "includes/ublas_interface.h"
 #include "modified_shape_functions/modified_shape_functions.h"
+#include "processes/find_nodal_neighbours_process.h"
 
 // Application includes
 #include "rom_application_variables.h"
@@ -106,6 +107,25 @@ class KRATOS_API(ROM_APPLICATION) RomAuxiliaryUtilities
     static std::vector<IndexType> GetHRomConditionParentsIds(
         const ModelPart& rModelPart,
         const std::map<std::string, std::map<IndexType, double>>& rHRomWeights);
+    
+    /**
+     * @brief Retrieve the IDs of elements neighboring nodes in a given sub-model part but not present in HRom weights.
+     *
+     * This function iterates over all the nodes in the provided sub-model part (`rGivenModelPart`) and collects 
+     * the IDs of the elements neighboring each node. The neighboring elements are determined using the 
+     * 'NEIGHBOUR_ELEMENTS' value attached to each node. The function then checks if these elements are already 
+     * present in `rHRomWeights`. If not, their IDs are added to the return vector. It's important to note that 
+     * this function assumes that the 'NEIGHBOUR_ELEMENTS' values are already computed for the nodes in the model part.
+     *
+     * @param rModelPart The complete model part which houses all the elements.
+     * @param rGivenModelPart The sub-model part with nodes for which neighboring elements should be fetched.
+     * @param rHRomWeights Map containing the original HROM conditions and elements weights.
+     * @return std::vector<IndexType> A list of unique neighboring element IDs not already present in HRom weights.
+     */
+    static std::vector<IndexType> GetNodalNeighbouringElementIdsNotInHRom(
+        ModelPart& rModelPart,
+        ModelPart& rGivenModelPart,
+        const std::map<std::string, std::map<IndexType, double>>& rHRomWeights);
 
     /**
      * @brief Identifies condition IDs from a given ModelPart that are not in the HROM weights

applications/RomApplication/python_scripts/hrom_training_utility.py

@@ -49,6 +49,7 @@ def __init__(self, solver, custom_settings):
         # Retrieve list of model parts from settings
         self.include_conditions_model_parts_list = settings["include_conditions_model_parts_list"].GetStringArray()
         self.include_elements_model_parts_list = settings["include_elements_model_parts_list"].GetStringArray()
+        self.include_nodal_neighbouring_elements_model_parts_list = settings["include_nodal_neighbouring_elements_model_parts_list"].GetStringArray()
 
         # Check if the model parts exist
         for model_part_name in self.include_conditions_model_parts_list:
@@ -132,7 +133,7 @@ def CreateHRomModelParts(self):
 
         # Output the HROM model part in mdpa format
         hrom_output_name = "{}HROM".format(model_part_output_name)
-        model_part_io = KratosMultiphysics.ModelPartIO(hrom_output_name, KratosMultiphysics.IO.WRITE | KratosMultiphysics.IO.MESH_ONLY)
+        model_part_io = KratosMultiphysics.ModelPartIO(hrom_output_name, KratosMultiphysics.IO.WRITE | KratosMultiphysics.IO.MESH_ONLY | KratosMultiphysics.IO.SCIENTIFIC_PRECISION)
         model_part_io.WriteModelPart(hrom_main_model_part)
         KratosMultiphysics.kratos_utilities.DeleteFileIfExisting("{}.time".format(hrom_output_name))
         if self.echo_level > 0:
@@ -153,7 +154,7 @@ def CreateHRomModelParts(self):
 
             # Write the HROM visualization mesh
             hrom_vis_output_name = "{}HROMVisualization".format(model_part_output_name)
-            model_part_io = KratosMultiphysics.ModelPartIO(hrom_vis_output_name, KratosMultiphysics.IO.WRITE | KratosMultiphysics.IO.MESH_ONLY)
+            model_part_io = KratosMultiphysics.ModelPartIO(hrom_vis_output_name, KratosMultiphysics.IO.WRITE | KratosMultiphysics.IO.MESH_ONLY | KratosMultiphysics.IO.SCIENTIFIC_PRECISION)
             model_part_io.WriteModelPart(hrom_visualization_model_part)
             KratosMultiphysics.kratos_utilities.DeleteFileIfExisting("{}.time".format(hrom_vis_output_name))
             if self.echo_level > 0:
@@ -170,6 +171,7 @@ def __GetHRomTrainingDefaultSettings(cls):
             "projection_strategy": "galerkin",
             "include_conditions_model_parts_list": [],
             "include_elements_model_parts_list": [],
+            "include_nodal_neighbouring_elements_model_parts_list":[],
             "include_minimum_condition": false,
             "include_condition_parents": true
         }""")
@@ -247,6 +249,25 @@ def AppendHRomWeightsToRomParameters(self):
 
                 # If needed, update your weights and indexes using __AddSelectedElementsWithZeroWeights function with the new_elements
                 weights, indexes = self.__AddSelectedElementsWithZeroWeights(weights, indexes, new_elements)
+        
+        # Add nodal neighbouring elements
+        for model_part_name in self.include_nodal_neighbouring_elements_model_parts_list:
+            # Check if the sub model part exists
+            if self.solver.model.HasModelPart(model_part_name):
+                nodal_neighbours_model_part = self.solver.model.GetModelPart(model_part_name)
+
+                # Call the GetNodalNeighbouringElementIdsNotInHRom function
+                new_nodal_neighbours = KratosROM.RomAuxiliaryUtilities.GetNodalNeighbouringElementIdsNotInHRom(
+                    root_model_part, # The complete model part
+                    nodal_neighbours_model_part, # The model part containing the nodal neighbouring elements to be included
+                    hrom_weights)
+
+                # Add the new nodal neighbouring elements to the elements dict with a null weight
+                for element_id in new_nodal_neighbours:
+                    hrom_weights["Elements"][element_id] = 0.0
+
+                # If needed, update your weights and indexes using __AddSelectedElementsWithZeroWeights function with the new_nodal_neighbours
+                weights, indexes = self.__AddSelectedElementsWithZeroWeights(weights, indexes, new_nodal_neighbours)
 
         # If required, keep at least one condition per submodelpart
         # This might be required by those BCs involving the faces (e.g. slip BCs)