Merge pull request #13032 from KratosMultiphysics/cut_domain_auxiliar…

…_utility

[FluidDynamicsApplication] Add a new utility to measure the volume of water/air in the cut elements.

applications/FluidDynamicsApplication/custom_python/add_custom_utilities_to_python.cpp

@@ -185,6 +185,8 @@ void  AddCustomUtilitiesToPython(pybind11::module& m)
         .def_static("CalculateFluidVolume", &FluidAuxiliaryUtilities::CalculateFluidVolume)
         .def_static("CalculateFluidPositiveVolume", &FluidAuxiliaryUtilities::CalculateFluidPositiveVolume)
         .def_static("CalculateFluidNegativeVolume", &FluidAuxiliaryUtilities::CalculateFluidNegativeVolume)
+        .def_static("CalculateFluidCutElementsNegativeVolume", &FluidAuxiliaryUtilities::CalculateFluidCutElementsNegativeVolume)
+        .def_static("CalculateFluidCutElementsPositiveVolume", &FluidAuxiliaryUtilities::CalculateFluidCutElementsPositiveVolume)
         .def_static("MapVelocityFromSkinToVolumeRBF", &FluidAuxiliaryUtilities::MapVelocityFromSkinToVolumeRBF)
         .def_static("FindMaximumEdgeLength", [](ModelPart& rModelPart){return FluidAuxiliaryUtilities::FindMaximumEdgeLength(rModelPart);})
         .def_static("FindMaximumEdgeLength", [](ModelPart& rModelPart, const bool CalculateNodalNeighbours){return FluidAuxiliaryUtilities::FindMaximumEdgeLength(rModelPart, CalculateNodalNeighbours);})

applications/FluidDynamicsApplication/custom_utilities/fluid_auxiliary_utilities.cpp

@@ -50,6 +50,96 @@ double FluidAuxiliaryUtilities::CalculateFluidVolume(const ModelPart& rModelPart
     return fluid_volume;
 }
 
