Merge pull request #11569 from KratosMultiphysics/Kratos_LSPG_HROM

Kratos lspg hrom

applications/RomApplication/custom_python/add_custom_utilities_to_python.cpp

@@ -43,10 +43,12 @@ void AddCustomUtilitiesToPython(pybind11::module& m)
     .def(init<ModelPart&, Parameters, BaseSchemeType::Pointer>()) //
     .def("GetProjectedResidualsOntoPhi",&RomResidualsUtility::GetProjectedResidualsOntoPhi) //
     .def("GetProjectedResidualsOntoPsi",&RomResidualsUtility::GetProjectedResidualsOntoPsi) //
+    .def("GetProjectedResidualsOntoJPhi",&RomResidualsUtility::GetProjectedResidualsOntoJPhi) //
     ;
 
     class_<RomAuxiliaryUtilities>(m, "RomAuxiliaryUtilities")
         .def_static("SetHRomComputingModelPart", &RomAuxiliaryUtilities::SetHRomComputingModelPart)
+        .def_static("SetHRomComputingModelPartWithNeighbours", &RomAuxiliaryUtilities::SetHRomComputingModelPartWithNeighbours)
         .def_static("SetHRomVolumetricVisualizationModelPart", &RomAuxiliaryUtilities::SetHRomVolumetricVisualizationModelPart)
         .def_static("GetHRomConditionParentsIds", &RomAuxiliaryUtilities::GetHRomConditionParentsIds)
         .def_static("GetNodalNeighbouringElementIdsNotInHRom", &RomAuxiliaryUtilities::GetNodalNeighbouringElementIdsNotInHRom)

applications/RomApplication/custom_utilities/rom_auxiliary_utilities.cpp

@@ -170,6 +170,170 @@ void RomAuxiliaryUtilities::RecursiveHRomModelPartCreation(
     }
 }
 
