Merge branch 'reset_ip_data' into 'master'

Reset integration point data by materialIDs See merge request ogs/ogs!4668
ufz · Jul 25, 2023 · 17cbb68 · 17cbb68
2 parents 7d33b1c + 58229e2
commit 17cbb68
Show file tree

Hide file tree

Showing 17 changed files with 406 additions and 21 deletions.
diff --git a/Applications/ApplicationsLib/CMakeLists.txt b/Applications/ApplicationsLib/CMakeLists.txt
@@ -9,7 +9,7 @@ target_link_libraries(
     ApplicationsLib
     PUBLIC BaseLib GeoLib NumLib
            $<$<TARGET_EXISTS:VTK::ParallelMPI>:VTK::ParallelMPI>
-    PRIVATE CMakeInfoLib MathLib MeshLib MeshGeoToolsLib pybind11::pybind11
+    PRIVATE CMakeInfoLib MathLib MeshLib MeshGeoToolsLib MeshToolsLib pybind11::pybind11
             $<$<BOOL:${OGS_USE_PETSC}>:petsc>
 )
 

diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
@@ -46,6 +46,7 @@
 #include "MeshGeoToolsLib/SearchLength.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/Utils/SetMeshSpaceDimension.h"
+#include "MeshToolsLib/ZeroMeshFieldDataByMaterialIDs.h"
 #include "NumLib/ODESolver/ConvergenceCriterion.h"
 #include "ProcessLib/CreateJacobianAssembler.h"
 #include "ProcessLib/DeactivatedSubdomain.h"
@@ -250,6 +251,21 @@ std::vector<std::unique_ptr<MeshLib::Mesh>> readMeshes(
         }
     }
 
+    auto const zero_mesh_field_data_by_material_ids =
+        //! \ogs_file_param{prj__zero_mesh_field_data_by_material_ids}
+        config.getConfigParameterOptional<std::vector<int>>(
+            "zero_mesh_field_data_by_material_ids");
+    if (zero_mesh_field_data_by_material_ids)
+    {
+        WARN(
+            "Tag 'zero_mesh_field_data_by_material_ids` is experimental. Its "
+            "name may be changed, or it may be removed due to its "
+            "corresponding feature becomes a single tool. Please use it with "
+            "care!");
+        MeshToolsLib::zeroMeshFieldDataByMaterialIDs(
+            *meshes[0], *zero_mesh_field_data_by_material_ids);
+    }
+
     MeshLib::setMeshSpaceDimension(meshes);
 
     return meshes;

diff --git a/Applications/Utils/MeshEdit/PVTU2VTU/CMakeLists.txt b/Applications/Utils/MeshEdit/PVTU2VTU/CMakeLists.txt
@@ -1,7 +1,6 @@
 ogs_add_executable(
     pvtu2vtu
     PVTU2VTU.cpp
-    ${PROJECT_SOURCE_DIR}/Applications/Utils/ModelPreparation/PartitionMesh/IntegrationPointDataTools.cpp
 )
 target_link_libraries(pvtu2vtu GitInfoLib tclap MeshToolsLib NumLib)
 install(TARGETS pvtu2vtu RUNTIME DESTINATION bin)
diff --git a/Applications/Utils/MeshEdit/PVTU2VTU/PVTU2VTU.cpp b/Applications/Utils/MeshEdit/PVTU2VTU/PVTU2VTU.cpp
@@ -21,7 +21,6 @@
 #include <unordered_set>
 #include <vector>
 
-#include "Applications/Utils/ModelPreparation/PartitionMesh/IntegrationPointDataTools.h"
 #include "BaseLib/FileTools.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