+
+double FluidAuxiliaryUtilities::CalculateFluidCutElementsPositiveVolume(const ModelPart& rModelPart)
+{
+
+    // Check that there are elements and distance variable in the nodal database
+    KRATOS_ERROR_IF(rModelPart.GetCommunicator().GlobalNumberOfElements() == 0) << "There are no elements in the provided model part. Fluid cut elements volume cannot be computed." << std::endl;
+    const auto &r_communicator = rModelPart.GetCommunicator();
+    if (r_communicator.LocalMesh().NumberOfNodes() != 0)
+    {
+        KRATOS_ERROR_IF_NOT(r_communicator.LocalMesh().NodesBegin()->SolutionStepsDataHas(DISTANCE)) << "Nodal solution step data has no \'DISTANCE\' variable. Positive cut elements volume cannot be computed" << std::endl;
+    }
+
+    double cut_positive_fluid_volume = 0.0;
+    if (r_communicator.LocalMesh().NumberOfElements() != 0)
+    {
+        // Set the modified shape functions fatory
+        const auto &r_geom_begin = r_communicator.LocalMesh().ElementsBegin()->GetGeometry();
+        auto mod_sh_func_factory = GetStandardModifiedShapeFunctionsFactory(r_geom_begin);
+
+        // Calculate the positive volume
+        Vector nodal_distances(r_geom_begin.PointsNumber());
+        cut_positive_fluid_volume = block_for_each<SumReduction<double>>(rModelPart.GetCommunicator().LocalMesh().Elements(), nodal_distances, [&mod_sh_func_factory](Element &rElement, Vector &rNodalDistancesTLS)
+                                                                         {
+            // Set the distances vector to check if the element is split
+            const auto& r_geom = rElement.GetGeometry();
+            const std::size_t n_nodes = r_geom.PointsNumber();
+            for (std::size_t i_node = 0; i_node < n_nodes; ++i_node) {
+                rNodalDistancesTLS[i_node] = r_geom[i_node].FastGetSolutionStepValue(DISTANCE);
+            }
+            // Split element check
+            double elem_volume = 0.0;
+            if (IsSplit(rNodalDistancesTLS)) {
+                // Compute positive volume fraction with the modified shape functions
+                auto p_mod_sh_func = mod_sh_func_factory(rElement.pGetGeometry(), rNodalDistancesTLS);
+                elem_volume = p_mod_sh_func->ComputePositiveSideDomainSize();
+            }
+
+            // Return the value to be reduced
+            return elem_volume; });
+    }
+
+    // Synchronize among processors
+    r_communicator.GetDataCommunicator().SumAll(cut_positive_fluid_volume);
+
+    return cut_positive_fluid_volume;
+}
+
+double FluidAuxiliaryUtilities::CalculateFluidCutElementsNegativeVolume(const ModelPart& rModelPart)
+{
+    // Check that there are elements and distance variable in the nodal database
+    KRATOS_ERROR_IF(rModelPart.GetCommunicator().GlobalNumberOfElements() == 0) << "There are no elements in the provided model part. Fluid cut elements volume cannot be computed." << std::endl;
+    const auto &r_communicator = rModelPart.GetCommunicator();
+    if (r_communicator.LocalMesh().NumberOfNodes() != 0)
+    {
+        KRATOS_ERROR_IF_NOT(r_communicator.LocalMesh().NodesBegin()->SolutionStepsDataHas(DISTANCE)) << "Nodal solution step data has no \'DISTANCE\' variable. Positive cut elements volume cannot be computed" << std::endl;
+    }
+
+    double cut_negative_fluid_volume = 0.0;
+    if (r_communicator.LocalMesh().NumberOfElements() != 0)
+    {
+        // Set the modified shape functions fatory
+        const auto &r_geom_begin = r_communicator.LocalMesh().ElementsBegin()->GetGeometry();
+        auto mod_sh_func_factory = GetStandardModifiedShapeFunctionsFactory(r_geom_begin);
+
+        // Calculate the positive volume
+        Vector nodal_distances(r_geom_begin.PointsNumber());
+        cut_negative_fluid_volume = block_for_each<SumReduction<double>>(rModelPart.GetCommunicator().LocalMesh().Elements(), nodal_distances, [&mod_sh_func_factory](Element &rElement, Vector &rNodalDistancesTLS)
+                                                                         {
+            // Set the distances vector to check if the element is split
+            const auto& r_geom = rElement.GetGeometry();
+            const std::size_t n_nodes = r_geom.PointsNumber();
+            for (std::size_t i_node = 0; i_node < n_nodes; ++i_node) {
+                rNodalDistancesTLS[i_node] = r_geom[i_node].FastGetSolutionStepValue(DISTANCE);
+            }
+            // Split element check
+            double elem_volume = 0.0;
+            if (IsSplit(rNodalDistancesTLS)) {
+                // Compute positive volume fraction with the modified shape functions
+                auto p_mod_sh_func = mod_sh_func_factory(rElement.pGetGeometry(), rNodalDistancesTLS);
+                elem_volume = p_mod_sh_func->ComputeNegativeSideDomainSize();
+            }
+            // Return the value to be reduced
+            return elem_volume; });
+    }
+
+    // Synchronize among processors
+    r_communicator.GetDataCommunicator().SumAll(cut_negative_fluid_volume);
+
+    return cut_negative_fluid_volume;
+}
 double FluidAuxiliaryUtilities::CalculateFluidPositiveVolume(const ModelPart& rModelPart)
 {
     // Check that there are elements and distance variable in the nodal database

applications/FluidDynamicsApplication/custom_utilities/fluid_auxiliary_utilities.h

@@ -130,6 +130,30 @@ class KRATOS_API(FLUID_DYNAMICS_APPLICATION) FluidAuxiliaryUtilities
      */
     static double CalculateFluidPositiveVolume(const ModelPart& rModelPart);
 
+    /**
+     * @brief Calculate the fluid negative volume in cut elements
+     * For the given model part, this function calculates the total negative fluid volume fraction in the cut elements.
+     * It is assumed that a unique element geometry type it is present in the mesh.
+     * Only simplex geometries (linear triangle and tetrahedron) are supported.
+     * The negative fraction must be described in terms of a continuous level set function stored
+     * in the variable DISTANCE of the historical database
+     * @param rModelPart The model part to calculate the negative volume
+     * @return double Fluid negative volume over the cut elements
+     */
+    static double CalculateFluidCutElementsNegativeVolume(const ModelPart &rModelPart);
+
+    /**
+     * @brief Calculate the fluid positive volume in cut elements
+     * For the given model part, this function calculates the total positive fluid volume fraction in the cut elements.
+     * It is assumed that a unique element geometry type it is present in the mesh.
+     * Only simplex geometries (linear triangle and tetrahedron) are supported.
+     * The positive fraction must be described in terms of a continuous level set function stored
+     * in the variable DISTANCE of the historical database
+     * @param rModelPart The model part to calculate the positive volume
+     * @return double Fluid positive volume over the cut elements
+     */
+    static double CalculateFluidCutElementsPositiveVolume(const ModelPart &rModelPart);
+
     /**
      * @brief Calculate the fluid negative volume
      * For the given model part, this function calculates the total negative fluid volume fraction.

applications/FluidDynamicsApplication/tests/test_fluid_auxiliary_utilities.py

@@ -150,6 +150,28 @@ def testFindMaximumEdgeLength(self):
         max_edge_length = KratosFluid.FluidAuxiliaryUtilities.FindMaximumEdgeLength(fluid_model_part, calculate_nodal_neighbours)
         ref_max_edge_length = 0.2828427124746191
         self.assertAlmostEqual(max_edge_length, ref_max_edge_length, 12)
+
+    def testCalculateFluidNegativeCutVolume(self):
+        # Set fluid level set
+        level_set_y = 1.0/2.0
+        fluid_model_part = self.model.GetModelPart("FluidModelPart")
+        for node in fluid_model_part.Nodes:
+            node.SetSolutionStepValue(Kratos.DISTANCE, 0, node.Y - level_set_y)
+
+        # Calculate the fluid negative volume
+        fluid_negative_cut_volume = KratosFluid.FluidAuxiliaryUtilities.CalculateFluidCutElementsNegativeVolume(fluid_model_part)
+        self.assertAlmostEqual(fluid_negative_cut_volume,0.1, 12)
+
+    def testCalculateFluidPositiveCutVolume(self):
+        # Set fluid level set
+        level_set_y = 1.0/2.0
+        fluid_model_part = self.model.GetModelPart("FluidModelPart")
+        for node in fluid_model_part.Nodes:
+            node.SetSolutionStepValue(Kratos.DISTANCE, 0, node.Y - level_set_y)
+        # Calculate the fluid negative volume
+        fluid_positive_cut_volume = KratosFluid.FluidAuxiliaryUtilities.CalculateFluidCutElementsPositiveVolume(fluid_model_part)
+        print(fluid_positive_cut_volume)
+        self.assertAlmostEqual(fluid_positive_cut_volume, 0.1, 12)
 
     def tearDown(self):
         KratosUtils.DeleteFileIfExisting("Cavity/square5.time")