+void RomAuxiliaryUtilities::SetHRomComputingModelPartWithNeighbours(
+    const Parameters HRomWeights,
+    ModelPart& rOriginModelPart,
+    ModelPart& rHRomComputingModelPart)
+{
+    // Ensure that the provided destination model part is empty
+    rHRomComputingModelPart.Clear();
+
+    // Auxiliary containers to save the entities involved in the HROM mesh
+    NodesPointerSetType hrom_nodes_set;
+    std::vector<Element::Pointer> hrom_elems_vect;
+    std::vector<Condition::Pointer> hrom_conds_vect;
+
+    const auto& r_elem_weights = HRomWeights["Elements"];
+    const auto& r_cond_weights = HRomWeights["Conditions"];
+
+    hrom_elems_vect.reserve(rOriginModelPart.NumberOfElements());
+    hrom_conds_vect.reserve(rOriginModelPart.NumberOfConditions());
+
+    FindGlobalNodalEntityNeighboursProcess<ModelPart::ElementsContainerType> find_nodal_elements_neighbours_process(rOriginModelPart);
+    find_nodal_elements_neighbours_process.Execute();
+    FindGlobalNodalEntityNeighboursProcess<ModelPart::ConditionsContainerType> find_nodal_conditions_neighbours_process(rOriginModelPart);
+    find_nodal_conditions_neighbours_process.Execute();
+
+    for (auto it = r_elem_weights.begin(); it != r_elem_weights.end(); ++it) {
+        // Get element from origin mesh
+        const IndexType elem_id = stoi(it.name());
+        const auto p_elem = rOriginModelPart.pGetElement(elem_id + 1);
+
+        // Add the element to the auxiliary container and to the main HROM model part
+        if(std::find(hrom_elems_vect.begin(), hrom_elems_vect.end(), p_elem) == hrom_elems_vect.end()) {
+            hrom_elems_vect.push_back(p_elem);
+            rHRomComputingModelPart.AddElement(p_elem);
+        }
+
+        // Add the element nodes to the auxiliary set and to the main HROM model part
+        const auto& r_geom = p_elem->GetGeometry();
+        const SizeType n_nodes = r_geom.PointsNumber();
+        for (IndexType i_node = 0; i_node < n_nodes; ++i_node) {
+            NodeType::Pointer p_node = r_geom(i_node);
+            hrom_nodes_set.insert(hrom_nodes_set.end(), p_node);
+            rHRomComputingModelPart.AddNode(p_node);
+
+            // Get the neighbor elements from each node and add them to the model part
+            for (auto& neighbor_elem : p_node->GetValue(NEIGHBOUR_ELEMENTS)) {
+                Kratos::Element::Pointer p_neighbor_elem = &neighbor_elem;
+                // Be careful to not add an element that is already in the model part.
+                // For this, we will check if the element is already in our vector of elements.
+                if(std::find(hrom_elems_vect.begin(), hrom_elems_vect.end(), p_neighbor_elem) == hrom_elems_vect.end()) {
+                    hrom_elems_vect.push_back(p_neighbor_elem);
+                    rHRomComputingModelPart.AddElement(p_neighbor_elem);
+
+                    // Add the nodes of the neighboring element to the model part
+                    const auto& neighbor_r_geom = p_neighbor_elem->GetGeometry();
+                    const SizeType neighbor_n_nodes = neighbor_r_geom.PointsNumber();
+                    for (IndexType i_neighbor_node = 0; i_neighbor_node < neighbor_n_nodes; ++i_neighbor_node) {
+                        NodeType::Pointer p_neighbor_node = neighbor_r_geom(i_neighbor_node);
+                        hrom_nodes_set.insert(hrom_nodes_set.end(), p_neighbor_node);
+                        rHRomComputingModelPart.AddNode(p_neighbor_node);
+                    }
+                }
+            }
+
+            // Get the neighbor conditions from each node and add them to the model part
+            for (auto& neighbor_cond : p_node->GetValue(NEIGHBOUR_CONDITIONS)) {
+                Kratos::Condition::Pointer p_neighbor_cond = &neighbor_cond;
+                // Be careful to not add a condition that is already in the model part.
+                // For this, we will check if the condition is already in our vector of conditions.
+                if(std::find(hrom_conds_vect.begin(), hrom_conds_vect.end(), p_neighbor_cond) == hrom_conds_vect.end()) {
+                    hrom_conds_vect.push_back(p_neighbor_cond);
+                    rHRomComputingModelPart.AddCondition(p_neighbor_cond);
+
+                    // Add the nodes of the neighboring condition to the model part
+                    const auto& neighbor_r_geom = p_neighbor_cond->GetGeometry();
+                    const SizeType neighbor_n_nodes = neighbor_r_geom.PointsNumber();
+                    for (IndexType i_neighbor_node = 0; i_neighbor_node < neighbor_n_nodes; ++i_neighbor_node) {
+                        NodeType::Pointer p_neighbor_node = neighbor_r_geom(i_neighbor_node);
+                        hrom_nodes_set.insert(hrom_nodes_set.end(), p_neighbor_node);
+                        rHRomComputingModelPart.AddNode(p_neighbor_node);
+                    }
+                }
+            }
+        }
+    }
+
+    for (auto it = r_cond_weights.begin(); it != r_cond_weights.end(); ++it) {
+        // Get the condition from origin mesh
+        const IndexType cond_id = stoi(it.name());
+        auto p_cond = rOriginModelPart.pGetCondition(cond_id + 1);
+
+        // Add the condition to the auxiliary container and to the main HROM model part
+        if(std::find(hrom_conds_vect.begin(), hrom_conds_vect.end(), p_cond) == hrom_conds_vect.end()) {
+            hrom_conds_vect.push_back(p_cond);
+            rHRomComputingModelPart.AddCondition(p_cond);
+        }
+
+        // Add the condition nodes to the auxiliary set and to the main HROM model part
+        const auto& r_geom = p_cond->GetGeometry();
+        const SizeType n_nodes = r_geom.PointsNumber();
+        for (IndexType i_node = 0; i_node < n_nodes; ++i_node) {
+            auto p_node = r_geom(i_node);
+            hrom_nodes_set.insert(hrom_nodes_set.end(), p_node);
+            rHRomComputingModelPart.AddNode(p_node);
+
+            // Get the neighbor elements from each node and add them to the model part
+            for (auto& neighbor_elem : p_node->GetValue(NEIGHBOUR_ELEMENTS)) {
+                Kratos::Element::Pointer p_neighbor_elem = &neighbor_elem;
+                // Be careful to not add an element that is already in the model part.
+                // For this, we will check if the element is already in our vector of elements.
+                if(std::find(hrom_elems_vect.begin(), hrom_elems_vect.end(), p_neighbor_elem) == hrom_elems_vect.end()) {
+                    hrom_elems_vect.push_back(p_neighbor_elem);
+                    rHRomComputingModelPart.AddElement(p_neighbor_elem);
+
+                    // Add the nodes of the neighboring element to the model part
+                    const auto& neighbor_r_geom = p_neighbor_elem->GetGeometry();
+                    const SizeType neighbor_n_nodes = neighbor_r_geom.PointsNumber();
+                    for (IndexType i_neighbor_node = 0; i_neighbor_node < neighbor_n_nodes; ++i_neighbor_node) {
+                        NodeType::Pointer p_neighbor_node = neighbor_r_geom(i_neighbor_node);
+                        hrom_nodes_set.insert(hrom_nodes_set.end(), p_neighbor_node);
+                        rHRomComputingModelPart.AddNode(p_neighbor_node);
+                    }
+                }
+            }
+
+            // Get the neighbor conditions from each node and add them to the model part
+            for (auto& neighbor_cond : p_node->GetValue(NEIGHBOUR_CONDITIONS)) {
+                Kratos::Condition::Pointer p_neighbor_cond = &neighbor_cond;
+                // Be careful to not add a condition that is already in the model part.
+                // For this, we will check if the condition is already in our vector of conditions.
+                if(std::find(hrom_conds_vect.begin(), hrom_conds_vect.end(), p_neighbor_cond) == hrom_conds_vect.end()) {
+                    hrom_conds_vect.push_back(p_neighbor_cond);
+                    rHRomComputingModelPart.AddCondition(p_neighbor_cond);
+
+                    // Add the nodes of the neighboring condition to the model part
+                    const auto& neighbor_r_geom = p_neighbor_cond->GetGeometry();
+                    const SizeType neighbor_n_nodes = neighbor_r_geom.PointsNumber();
+                    for (IndexType i_neighbor_node = 0; i_neighbor_node < neighbor_n_nodes; ++i_neighbor_node) {
+                        NodeType::Pointer p_neighbor_node = neighbor_r_geom(i_neighbor_node);
+                        hrom_nodes_set.insert(hrom_nodes_set.end(), p_neighbor_node);
+                        rHRomComputingModelPart.AddNode(p_neighbor_node);
+                    }
+                }
+            }
+        }
+    }
+    hrom_elems_vect.shrink_to_fit();
+    hrom_conds_vect.shrink_to_fit();
+
+    //TODO: ADD MPC'S
+
+    // Add properties to the HROM mesh
+    // Note that we add all the properties although some of them might note be used in the HROM mesh
+    auto& r_root_model_part = const_cast<ModelPart&>(rOriginModelPart).GetRootModelPart();
+    auto& r_properties = r_root_model_part.rProperties();
+    for (auto it_p_prop = r_properties.ptr_begin(); it_p_prop < r_properties.ptr_end(); ++it_p_prop) {
+        rHRomComputingModelPart.AddProperties(*it_p_prop);
+    }
+
+    // Create and fill the HROM calculation sub model parts
+    for (auto& r_orig_sub_mp : rOriginModelPart.SubModelParts()) {
+        RecursiveHRomModelPartCreation(hrom_nodes_set, hrom_elems_vect, hrom_conds_vect, r_orig_sub_mp, rHRomComputingModelPart);
+    }
+}
+
 //TODO: Make it thin walled and beam compatible
 void RomAuxiliaryUtilities::SetHRomVolumetricVisualizationModelPart(
     const ModelPart& rOriginModelPart,
@@ -615,4 +779,16 @@ void RomAuxiliaryUtilities::GetPsiElemental(
         }
     }
 
