[MappingApplication] Proper Implementation of a Radial Basis Function Mapper for FEM-FEM and IBRA-FEM

applications/MappingApplication/custom_searching/interface_communicator.cpp

@@ -75,7 +75,9 @@ InterfaceCommunicator::InterfaceCommunicator(
         "max_num_search_iterations"     : 0,
         "print_bounding_boxes_to_file"  : false,
         "bounding_boxes_file_path"      : "",
-        "echo_level"                    : 0
+        "echo_level"                    : 0,
+        "required_rbf_support_points"  : 0,
+        "use_all_rbf_support_points": true
     })");
     // deliberately only validating defaults, but not assigning, since computing the defaults is expensive
     // see "ExchangeInterfaceData"

applications/MappingApplication/custom_utilities/mapper_interface_info.h

@@ -25,6 +25,8 @@
 
 namespace Kratos
 {
+class RBFSupportAccumulator;
+
 ///@addtogroup MappingApplication
 ///@{
 
@@ -118,6 +120,15 @@ class MapperInterfaceInfo
     virtual MapperInterfaceInfo::Pointer Create(const CoordinatesArrayType& rCoordinates,
                                                 const IndexType SourceLocalSystemIndex,
                                                 const IndexType SourceRank) const = 0;
+
+    /**
+    * @brief Returning the results of the neighbour search for the RBF mapper
+    * @author Juan I. Camarotti
+    */
+    virtual const RBFSupportAccumulator& GetRBFSupportAccumulator() const
+    {
+        KRATOS_ERROR << "Base class function called!" << std::endl;
+    }
 
     // needed for serialization
     virtual MapperInterfaceInfo::Pointer Create() const = 0;

applications/MappingApplication/custom_utilities/mapper_local_system.h

@@ -133,7 +133,7 @@ class MapperLocalSystem
     virtual CoordinatesArrayType& Coordinates() const = 0;
 
 
-    void AddInterfaceInfo(MapperInterfaceInfoPointerType pInterfaceInfo) // TODO pass by const ref?
+    virtual void AddInterfaceInfo(MapperInterfaceInfoPointerType pInterfaceInfo) // TODO pass by const ref?
     {
         mInterfaceInfos.push_back(pInterfaceInfo);
     }

applications/MappingApplication/custom_utilities/mapper_utilities.cpp

@@ -51,6 +51,26 @@ void AssignInterfaceEquationIds(Communicator& rModelPartCommunicator)
     rModelPartCommunicator.SynchronizeNonHistoricalVariable(INTERFACE_EQUATION_ID);
 }
 