@@ -32,6 +31,7 @@
 #include "MeshLib/Properties.h"
 #include "MeshLib/Utils/IntegrationPointWriter.h"
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
+#include "MeshToolsLib/IntegrationPointDataTools.h"
 
 struct MeshEntityMapInfo
 {
@@ -125,7 +125,7 @@ bool createPropertyVector(
         for (auto const& mesh : partitioned_meshes)
         {
             partition_element_offsets.emplace_back(
-                ApplicationUtils::getIntegrationPointDataOffsetsOfMeshElements(
+                MeshToolsLib::getIntegrationPointDataOffsetsOfMeshElements(
                     mesh->getElements(), *pv, properties));
         }
 
@@ -140,14 +140,14 @@ bool createPropertyVector(
         for (auto const element : merged_mesh.getElements())
         {
             int const number_of_integration_points =
-                ApplicationUtils::getNumberOfElementIntegrationPoints(
-                    ip_meta_data, *element);
+                MeshToolsLib::getNumberOfElementIntegrationPoints(ip_meta_data,
+                                                                  *element);
             counter += number_of_integration_points;
         }
         new_pv->resize(counter * pv_num_components);
 
         auto const global_ip_offsets =
-            ApplicationUtils::getIntegrationPointDataOffsetsOfMeshElements(
+            MeshToolsLib::getIntegrationPointDataOffsetsOfMeshElements(
                 merged_mesh.getElements(), *pv, properties);
 
         std::size_t element_counter = 0;

diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/CMakeLists.txt b/Applications/Utils/ModelPreparation/PartitionMesh/CMakeLists.txt
@@ -1,6 +1,6 @@
 if(OGS_USE_MPI)
     ogs_add_executable(
-        binaryToPVTU BinaryToPVTU.cpp IntegrationPointDataTools.cpp NodeWiseMeshPartitioner.cpp
+        binaryToPVTU BinaryToPVTU.cpp NodeWiseMeshPartitioner.cpp
     )
     target_link_libraries(
            binaryToPVTU GitInfoLib ProcessLib MeshLib MeshToolsLib tclap
@@ -10,7 +10,7 @@ if(OGS_USE_MPI)
 endif()
 
 ogs_add_executable(
-    partmesh PartitionMesh.cpp Metis.cpp IntegrationPointDataTools.cpp NodeWiseMeshPartitioner.cpp
+    partmesh PartitionMesh.cpp Metis.cpp NodeWiseMeshPartitioner.cpp
 )
 target_link_libraries(partmesh GitInfoLib MeshLib MeshToolsLib tclap NumLib)
 add_dependencies(partmesh mpmetis)

diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp b/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp
@@ -21,12 +21,12 @@
 #include "BaseLib/Error.h"
 #include "BaseLib/FileTools.h"
 #include "BaseLib/Logging.h"
-#include "IntegrationPointDataTools.h"
 #include "MeshLib/Elements/Elements.h"
 #include "MeshLib/IO/NodeData.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
 #include "MeshLib/MeshEnums.h"
 #include "MeshLib/Utils/IntegrationPointWriter.h"
+#include "MeshToolsLib//IntegrationPointDataTools.h"
 
 namespace ApplicationUtils
 {
@@ -383,7 +383,8 @@ std::size_t copyFieldPropertyDataToPartitions(
         for (auto const element : elements)
         {
             int const number_of_element_field_data =
-                getNumberOfElementIntegrationPoints(ip_meta_data, *element) *
+                MeshToolsLib::getNumberOfElementIntegrationPoints(ip_meta_data,
+                                                                  *element) *
                 n_components;
             // The original element ID is not changed.
             auto const element_id = element->getID();
@@ -431,7 +432,8 @@ void setIntegrationPointNumberOfPartition(MeshLib::Properties const& properties,
             for (auto const element : elements)
             {
                 int const number_of_integration_points =
-                    getNumberOfElementIntegrationPoints(ip_meta_data, *element);
+                    MeshToolsLib::getNumberOfElementIntegrationPoints(
+                        ip_meta_data, *element);
                 counter += number_of_integration_points;
             }
             return counter;
@@ -476,8 +478,9 @@ bool copyPropertyVector(
             partitioned_pv_size = pv->size() * partitions.size();
         }
 
-        element_ip_data_offsets = getIntegrationPointDataOffsetsOfMeshElements(
-            global_mesh_elements, *pv, properties);
+        element_ip_data_offsets =
+            MeshToolsLib::getIntegrationPointDataOffsetsOfMeshElements(
+                global_mesh_elements, *pv, properties);
     }
 
     auto partitioned_pv = partitioned_properties.createNewPropertyVector<T>(
@@ -652,7 +655,8 @@ void checkFieldPropertyVectorSize(
         for (auto const element : global_mesh_elements)
         {
             int const number_of_integration_points =
-                getNumberOfElementIntegrationPoints(ip_meta_data, *element);
+                MeshToolsLib::getNumberOfElementIntegrationPoints(ip_meta_data,
+                                                                  *element);
             number_of_total_integration_points += number_of_integration_points;
         }
 

diff --git a/Documentation/ProjectFile/prj/t_zero_mesh_field_data_by_material_ids.md b/Documentation/ProjectFile/prj/t_zero_mesh_field_data_by_material_ids.md
@@ -0,0 +1,6 @@
+It specifies material IDs for resetting the integration point data. For the
+ elements with the specified material IDs, their double-type integration point
+ data are set to zero.
+
+It can be used for a restarted simulation with reactivated elements, e.g.
+a back-filling simulation.
diff --git a/...rtitionMesh/IntegrationPointDataTools.cpp → MeshToolsLib/IntegrationPointDataTools.cpp b/...rtitionMesh/IntegrationPointDataTools.cpp → MeshToolsLib/IntegrationPointDataTools.cpp
@@ -18,7 +18,7 @@
 #include "MeshLib/Utils/IntegrationPointWriter.h"
 #include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
 
-namespace ApplicationUtils
+namespace MeshToolsLib
 {
 template <typename ElementType>
 int getNumberOfElementIntegrationPointsGeneral(
@@ -127,4 +127,4 @@ std::vector<std::size_t> getIntegrationPointDataOffsetsOfMeshElements(
     return element_ip_data_offsets;
 }
 
-}  // namespace ApplicationUtils
+}  // namespace MeshToolsLib
diff --git a/...PartitionMesh/IntegrationPointDataTools.h → MeshToolsLib/IntegrationPointDataTools.h b/...PartitionMesh/IntegrationPointDataTools.h → MeshToolsLib/IntegrationPointDataTools.h
@@ -21,9 +21,8 @@ class PropertyVectorBase;
 struct IntegrationPointMetaData;
 }  // namespace MeshLib
 
-namespace ApplicationUtils
+namespace MeshToolsLib
 {
-
 int getNumberOfElementIntegrationPoints(
     MeshLib::IntegrationPointMetaData const& ip_meta_data,
     MeshLib::Element const& e);
@@ -32,4 +31,4 @@ std::vector<std::size_t> getIntegrationPointDataOffsetsOfMeshElements(
     std::vector<MeshLib::Element*> const& mesh_elements,
     MeshLib::PropertyVectorBase const& pv,
     MeshLib::Properties const& properties);
-}  // namespace ApplicationUtils
+}  // namespace MeshToolsLib
diff --git a/MeshToolsLib/ZeroMeshFieldDataByMaterialIDs.cpp b/MeshToolsLib/ZeroMeshFieldDataByMaterialIDs.cpp
@@ -0,0 +1,100 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on July 6, 2023, 4:21 PM
+ */
+
+#include "ZeroMeshFieldDataByMaterialIDs.h"
+
+#include <algorithm>
+#include <functional>
+
+#include "IntegrationPointDataTools.h"
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/Utils/IntegrationPointWriter.h"
+
+namespace MeshToolsLib
+{
+void zeroMeshFieldDataByMaterialIDs(
+    MeshLib::Mesh& mesh, std::vector<int> const& selected_material_ids)
+{
+    auto const materialIds = materialIDs(mesh);
+    if (!materialIds)
+    {
+        OGS_FATAL(
+            "Mesh contains no int-property vector named 'MaterialIDs' needed "
+            "by the 'zero_mesh_field_data_by_materialIDs'");
+    }
+
+    std::vector<std::size_t> element_ids_for_selected_materials;
+    for (std::size_t i = 0; i < materialIds->size(); ++i)
+    {
+        if (std::find(selected_material_ids.begin(),
+                      selected_material_ids.end(),
+                      (*materialIds)[i]) != selected_material_ids.end())
+        {
+            element_ids_for_selected_materials.push_back(i);
+        }
+    }
+
+    MeshLib::Properties& properties = mesh.getProperties();
+
+    std::vector<std::size_t> element_ip_data_offsets;
+
+    for (auto const& [name, property] : properties)
+    {
+        if (auto const item_type = property->getMeshItemType();
+            item_type != MeshLib::MeshItemType::IntegrationPoint)
+        {
+            continue;
+        }
+
+        // For special field data such as OGS_VERSION,
+        // IntegrationPointMetaData,
+        // etc., which are not "real" integration points:
+        if (!property->getPropertyName().ends_with("_ip"))
+        {
+            continue;
+        }
+
+        if (properties.template hasPropertyVector<double>(
+                name, MeshLib::MeshItemType::IntegrationPoint))
+        {
+            auto& pv = *properties.template getPropertyVector<double>(name);
+            const int n_components = pv.getNumberOfGlobalComponents();
+
+            if (element_ip_data_offsets.empty())
+            {
+                // The returned values has already been multiplied with
+                // pv.getNumberOfGlobalComponents()
+                element_ip_data_offsets =
+                    MeshToolsLib::getIntegrationPointDataOffsetsOfMeshElements(
+                        mesh.getElements(), pv, properties);
+
+                // element_ip_data_offsets / pv.getNumberOfGlobalComponents()
+                std::transform(element_ip_data_offsets.begin(),
+                               element_ip_data_offsets.end(),
+                               element_ip_data_offsets.begin(),
+                               [n = n_components](double const v)
+                               { return v / n; });
+            }
+
+            for (auto const element_id : element_ids_for_selected_materials)
+            {
+                std::fill(
+                    pv.begin() +
+                        n_components * element_ip_data_offsets[element_id],
+                    pv.begin() +
+                        n_components * element_ip_data_offsets[element_id + 1],
+                    0.0);
+            }
+        }
+    }
+}
+}  // namespace MeshToolsLib
diff --git a/MeshToolsLib/ZeroMeshFieldDataByMaterialIDs.h b/MeshToolsLib/ZeroMeshFieldDataByMaterialIDs.h
@@ -0,0 +1,27 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on July 6, 2023, 4:21 PM
+ */
+
+#pragma once
+
+#include <vector>
+
+namespace MeshLib
+{
+class Mesh;
+
+}
+
+namespace MeshToolsLib
+{
+
+void zeroMeshFieldDataByMaterialIDs(
+    MeshLib::Mesh& mesh, std::vector<int> const& selected_material_ids);
+}  // namespace MeshToolsLib
diff --git a/Tests/Data/MeshLib/2D_mesh_with_reset_sigma_ip.vtu b/Tests/Data/MeshLib/2D_mesh_with_reset_sigma_ip.vtu