+void RomAuxiliaryUtilities::GetJPhiElemental(
+    Matrix &rJPhiElemental,
+    const Element::DofsVectorType& rDofs,
+    const Matrix &rJPhi)
+    {
+        for(std::size_t i = 0; i < rDofs.size(); ++i)
+        {
+            const Dof<double>& r_dof = *rDofs[i];
+            noalias(row(rJPhiElemental, i)) = row(rJPhi, r_dof.EquationId());
+        }
+    }
+
 } // namespace Kratos

applications/RomApplication/custom_utilities/rom_auxiliary_utilities.h

@@ -82,6 +82,24 @@ class KRATOS_API(ROM_APPLICATION) RomAuxiliaryUtilities
         const ModelPart& rOriginModelPart,
         ModelPart& rHRomComputingModelPart);
 
+    /**
+     * @brief Sets the HROM model part including neighboring entities based on the nodal weights
+     *
+     * Provided an origin model part and a parameters object containing the HROM weights,
+     * this method stores not only the elements and conditions required by the HROM in the destination
+     * model part, but also includes neighboring elements and conditions of the nodes involved.
+     * The resulting model part features the same submodelpart hierarchy that the origin has, but is
+     * augmented with the neighboring elements and conditions.
+     *
+     * @param HRomWeights Parameters object containing the HROM elemental and condition weights
+     * @param rOriginModelPart Origin model part (this is likely the computing model part)
+     * @param rHRomComputingModelPart Destination model part to store the HROM mesh augmented with neighbours
+     */
+    static void SetHRomComputingModelPartWithNeighbours(
+        const Parameters HRomWeights,
+        ModelPart& rOriginModelPart,
+        ModelPart& rHRomComputingModelPart);
+
     /**
      * @brief Sets the HROM skin visualization model part for a volumetric body
      * This function detects the skin of the origin modelpart and creates the corresponding skin
@@ -247,6 +265,21 @@ class KRATOS_API(ROM_APPLICATION) RomAuxiliaryUtilities
         const Element::GeometryType& rGeom,
         const std::unordered_map<Kratos::VariableData::KeyType, Matrix::size_type>& rVarToRowMapping);
 
+    /**
+     * @brief Obtain the JPhi elemental matrix for a particular element. 
+     * JPhi represents the projection of the Jacobian onto the ROM_BASIS.
+     * @param rJPhiElemental The matrix to store the result in. Must have the appropriate size already.
+     * @param rDofs The set of degrees of freedom (DoFs) of the element.
+     * @param rJPhi The JPhi matrix, from which rows are extracted according to the equation ID of each DoF.
+     * 
+     * This function loops over all the DoFs for the given element. For each DoF, it uses its equation ID to extract a 
+     * corresponding row from the rJPhi matrix, which is then stored in the corresponding row of rJPhiElemental.
+     */
+    static void GetJPhiElemental(
+        Matrix &rJPhiElemental,
+        const Element::DofsVectorType& rDofs,
+        const Matrix &rJPhi);
+
     ///@}
 
     private:

