Move BC meshes post refactoring

Applications/ApplicationsLib/ProjectData.cpp

@@ -21,8 +21,8 @@
 #include "BaseLib/FileTools.h"
 #include "BaseLib/uniqueInsert.h"
 
+#include "GeoLib/GEOObjects.h"
 #include "MathLib/Curve/CreatePiecewiseLinearCurve.h"
-#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
 #include "MeshGeoToolsLib/ConstructMeshesFromGeometries.h"
 #include "MeshGeoToolsLib/CreateSearchLength.h"
 #include "MeshGeoToolsLib/SearchLength.h"
@@ -90,64 +90,77 @@
 #include "ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowWithPrhoProcess.h"
 #endif
 
-namespace detail
+namespace
 {
-static void readGeometry(std::string const& fname,
-                         GeoLib::GEOObjects& geo_objects)
+void readGeometry(std::string const& fname, GeoLib::GEOObjects& geo_objects)
 {
     DBUG("Reading geometry file \'%s\'.", fname.c_str());
     GeoLib::IO::BoostXmlGmlInterface gml_reader(geo_objects);
     gml_reader.readFile(fname);
 }
-}
-
-ProjectData::ProjectData() = default;
 
-ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
-                         std::string const& project_directory,
-                         std::string const& output_directory)
+std::unique_ptr<MeshLib::Mesh> readSingleMesh(
+    BaseLib::ConfigTree const& mesh_config_parameter,
+    std::string const& project_directory)
 {
-    std::string const geometry_file = BaseLib::copyPathToFileName(
-        //! \ogs_file_param{prj__geometry}
-        project_config.getConfigParameter<std::string>("geometry"),
-        project_directory);
-    detail::readGeometry(geometry_file, *_geoObjects);
+    std::string const mesh_file = BaseLib::copyPathToFileName(
+        mesh_config_parameter.getValue<std::string>(), project_directory);
 
+    auto mesh = std::unique_ptr<MeshLib::Mesh>(
+        MeshLib::IO::readMeshFromFile(mesh_file));
+    if (!mesh)
     {
-        //! \ogs_file_param{prj__mesh}
-        auto const mesh_param = project_config.getConfigParameter("mesh");
+        OGS_FATAL("Could not read mesh from \'%s\' file. No mesh added.",
+                  mesh_file.c_str());
+    }
 
-        std::string const mesh_file = BaseLib::copyPathToFileName(
-            mesh_param.getValue<std::string>(), project_directory);
+    if (auto const axially_symmetric =
+            //! \ogs_file_attr{prj__mesh__axially_symmetric}
+        mesh_config_parameter.getConfigAttributeOptional<bool>(
+            "axially_symmetric"))
+    {
+        mesh->setAxiallySymmetric(*axially_symmetric);
+    }
 
-        MeshLib::Mesh* const mesh = MeshLib::IO::readMeshFromFile(mesh_file);
-        if (!mesh)
-        {
-            OGS_FATAL("Could not read mesh from \'%s\' file. No mesh added.",
-                      mesh_file.c_str());
-        }
+    return mesh;
+}
 
-        if (auto const axially_symmetric =
-                //! \ogs_file_attr{prj__mesh__axially_symmetric}
-            mesh_param.getConfigAttributeOptional<bool>("axially_symmetric"))
-        {
-            mesh->setAxiallySymmetric(*axially_symmetric);
-        }
-        _mesh_vec.push_back(mesh);
-    }
+std::vector<std::unique_ptr<MeshLib::Mesh>> readMeshes(
+    BaseLib::ConfigTree const& config, std::string const& project_directory)
+{
+    std::vector<std::unique_ptr<MeshLib::Mesh>> meshes;
+
+    //! \ogs_file_param{prj__mesh}
+    meshes.push_back(
+        readSingleMesh(config.getConfigParameter("mesh"), project_directory));
+
+    std::string const geometry_file = BaseLib::copyPathToFileName(
+        //! \ogs_file_param{prj__geometry}
+        config.getConfigParameter<std::string>("geometry"),
+        project_directory);
+    GeoLib::GEOObjects geoObjects;
+    readGeometry(geometry_file, geoObjects);
 
     std::unique_ptr<MeshGeoToolsLib::SearchLength> search_length_algorithm =
-        MeshGeoToolsLib::createSearchLengthAlgorithm(project_config,
-                                                     *_mesh_vec[0]);
+        MeshGeoToolsLib::createSearchLengthAlgorithm(config, *meshes[0]);
     auto additional_meshes =
         MeshGeoToolsLib::constructAdditionalMeshesFromGeoObjects(
-            *_geoObjects, *_mesh_vec[0], std::move(search_length_algorithm));
+            geoObjects, *meshes[0], std::move(search_length_algorithm));
+
+    std::move(begin(additional_meshes), end(additional_meshes),
+              std::back_inserter(meshes));
+
+    return meshes;
+}
+}  // namespace
+
+ProjectData::ProjectData() = default;
 