+void AssignInterfaceEquationIdsOnConditions(Communicator& rModelPartCommunicator)
+{
+    if (rModelPartCommunicator.GetDataCommunicator().IsNullOnThisRank()) {
+        return;
+    }
+
+    const int num_conditions_local = rModelPartCommunicator.LocalMesh().NumberOfConditions();
+    int num_conditions_accumulated = rModelPartCommunicator.GetDataCommunicator().ScanSum(num_conditions_local);
+    const int start_equation_id = num_conditions_accumulated - num_conditions_local;
+    const auto conditions_begin = rModelPartCommunicator.LocalMesh().ConditionsBegin();
+
+    IndexPartition<unsigned int>(num_conditions_local).for_each(
+        [conditions_begin, start_equation_id](unsigned int i){
+            (conditions_begin + i)->SetValue(INTERFACE_EQUATION_ID, start_equation_id + i);
+        }
+    );
+
+    rModelPartCommunicator.SynchronizeNonHistoricalVariable(INTERFACE_EQUATION_ID);
+}
+
 template <typename TContainer>
 double ComputeMaxEdgeLengthLocal(const TContainer& rEntityContainer)
 {

applications/MappingApplication/custom_utilities/mapper_utilities.h

@@ -195,6 +195,16 @@ void UpdateModelPartFromSystemVector(
 */
 void KRATOS_API(MAPPING_APPLICATION) AssignInterfaceEquationIds(Communicator& rModelPartCommunicator);
 
+/**
+* @brief Assigning INTERFACE_EQUATION_IDs to the conditions, with and without MPI
+* This function assigns the INTERFACE_EQUATION_IDs to the conditions, which
+* act as EquationIds for the MappingMatrix in the case of IGA. This work with and without MPI,
+* in MPI a ScanSum is performed with the local number of nodes
+* @param rModelPartCommunicator The Modelpart-Communicator to be used
+* @author Juan I. Camarotti
+*/
+void KRATOS_API(MAPPING_APPLICATION) AssignInterfaceEquationIdsOnConditions(Communicator& rModelPartCommunicator);
+
 void KRATOS_API(MAPPING_APPLICATION) CreateMapperLocalSystemsFromNodes(const MapperLocalSystem& rMapperLocalSystemPrototype,
                                        const Communicator& rModelPartCommunicator,
                                        std::vector<Kratos::unique_ptr<MapperLocalSystem>>& rLocalSystems);

applications/MappingApplication/custom_utilities/mapping_matrix_utilities.cpp

@@ -206,12 +206,59 @@ void MappingMatrixUtilitiesType::BuildMappingMatrix(
 
     BuildMatrix(rpMappingMatrix, rMapperLocalSystems);
 
-    // refactor to be used from the mapper directly
-    // if (EchoLevel > 2) {
-    //     const std::string base_file_name = "O_" + rModelPartOrigin.Name() + "__D_" + rModelPartDestination.Name() +".mm";
-    //     MappingSparseSpaceType::WriteMatrixMarketMatrix(("MappingMatrix_"+base_file_name).c_str(), *rpMappingMatrix, false);
-    //     CheckRowSum<MappingSparseSpaceType, DenseSpaceType>(*rpMappingMatrix, base_file_name);
-    // }
+    MappingMatrixUtilitiesType::InitializeSystemVector(rpInterfaceVectorOrigin, num_nodes_origin);
+    MappingMatrixUtilitiesType::InitializeSystemVector(rpInterfaceVectorDestination, num_nodes_destination);
+
+    KRATOS_CATCH("")
+}
+
+template<>
+void MappingMatrixUtilitiesType::BuildMappingMatrixRBFMapper(
+    Kratos::unique_ptr<typename MappingSparseSpaceType::MatrixType>& rpMappingMatrix,
+    Kratos::unique_ptr<typename MappingSparseSpaceType::VectorType>& rpInterfaceVectorOrigin,
+    Kratos::unique_ptr<typename MappingSparseSpaceType::VectorType>& rpInterfaceVectorDestination,
+    const ModelPart& rModelPartOrigin,
+    const ModelPart& rModelPartDestination,
+    std::vector<Kratos::unique_ptr<MapperLocalSystem>>& rMapperLocalSystems,
+    const IndexType NumberOfPolynomialTerms,
+    const bool BuildOriginInterpolationMatrix,
+    const bool OriginIsIga,
+    const int EchoLevel)
+{
+    KRATOS_TRY
+
+    static_assert(!MappingSparseSpaceType::IsDistributed(), "Using a distributed Space!");
+
+    const SizeType num_nodes_origin = rModelPartOrigin.NumberOfNodes();
+    const SizeType num_conditions_origin = rModelPartOrigin.NumberOfConditions();
+    const SizeType num_nodes_destination = rModelPartDestination.NumberOfNodes();
+    const SizeType num_conditions_destination = rModelPartDestination.NumberOfConditions();
+
+    IndexType origin_size;
+    IndexType destination_size;
+
+    if (!OriginIsIga){
+        origin_size = num_nodes_origin;
+        destination_size = num_nodes_destination;
+    } else if (OriginIsIga && BuildOriginInterpolationMatrix) {
+        origin_size = num_conditions_origin;
+        destination_size = num_conditions_destination;
+    } else if (OriginIsIga && !BuildOriginInterpolationMatrix) {
+        origin_size = num_conditions_origin;
+        destination_size = num_nodes_destination;
+    }
+
+    // Initialize the Matrix
+    // This has to be done always since the Graph has changed if the Interface is updated!
+    if (BuildOriginInterpolationMatrix){
+        ConstructMatrixStructure(rpMappingMatrix, rMapperLocalSystems,
+                                origin_size + NumberOfPolynomialTerms, destination_size + NumberOfPolynomialTerms);
+    } else {
+        ConstructMatrixStructure(rpMappingMatrix, rMapperLocalSystems,
+                                origin_size + NumberOfPolynomialTerms, destination_size);
+    }
+
+    BuildMatrix(rpMappingMatrix, rMapperLocalSystems);
 
     MappingMatrixUtilitiesType::InitializeSystemVector(rpInterfaceVectorOrigin, num_nodes_origin);
     MappingMatrixUtilitiesType::InitializeSystemVector(rpInterfaceVectorDestination, num_nodes_destination);

applications/MappingApplication/custom_utilities/mapping_matrix_utilities.h

@@ -42,6 +42,18 @@ static void BuildMappingMatrix(
     std::vector<Kratos::unique_ptr<MapperLocalSystem>>& rMapperLocalSystems,
     const int EchoLevel);
 
+static void BuildMappingMatrixRBFMapper(
+    Kratos::unique_ptr<typename TSparseSpace::MatrixType>& rpMappingMatrix,
+    Kratos::unique_ptr<typename TSparseSpace::VectorType>& rpInterfaceVectorOrigin,
+    Kratos::unique_ptr<typename TSparseSpace::VectorType>& rpInterfaceVectorDestination,
+    const ModelPart& rModelPartOrigin,
+    const ModelPart& rModelPartDestination,
+    std::vector<Kratos::unique_ptr<MapperLocalSystem>>& rMapperLocalSystems,
+    const IndexType NumberOfPolynomialTerms,
+    const bool BuildOriginInterpolationMatrix,
+    const bool OriginIsIga,
+    const int EchoLevel);
+
 static void CheckRowSum(
     const typename TSparseSpace::MatrixType& rM,
     const std::string& rBaseFileName,

applications/MappingApplication/custom_utilities/radial_basis_functions_utilities.cpp

@@ -21,7 +21,7 @@ namespace Kratos
 namespace RadialBasisFunctionsUtilities
 {
     // Extracted from: Review of coupling methods for non-matching meshes (https://doi.org/10.1016/j.cma.2006.03.017)
-    double CalculateWendlandC2SupportRadius(const Matrix& rPoints, const double k = 2.5)
+    double CalculateWendlandC2SupportRadius(const Matrix& rPoints, const double k)
     {
         const SizeType n_points = rPoints.size1();
         KRATOS_ERROR_IF(n_points < 2) << "At least two points are required to estimate spacing." << std::endl;
@@ -42,4 +42,4 @@ namespace RadialBasisFunctionsUtilities
     }
 
 } // namespace RadialBasisFunctionsUtilities
-} // namespace Kratos.
+} // namespace Kratos.

applications/MappingApplication/custom_utilities/radial_basis_functions_utilities.h

@@ -81,8 +81,8 @@ struct WendlandC2 {
     }
 };
 
-double KRATOS_API(MAPPING_APPLICATION) CalculateWendlandC2SupportRadius(const Matrix& rPoints, const double k);
+double KRATOS_API(MAPPING_APPLICATION) CalculateWendlandC2SupportRadius(const Matrix& rPoints, const double k = 2.5);
 
 }  // namespace RadialBasisFunctionsUtilities.
 