applications/RomApplication/custom_utilities/rom_residuals_utility.h

@@ -72,6 +72,18 @@ namespace Kratos
 
         ~RomResidualsUtility()= default;
 
+        /**
+         * Resizes a Matrix if it's not the right size
+         */
+        template<typename TMatrix>
+        static void ResizeIfNeeded(TMatrix& rMat, const std::size_t Rows, const std::size_t Cols)
+
+        {
+            if(rMat.size1() != Rows || rMat.size2() != Cols) {
+                rMat.resize(Rows, Cols, false);
+            }
+        };
+
         Matrix GetProjectedResidualsOntoPhi()
         {
             // Getting the number of elements and conditions from the model
@@ -208,6 +220,70 @@ namespace Kratos
         return matrix_residuals;
         }
 
+        Matrix GetProjectedResidualsOntoJPhi(
+            Matrix& rJPhi
+        )
+        {
+            const int n_elements = static_cast<int>(mrModelPart.Elements().size());
+            const int n_conditions = static_cast<int>(mrModelPart.Conditions().size());
+
+            const auto& r_current_process_info = mrModelPart.GetProcessInfo();
+
+            //contributions to the system
+            Vector rhs_contribution;
+
+            //vector containing the localization in the system of the different terms
+            Element::EquationIdVectorType equation_id;
+            Matrix matrix_residuals( (n_elements + n_conditions), mRomDofs);
+            Matrix phi_j_elemental;
+
+            const auto el_begin = mrModelPart.ElementsBegin();
+            const auto cond_begin = mrModelPart.ConditionsBegin();
+
+            //dofs container initialization
+            Element::DofsVectorType elem_dofs;
+            Condition::DofsVectorType cond_dofs;
+            #pragma omp parallel firstprivate(n_elements, n_conditions, rhs_contribution, equation_id, phi_j_elemental, el_begin, cond_begin, elem_dofs, cond_dofs)
+            {
+                #pragma omp for
+                for (int k = 0; k < n_elements; k++){
+                    auto r_element = el_begin + k;
+                    if (r_element->IsDefined(ACTIVE) && r_element->IsNot(ACTIVE)) continue;
+
+                    mpScheme->CalculateRHSContribution(*r_element, rhs_contribution, equation_id, r_current_process_info);
+                    r_element->GetDofList(elem_dofs, r_current_process_info);
+
+                    const std::size_t ndofs = elem_dofs.size();
+                    ResizeIfNeeded(phi_j_elemental, ndofs, mRomDofs);
+                    RomAuxiliaryUtilities::GetJPhiElemental(phi_j_elemental, elem_dofs, rJPhi);
+
+                    #pragma omp critical
+                    {
+                        noalias(row(matrix_residuals, k)) = prod(trans(phi_j_elemental), rhs_contribution);
+                    }
+                }
+
+                #pragma omp for
+                for (int k = 0; k < n_conditions; k++){
+                    auto r_condition = cond_begin + k;
+                    if (r_condition->IsDefined(ACTIVE) && r_condition->IsNot(ACTIVE)) continue;
+
+                    mpScheme->CalculateRHSContribution(*r_condition, rhs_contribution, equation_id, r_current_process_info);
+                    r_condition->GetDofList(cond_dofs, r_current_process_info);
+
+                    const std::size_t ndofs = cond_dofs.size();
+                    ResizeIfNeeded(phi_j_elemental, ndofs, mRomDofs);
+                    RomAuxiliaryUtilities::GetJPhiElemental(phi_j_elemental, cond_dofs, rJPhi);
+
+                    #pragma omp critical
+                    {
+                        noalias(row(matrix_residuals, n_elements + k)) = prod(trans(phi_j_elemental), rhs_contribution);
+                    }
+                }
+            }
+            return matrix_residuals;
+        }
+
     protected:
         std::vector< std::string > mNodalVariablesNames;
         int mNodalDofs;