[RomApp] Adding candidate elements (and/or conditions) to HROM

applications/RomApplication/custom_python/add_custom_utilities_to_python.cpp

@@ -54,6 +54,8 @@ void AddCustomUtilitiesToPython(pybind11::module& m)
         .def_static("GetElementIdsNotInHRomModelPart", &RomAuxiliaryUtilities::GetElementIdsNotInHRomModelPart)
         .def_static("GetHRomMinimumConditionsIds", &RomAuxiliaryUtilities::GetHRomMinimumConditionsIds)
         .def_static("ProjectRomSolutionIncrementToNodes", &RomAuxiliaryUtilities::ProjectRomSolutionIncrementToNodes)
+        .def_static("GetElementIdsInModelPart", &RomAuxiliaryUtilities::GetElementIdsInModelPart)
+        .def_static("GetConditionIdsInModelPart", &RomAuxiliaryUtilities::GetConditionIdsInModelPart)
         ;
 }
 

applications/RomApplication/custom_utilities/rom_auxiliary_utilities.cpp

@@ -342,7 +342,7 @@ std::vector<IndexType> RomAuxiliaryUtilities::GetNodalNeighbouringElementIdsNotI
 {
     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();
 
@@ -373,7 +373,7 @@ std::vector<IndexType> RomAuxiliaryUtilities::GetElementIdsNotInHRomModelPart(
 
     for (const auto& r_elem : rModelPartWithElementsToInclude.Elements()) {
         IndexType element_id = r_elem.Id();
-        
+
         // Check if the element is already added
         if (r_elem_weights.find(element_id - 1) == r_elem_weights.end()) {
             new_element_ids.push_back(element_id - 1);
@@ -383,6 +383,7 @@ std::vector<IndexType> RomAuxiliaryUtilities::GetElementIdsNotInHRomModelPart(
     return new_element_ids;
 }
 
+
 std::vector<IndexType> RomAuxiliaryUtilities::GetConditionIdsNotInHRomModelPart(
     const ModelPart& rModelPart,
     const ModelPart& rModelPartWithConditionsToInclude,
@@ -393,7 +394,7 @@ std::vector<IndexType> RomAuxiliaryUtilities::GetConditionIdsNotInHRomModelPart(
 
     for (const auto& r_cond : rModelPartWithConditionsToInclude.Conditions()) {
         IndexType condition_id = r_cond.Id();
-        
+
         // Check if the condition is already added
         if (r_cond_weights.find(condition_id - 1) == r_cond_weights.end()) {
             new_condition_ids.push_back(condition_id - 1);
@@ -403,6 +404,28 @@ std::vector<IndexType> RomAuxiliaryUtilities::GetConditionIdsNotInHRomModelPart(
     return new_condition_ids;
 }
 
+std::vector<IndexType> RomAuxiliaryUtilities::GetElementIdsInModelPart(
+    const ModelPart& rModelPart)
+{
+    std::vector<IndexType> element_ids;
+
+    for (const auto& r_elem : rModelPart.Elements()) {
+        element_ids.push_back(r_elem.Id());
+    }
+    return element_ids;
+}
+
+std::vector<IndexType> RomAuxiliaryUtilities::GetConditionIdsInModelPart(
+    const ModelPart& rModelPart)
+{
+    std::vector<IndexType> condition_ids;
+
+    for (const auto& r_cond : rModelPart.Conditions()) {
+        condition_ids.push_back(r_cond.Id());
+    }
+    return condition_ids;
+}
+
 std::vector<IndexType> RomAuxiliaryUtilities::GetHRomMinimumConditionsIds(
     const ModelPart& rModelPart,
     const std::map<IndexType, double>& rHRomConditionWeights)
@@ -567,6 +590,6 @@ void RomAuxiliaryUtilities::GetPsiElemental(
                 noalias(row(rPsiElemental, i)) = row(r_nodal_rom_basis, row_id);
             }
         }
-    } 
+    }
 
 } // namespace Kratos

applications/RomApplication/custom_utilities/rom_auxiliary_utilities.h

@@ -107,14 +107,14 @@ 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 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.
@@ -168,6 +168,15 @@ class KRATOS_API(ROM_APPLICATION) RomAuxiliaryUtilities
         const ModelPart& rModelPart,
         const std::map<IndexType, double>& rHRomConditionWeights);
 
+
+    static std::vector<IndexType> GetElementIdsInModelPart(
+        const ModelPart& rModelPart
+    );
+
+    static std::vector<IndexType> GetConditionIdsInModelPart(
+        const ModelPart& rModelPart
+    );
+
     /**
      * @brief Project the ROM solution increment onto the nodal basis
      * For a given model part this function takes the ROM_SOLUTION_INCREMENT, which is assumed to be

applications/RomApplication/python_scripts/empirical_cubature_method.py

@@ -14,31 +14,42 @@ class EmpiricalCubatureMethod():
     Reference: Hernandez 2020. "A multiscale method for periodic structures using domain decomposition and ECM-hyperreduction"
     """
 
-
-    """
-    Constructor setting up the parameters for the Element Selection Strategy
-        ECM_tolerance: approximation tolerance for the element selection algorithm
-        Filter_tolerance: parameter limiting the number of candidate points (elements) to those above this tolerance
-        Plotting: whether to plot the error evolution of the element selection algorithm
-    """
-    def __init__(self, ECM_tolerance = 0, Filter_tolerance = 1e-16, Plotting = False):
+    def __init__(
+        self,
+        ECM_tolerance = 0,
+        Filter_tolerance = 1e-16,
+        Plotting = False,
+        MaximumNumberUnsuccesfulIterations = 100
+    ):
+        """
+        Constructor setting up the parameters for the Element Selection Strategy
+            ECM_tolerance: approximation tolerance for the element selection algorithm
+            Filter_tolerance: parameter limiting the number of candidate points (elements) to those above this tolerance
+            Plotting: whether to plot the error evolution of the element selection algorithm
+        """
         self.ECM_tolerance = ECM_tolerance
         self.Filter_tolerance = Filter_tolerance
         self.Name = "EmpiricalCubature"
         self.Plotting = Plotting
-
-
-
-    """
-    Method for setting up the element selection
-    input:  - ResidualsBasis: numpy array containing a basis to the residuals projected
-            - constrain_sum_of_weights: enable the user to constrain weights to be the sum of the number of entities.
-            - constrain_conditions: enable the user to enforce weights to consider conditions (for specific boundary conditions).
-    """
-    def SetUp(self, ResidualsBasis, constrain_sum_of_weights=True, constrain_conditions = False, number_of_conditions = 0):
-
+        self.MaximumNumberUnsuccesfulIterations = MaximumNumberUnsuccesfulIterations
+
+    def SetUp(
+        self,
+        ResidualsBasis,
+        InitialCandidatesSet = None,
+        constrain_sum_of_weights=True,
+        constrain_conditions = False,
+        number_of_conditions = 0
+    ):
+        """
+        Method for setting up the element selection
+        input:  - ResidualsBasis: numpy array containing a basis to the residuals projected
+                - constrain_sum_of_weights: enable the user to constrain weights to be the sum of the number of entities.
+                - constrain_conditions: enable the user to enforce weights to consider conditions (for specific boundary conditions).
+        """
         self.W = np.ones(np.shape(ResidualsBasis)[0])
         self.G = ResidualsBasis.T
+        self.y = InitialCandidatesSet
         self.add_constrain_count = None
         total_number_of_entities = np.shape(self.G)[1]
         elements_constraint = np.ones(total_number_of_entities)
@@ -66,47 +77,72 @@ def SetUp(self, ResidualsBasis, constrain_sum_of_weights=True, constrain_conditi
             self.G = np.vstack([ self.G , projection_of_constant_vector_conditions ] )
             self.add_constrain_count = -2
         self.b = self.G @ self.W
+        self.UnsuccesfulIterations = 0
 
-
-    """
-    Method performing calculations required before launching the Calculate method
-    """
     def Initialize(self):
-        self.Gnorm = np.sqrt(sum(np.multiply(self.G, self.G), 0))
+        """
+        Method performing calculations required before launching the Calculate method
+        """
+        self.GnormNOONE = np.linalg.norm(self.G[:self.add_constrain_count,:], axis = 0)
         M = np.shape(self.G)[1]
         normB = np.linalg.norm(self.b)
-        self.y = np.arange(0,M,1) # Set of candidate points (those whose associated column has low norm are removed)
-        GnormNOONE = np.sqrt(sum(np.multiply(self.G[:self.add_constrain_count,:], self.G[:self.add_constrain_count,:]), 0))
-        if self.Filter_tolerance > 0:
-            TOL_REMOVE = self.Filter_tolerance * normB
-            rmvpin = np.where(GnormNOONE[self.y] < TOL_REMOVE)
-            self.y = np.delete(self.y,rmvpin)
+
+        if self.y is None:
+            self.y = np.arange(0,M,1) # Set of candidate points (those whose associated column has low norm are removed)
+
+            if self.Filter_tolerance > 0:
+                TOL_REMOVE = self.Filter_tolerance * normB
+                rmvpin = np.where(self.GnormNOONE[self.y] < TOL_REMOVE)
+                #self.y_complement = self.y[rmvpin]
+                self.y = np.delete(self.y,rmvpin)
+        else:
+            self.y_complement = np.arange(0,M,1)
+            self.y_complement = np.delete(self.y_complement, self.y)# Set of candidate points (those whose associated column has low norm are removed)
+            if self.Filter_tolerance > 0:
+                TOL_REMOVE = self.Filter_tolerance * normB
+                rmvpin = np.where(self.GnormNOONE[self.y_complement] < TOL_REMOVE)
+                self.y_complement = np.delete(self.y_complement,rmvpin)
+
         self.z = {}  # Set of intergration points
         self.mPOS = 0 # Number of nonzero weights
-        self.r = self.b # residual vector
+        self.r = self.b.copy() # residual vector
         self.m = len(self.b) # Default number of points
         self.nerror = np.linalg.norm(self.r)/normB
         self.nerrorACTUAL = self.nerror
 
-
     def Run(self):
         self.Initialize()
         self.Calculate()
 
+    def expand_candidates_with_complement(self):
+        self.y = np.r_[self.y,self.y_complement]
+        print('expanding set to include the complement...')
+        ExpandedSetFlag = True
+        return ExpandedSetFlag
 
-    """
-    Method launching the element selection algorithm to find a set of elements: self.z, and wiegths: self.w
-    """
     def Calculate(self):
-
+        """
+        Method launching the element selection algorithm to find a set of elements: self.z, and wiegths: self.w
+        """
+        MaximumLengthZ = 0
+        ExpandedSetFlag = False
         k = 1 # number of iterations
-        while self.nerrorACTUAL > self.ECM_tolerance and self.mPOS < self.m and len(self.y) != 0:
+        self.success = True
+        while self.nerrorACTUAL > self.ECM_tolerance and self.mPOS < self.m and np.size(self.y) != 0:
+
+            if  self.UnsuccesfulIterations >  self.MaximumNumberUnsuccesfulIterations and not ExpandedSetFlag:
+                ExpandedSetFlag = self.expand_candidates_with_complement()
 
             #Step 1. Compute new point
-            ObjFun = self.G[:,self.y].T @ self.r.T
-            ObjFun = ObjFun.T / self.Gnorm[self.y]
-            indSORT = np.argmax(ObjFun)
-            i = self.y[indSORT]
+            if np.size(self.y)==1:
+                #candidate set consists of a single element
+                indSORT = 0
+                i = int(self.y)
+            else:
+                ObjFun = self.G[:,self.y].T @ self.r.T
+                ObjFun = ObjFun.T / self.GnormNOONE[self.y]
+                indSORT = np.argmax(ObjFun)
+                i = self.y[indSORT]
             if k==1:
                 alpha = np.linalg.lstsq(self.G[:, [i]], self.b)[0]
                 H = 1/(self.G[:,i] @ self.G[:,i].T)
@@ -118,7 +154,13 @@ def Calculate(self):
                 self.z = i
             else:
                 self.z = np.r_[self.z,i]
-            self.y = np.delete(self.y,indSORT)
+
+            #self.y = np.delete(self.y,indSORT)
+            if np.size(self.y)==1:
+                self.expand_candidates_with_complement()
+                self.y = np.delete(self.y,indSORT)
+            else:
+                self.y = np.delete(self.y,indSORT)
 
             # Step 4. Find possible negative weights
             if any(alpha < 0):
@@ -130,9 +172,15 @@ def Calculate(self):
                 alpha = H @ (self.G[:, self.z].T @ self.b)
                 alpha = alpha.reshape(len(alpha),1)
 
+            if np.size(self.z) > MaximumLengthZ :
+                self.UnsuccesfulIterations = 0
+            else:
+                self.UnsuccesfulIterations += 1
+
             #Step 6 Update the residual
-            if len(alpha)==1:
-                self.r = self.b - (self.G[:,self.z] * alpha)
+            if np.size(alpha)==1:
+                self.r = self.b.reshape(-1,1) - (self.G[:,self.z] * alpha).reshape(-1,1)
+                self.r = np.squeeze(self.r)
             else:
                 Aux = self.G[:,self.z] @ alpha
                 self.r = np.squeeze(self.b - Aux.T)
@@ -150,8 +198,18 @@ def Calculate(self):
                 ERROR_GLO = np.c_[ ERROR_GLO , self.nerrorACTUAL]
                 NPOINTS = np.c_[ NPOINTS , np.size(self.z)]
 
+            MaximumLengthZ = max(MaximumLengthZ, np.size(self.z))
             k = k+1
 
+            if k>1000:
+                """
+                this means using the initial candidate set, it was impossible to obtain a set of positive weights.
+                Try again without constraints!!!
+                TODO: incorporate this into greater workflow
+                """
+                self.success = False
+                break
+
         self.w = alpha.T * np.sqrt(self.W[self.z]) #TODO FIXME cope with weights vectors different from 1
 
         print(f'Total number of iterations = {k}')
@@ -163,17 +221,16 @@ def Calculate(self):
             plt.ylabel('Error %')
             plt.show()
 
-
-
-    """
-    Method for the quick update of weights (self.w), whenever a negative weight is found
-    """
     def _UpdateWeightsInverse(self, A,Aast,a,xold):
+        """
+        Method for the quick update of weights (self.w), whenever a negative weight is found
+        """
         c = np.dot(A.T, a)
         d = np.dot(Aast, c).reshape(-1, 1)
         s = np.dot(a.T, a) - np.dot(c.T, d)
         aux1 = np.hstack([Aast + np.outer(d, d) / s, -d / s])
         if np.shape(-d.T / s)[1]==1:
+            s = s.reshape(1,-1)
             aux2 = np.squeeze(np.hstack([-d.T / s, 1 / s]))
         else:
             aux2 = np.hstack([np.squeeze(-d.T / s), 1 / s])
@@ -182,24 +239,20 @@ def _UpdateWeightsInverse(self, A,Aast,a,xold):
         x = np.vstack([(xold - d * v), v])
         return Bast, x
 
-
-
-    """
-    Method for the quick update of weights (self.w), whenever a negative weight is found
-    """
     def _MultiUpdateInverseHermitian(self, invH, neg_indexes):
+        """
+        Method for the quick update of weights (self.w), whenever a negative weight is found
+        """
         neg_indexes = np.sort(neg_indexes)
         for i in range(np.size(neg_indexes)):
             neg_index = neg_indexes[i] - i
             invH = self._UpdateInverseHermitian(invH, neg_index)
         return invH
 
-
-
-    """
-    Method for the quick update of weights (self.w), whenever a negative weight is found
-    """
     def _UpdateInverseHermitian(self, invH, neg_index):
+        """
+        Method for the quick update of weights (self.w), whenever a negative weight is found
+        """
         if neg_index == np.shape(invH)[1]:
             aux = (invH[0:-1, -1] * invH[-1, 0:-1]) / invH(-1, -1)
             invH_new = invH[:-1, :-1] - aux