-    // release the unique_ptr's while copying to the raw pointers storage.
-    // TODO (naumov) Store unique_ptr's in _mesh_vec.
-    std::transform(begin(additional_meshes), end(additional_meshes),
-                   std::back_inserter(_mesh_vec),
-                   [](auto&& mesh) { return mesh.release(); });
+ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
+                         std::string const& project_directory,
+                         std::string const& output_directory)
+{
+    _mesh_vec = readMeshes(project_config, project_directory);
 
     //! \ogs_file_param{prj__curves}
     parseCurves(project_config.getConfigSubtreeOptional("curves"));
@@ -173,14 +186,6 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
                   output_directory);
 }
 
-ProjectData::~ProjectData()
-{
-    delete _geoObjects;
-
-    for (MeshLib::Mesh* m : _mesh_vec)
-        delete m;
-}
-
 void ProjectData::parseProcessVariables(
     BaseLib::ConfigTree const& process_variables_config)
 {

Applications/ApplicationsLib/ProjectData.h

@@ -18,17 +18,11 @@
 #include <string>
 
 #include "BaseLib/ConfigTree.h"
-
-#include "GeoLib/GEOObjects.h"
-
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
 #include "ProcessLib/Parameter/Parameter.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/ProcessVariable.h"
 
-namespace MathLib
-{
-class PiecewiseLinearInterpolation;
-}
 namespace MeshLib
 {
 class Mesh;
@@ -71,10 +65,6 @@ class ProjectData final
 
     ProjectData(ProjectData&) = delete;
 
-    ~ProjectData();
-
-    /// Returns the GEOObjects containing all points, polylines and surfaces.
-    GeoLib::GEOObjects* getGEOObjects() { return _geoObjects; }
 
     //
     // Process interface
@@ -116,8 +106,7 @@ class ProjectData final
 
     void parseCurves(boost::optional<BaseLib::ConfigTree> const& config);
 
-    GeoLib::GEOObjects* _geoObjects = new GeoLib::GEOObjects();
-    std::vector<MeshLib::Mesh*> _mesh_vec;
+    std::vector<std::unique_ptr<MeshLib::Mesh>> _mesh_vec;
     std::map<std::string, std::unique_ptr<ProcessLib::Process>> _processes;
     std::vector<ProcessLib::ProcessVariable> _process_variables;
 

Applications/Utils/MeshGeoTools/ConstructMeshesFromGeometry.cpp

@@ -37,6 +37,16 @@ int main(int argc, char* argv[])
         ' ',
         "0.1");
 
+    TCLAP::ValueArg<double> search_length_arg(
+        "s",
+        "searchlength",
+        "search length determining radius for the node search algorithm. "
+        "Non-negative floating point number (default 1e-16) ",
+        false,
+        1e-16,
+        "float");
+    cmd.add(search_length_arg);
+
     TCLAP::ValueArg<std::string> geometry_arg("g",
                                               "geometry",
                                               "the file name the geometry",
@@ -60,8 +70,12 @@ int main(int argc, char* argv[])
 
     auto const geo_objects = readGeometry(geometry_arg.getValue());
 
+    double const search_length = search_length_arg.getValue();
+
     auto const extracted_meshes = constructAdditionalMeshesFromGeoObjects(
-        *geo_objects, *mesh, std::make_unique<MeshGeoToolsLib::SearchLength>());
+        *geo_objects,
+        *mesh,
+        std::make_unique<MeshGeoToolsLib::SearchLength>(search_length));
 
     for (auto const& m_ptr : extracted_meshes)
     {

MeshGeoToolsLib/ConstructMeshesFromGeometries.cpp

@@ -106,6 +106,12 @@ constructAdditionalMeshesFromGeoObjects(GeoLib::GEOObjects const& geo_objects,
                   std::back_inserter(additional_meshes));
     }
 
+    // Set axial symmetry for the additional meshes to the same value as the
+    // "bulk" mesh.
+    std::for_each(begin(additional_meshes), end(additional_meshes),
+                  [axial_symmetry = mesh.isAxiallySymmetric()](auto& m) {
+                      m->setAxiallySymmetric(axial_symmetry);
+                  });
     return additional_meshes;
 }
 }  // namespace MeshGeoToolsLib