-}  // namespace Kratos.
+}  // namespace Kratos.

applications/MappingApplication/mapping_application.cpp

@@ -38,6 +38,7 @@
 #include "custom_mappers/projection_3D_2D_mapper.h"
 #include "custom_mappers/beam_mapper.h"
 #include "custom_mappers/coupling_geometry_mapper.h"
+#include "custom_mappers/radial_basis_function_mapper.h"
 #include "custom_mappers/nearest_neighbor_mapper_iga.h"
 
 // Macros for registering mappers
@@ -95,6 +96,7 @@ void KratosMappingApplication::Register()
 
     KRATOS_REGISTER_MAPPER_WITH_BACKEND(NearestNeighborMapperIGA, "nearest_neighbor_iga");
 
+    KRATOS_REGISTER_MAPPER(RadialBasisFunctionMapper,  "radial_basis_function");
     KRATOS_REGISTER_MAPPER(CouplingGeometryMapper,  "coupling_geometry");
     KRATOS_REGISTER_MAPPER(BeamMapper,  "beam_mapper");
 

applications/MappingApplication/tests/cpp_tests/test_radial_basis_functions_utilities.cpp

@@ -121,4 +121,4 @@ KRATOS_TEST_CASE_IN_SUITE(RadialBasisFunctionsUtilities_WendlandC2SupportRadius,
     KRATOS_EXPECT_NEAR(support_radius, expected_radius, 1e-12);
 }
 
-} // namespace Kratos::Testing
+}

applications/MappingApplication/tests/test_coupling_geometry_mapper.py

@@ -3,6 +3,7 @@
 from KratosMultiphysics import KratosUnittest
 from KratosMultiphysics.testing.utilities import ReadModelPart
 import os
+from KratosMultiphysics.kratos_utilities import CheckIfApplicationsAvailable
 
 def GetFilePath(file_name):
     return os.path.join(os.path.dirname(os.path.realpath(__file__)), "mdpa_files", file_name)
@@ -30,10 +31,11 @@ def test_inverse_map_forces(self):
         reference_result = [0.2380991480071958, 0.2380991480071958, 0.2380991480071958, 1.3120351229689677, 1.3120351229689677, 1.3120351229689677, 0.6908309106360845, 0.6908309106360845, 0.6908309106360845, 0.9063686826513201, 0.9063686826513201, 0.9063686826513201, 0.9261336708771284, 0.9261336708771284, 0.9261336708771284, 0.9265324648593039, 0.9265324648593039, 0.9265324648593039]
         self.assertVectorAlmostEqual(mapped_results,reference_result)
 
-@KratosUnittest.skipIfApplicationsNotAvailable("IgaApplication")
 class TestIgaFEMCouplingGeometryMapper(KratosUnittest.TestCase):
     @classmethod
     def setUpClass(self):
