Merge branch 'reset_ip_data' into 'master'

Reset integration point data by materialIDs

See merge request ogs/ogs!4668

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>
 )
 

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;

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)

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;

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)

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;
  }
 

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.

Applications/Utils/ModelPreparation/PartitionMesh/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

Applications/Utils/ModelPreparation/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

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

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