NumLib/DOF/LocalToGlobalIndexMap.cpp

@@ -189,12 +189,12 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
     std::vector<MeshLib::MeshSubset>&& mesh_subsets,
     std::vector<int> const& global_component_ids,
+    std::vector<int> const& variable_component_offsets,
     std::vector<MeshLib::Element*> const& elements,
     NumLib::MeshComponentMap&& mesh_component_map)
     : _mesh_subsets(std::move(mesh_subsets)),
       _mesh_component_map(std::move(mesh_component_map)),
-      _variable_component_offsets(
-          to_cumulative(std::vector<int>(1, 1)))  // Single variable only.
+      _variable_component_offsets(variable_component_offsets)
 {
     // Each subset in the mesh_subsets represents a single component.
     if (_mesh_subsets.size() != global_component_ids.size())
@@ -247,9 +247,9 @@ LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
         all_mesh_subsets.emplace_back(new_mesh_subset);
     all_mesh_subsets.emplace_back(std::move(new_mesh_subset));
 
-    return new LocalToGlobalIndexMap(std::move(all_mesh_subsets),
-                                     global_component_ids, elements,
-                                     std::move(mesh_component_map));
+    return new LocalToGlobalIndexMap(
+        std::move(all_mesh_subsets), global_component_ids,
+        _variable_component_offsets, elements, std::move(mesh_component_map));
 }
 
 std::size_t LocalToGlobalIndexMap::dofSizeWithGhosts() const

NumLib/DOF/LocalToGlobalIndexMap.h

@@ -153,6 +153,7 @@ class LocalToGlobalIndexMap final
     explicit LocalToGlobalIndexMap(
         std::vector<MeshLib::MeshSubset>&& mesh_subsets,
         std::vector<int> const& global_component_ids,
+        std::vector<int> const& variable_component_offsets,
         std::vector<MeshLib::Element*> const& elements,
         NumLib::MeshComponentMap&& mesh_component_map);
 

NumLib/DOF/MeshComponentMap.cpp

@@ -172,7 +172,21 @@ MeshComponentMap MeshComponentMap::getSubset(
         {
             auto const global_index =
                 getGlobalIndex(bulk_location, component_id);
-            assert (global_index != nop);
+            if (global_index == nop)
+            {
+                OGS_FATAL(
+                    "Could not find a global index for global component %d for "
+                    "the mesh '%s', node %d, in the corresponding bulk mesh "
+                    "'%s' and node %d. This happens because the boundary mesh "
+                    "is larger then the definition region of the bulk "
+                    "component, usually because the geometry for the boundary "
+                    "condition is too large.",
+                    component_id,
+                    new_mesh_subset.getMesh().getName().c_str(),
+                    node_id,
+                    bulk_mesh_subsets.front().getMesh().getName().c_str(),
+                    bulk_node_ids_map[node_id]);
+            }
             subset_dict.insert({new_location, component_id, global_index});
         }
     }