+        if not CheckIfApplicationsAvailable("IgaApplication"):
+            raise KratosUnittest.SkipTest("The IgaApplication is not available!")
         import KratosMultiphysics.IgaApplication as Iga
         self.mapper_parameters = KM.Parameters("""{
             "mapper_type": "coupling_geometry",
@@ -63,10 +65,11 @@ def test_map_displacements(self):
         reference_result = [1.0, 1.0, 1.0, 0.9999999999999998, 0.9999999999999998, 0.9999999999999998]
         self.assertVectorAlmostEqual(mapped_results,reference_result)
 
-@KratosUnittest.skipIfApplicationsNotAvailable("IgaApplication")
 class TestDualMortarIgaFEMCouplingGeometryMapper(KratosUnittest.TestCase):
     @classmethod
     def setUpClass(self):
+        if not CheckIfApplicationsAvailable("IgaApplication"):
+            raise KratosUnittest.SkipTest("The IgaApplication is not available!")
         import KratosMultiphysics.IgaApplication as Iga
         self.mapper_parameters = KM.Parameters("""{
             "mapper_type": "coupling_geometry",

applications/MappingApplication/tests/test_nearest_element_mapper.py

@@ -5,6 +5,7 @@
 import blade_mapping_test
 import quadratic_mapper_tests
 from pathlib import Path
+from KratosMultiphysics.kratos_utilities import CheckIfApplicationsAvailable
 
 def GetFilePath(file_name):
     return Path(__file__).resolve().parent / "mdpa_files" / file_name
@@ -163,10 +164,11 @@ def setUpClass(cls):
         }""")
         super().setUpMapper(mapper_params, switch_sides=True)
 
-@KratosUnittest.skipIfApplicationsNotAvailable("IgaApplication")
 class BasicTestsLineMappingIGAFEM(basic_mapper_tests.BasicMapperTests):
     @classmethod
     def setUpClass(cls):
+        if not CheckIfApplicationsAvailable("IgaApplication"):
+            raise KratosUnittest.SkipTest("The IgaApplication is not available!")
         import KratosMultiphysics.IgaApplication as Iga
         mapper_params = KM.Parameters("""{
             "mapper_type": "nearest_element",
@@ -310,6 +312,8 @@ def test_Map_USE_TRANSPOSE_constant_scalar_both_non_historical(self):
 class BasicTestsSurfaceMappingIGAFEM(basic_mapper_tests.BasicMapperTests):
     @classmethod
     def setUpClass(cls):
+        if not CheckIfApplicationsAvailable("IgaApplication"):
+            raise KratosUnittest.SkipTest("The IgaApplication is not available!")
         import KratosMultiphysics.IgaApplication as Iga
         mapper_params = KM.Parameters("""{
             "mapper_type": "nearest_element",

applications/MappingApplication/tests/test_nearest_neighbor_mapper_iga.py

@@ -6,6 +6,7 @@
 data_comm = KM.Testing.GetDefaultDataCommunicator()
 if data_comm.IsDistributed():
     from KratosMultiphysics.MappingApplication import MPIExtension as MappingMPIExtension
+from KratosMultiphysics.kratos_utilities import CheckIfApplicationsAvailable
 
 # Additional imports
 import mapper_test_case
@@ -16,10 +17,11 @@
 def GetFilePath(file_name):
     return Path(__file__).resolve().parent / "mdpa_files" / file_name
 
-@KratosUnittest.skipIfApplicationsNotAvailable("IgaApplication")
 class BasicTestsLineMappingIGAFEM(basic_mapper_tests.BasicMapperTests):
     @classmethod
     def setUpClass(cls):
+        if not CheckIfApplicationsAvailable("IgaApplication"):
+            raise KratosUnittest.SkipTest("The IgaApplication is not available!")
         import KratosMultiphysics.IgaApplication as Iga
         mapper_params = KM.Parameters("""{
             "mapper_type": "nearest_neighbor_iga",
@@ -161,10 +163,11 @@ def test_Map_USE_TRANSPOSE_constant_scalar_FROM_NON_HISTORICAL(self):
     def test_Map_USE_TRANSPOSE_constant_scalar_both_non_historical(self):
         self.skipTest("Not implemented for this mapper")
 
-@KratosUnittest.skipIfApplicationsNotAvailable("IgaApplication")
 class BasicTestsSurfaceMappingIGAFEM(basic_mapper_tests.BasicMapperTests):
     @classmethod
     def setUpClass(cls):
+        if not CheckIfApplicationsAvailable("IgaApplication"):
+            raise KratosUnittest.SkipTest("The IgaApplication is not available!")
         import KratosMultiphysics.IgaApplication as Iga
         mapper_params = KM.Parameters("""{
             "mapper_type": "nearest_neighbor_iga",