Allow multiple meshes in project

Applications/ApplicationsLib/ProjectData.cpp

@@ -112,6 +112,7 @@ std::unique_ptr<MeshLib::Mesh> readSingleMesh(
 {
     std::string const mesh_file = BaseLib::copyPathToFileName(
         mesh_config_parameter.getValue<std::string>(), project_directory);
+    DBUG("Reading mesh file \'%s\'.", mesh_file.c_str());
 
     auto mesh = std::unique_ptr<MeshLib::Mesh>(
         MeshLib::IO::readMeshFromFile(mesh_file));
@@ -121,6 +122,11 @@ std::unique_ptr<MeshLib::Mesh> readSingleMesh(
                   mesh_file.c_str());
     }
 
+#ifdef DOXYGEN_DOCU_ONLY
+    //! \ogs_file_attr{prj__meshes__mesh__axially_symmetric}
+    mesh_config_parameter.getConfigAttributeOptional<bool>("axially_symmetric");
+#endif  // DOXYGEN_DOCU_ONLY
+
     if (auto const axially_symmetric =
             //! \ogs_file_attr{prj__mesh__axially_symmetric}
         mesh_config_parameter.getConfigAttributeOptional<bool>(
@@ -137,26 +143,45 @@ std::vector<std::unique_ptr<MeshLib::Mesh>> readMeshes(
 {
     std::vector<std::unique_ptr<MeshLib::Mesh>> meshes;
 
-    meshes.push_back(
+    //! \ogs_file_param{prj__meshes}
+    auto optional_meshes = config.getConfigSubtreeOptional("meshes");
+    if (optional_meshes)
+    {
+        DBUG("Reading multiple meshes.");
+        for (auto mesh_config :
+             //! \ogs_file_param{prj__meshes__mesh}
+             optional_meshes->getConfigParameterList("mesh"))
+        {
+            meshes.push_back(readSingleMesh(mesh_config, project_directory));
+        }
+    }
+    else
+    {  // Read single mesh with geometry.
+        WARN(
+            "Consider switching from mesh and geometry input to multiple "
+            "meshes input. See "
+            "https://www.opengeosys.org/docs/tools/model-preparation/"
+            "constructmeshesfromgeometry/ tool for conversion.");
         //! \ogs_file_param{prj__mesh}
-        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(config, *meshes[0]);
-    auto additional_meshes =
-        MeshGeoToolsLib::constructAdditionalMeshesFromGeoObjects(
-            geoObjects, *meshes[0], std::move(search_length_algorithm));
-
-    std::move(begin(additional_meshes), end(additional_meshes),
-              std::back_inserter(meshes));
-
+        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(config, *meshes[0]);
+        auto additional_meshes =
+            MeshGeoToolsLib::constructAdditionalMeshesFromGeoObjects(
+                geoObjects, *meshes[0], std::move(search_length_algorithm));
+
+        std::move(begin(additional_meshes), end(additional_meshes),
+                  std::back_inserter(meshes));
+    }
     return meshes;
 }
 }  // namespace
@@ -215,20 +240,26 @@ void ProjectData::parseProcessVariables(
 {
     DBUG("Parse process variables:");
 
-    if (_mesh_vec.empty() || _mesh_vec[0] == nullptr)
-    {
-        ERR("A mesh is required to define process variables.");
-        return;
-    }
-
     std::set<std::string> names;
 
     for (auto var_config
          //! \ogs_file_param{prj__process_variables__process_variable}
          : process_variables_config.getConfigSubtreeList("process_variable"))
     {
-        auto pv =
-            ProcessLib::ProcessVariable{var_config, _mesh_vec, _parameters};
+        // Either the mesh name is given, or the first mesh's name will be
+        // taken. Taking the first mesh's value is deprecated.
+        auto const mesh_name =
+            //! \ogs_file_param{prj__process_variables__process_variable__mesh}
+            var_config.getConfigParameter<std::string>("mesh",
+                                                       _mesh_vec[0]->getName());
+
+        auto& mesh = *BaseLib::findElementOrError(
+            begin(_mesh_vec), end(_mesh_vec),
+            [&mesh_name](auto const& m) { return m->getName() == mesh_name; },
+            "Expected to find a mesh named " + mesh_name + ".");
+
+        auto pv = ProcessLib::ProcessVariable{var_config, mesh, _mesh_vec,
+                                              _parameters};
         if (!names.insert(pv.getName()).second)
             OGS_FATAL("A process variable with name `%s' already exists.",
                       pv.getName().c_str());

Applications/FileIO/XmlIO/OpenGeoSysProject.xsd

@@ -15,6 +15,16 @@
     </xs:sequence>
   </xs:complexType>
 
+  <xs:complexType name="meshType">
+    <xs:complexType>
+      <xs:simpleContent>
+        <xs:extension base="xs:string">
+          <xs:attribute name="axially_symmetric" type="xs:boolean"/>
+        </xs:extension>
+      </xs:simpleContent>
+    </xs:complexType>
+  </xs:complexType>
+
   <xs:complexType name="parameterType">
     <xs:sequence>
       <xs:element ref="name"  minOccurs="1" maxOccurs="1" />
@@ -29,8 +39,8 @@
       <xs:element name="geometrical_set" type="xs:string" minOccurs="0" />
       <xs:element name="geometry" type="xs:string" minOccurs="0" />
       <xs:element name="type" type="xs:string" />
-      <xs:element name="field_name" type="xs:string" minOccurs="0" />
       <xs:element name="mesh" type="xs:string" minOccurs="0" />
+      <xs:element name="field_name" type="xs:string" minOccurs="0" />
       <xs:element name="parameter" type="xs:string" minOccurs="0" />
     </xs:sequence>
   </xs:complexType>
@@ -50,6 +60,7 @@
   <xs:complexType name="pvariableType">
     <xs:sequence>
       <xs:element ref="name"/>
+      <xs:element name="mesh" type="xs:string" minOccurs="0"/>
       <xs:element name="components" minOccurs="1" maxOccurs="1">
         <xs:simpleType>
           <xs:restriction base="xs:nonNegativeInteger">
@@ -76,16 +87,15 @@
   <xs:element name="OpenGeoSysProject">
     <xs:complexType>
       <xs:sequence>
-        <xs:element name="mesh" minOccurs="0">
+        <xs:element name="mesh" type="meshType" minOccurs="0"/>
+        <xs:element name="geometry" type="xs:string" minOccurs="0"/>
+        <xs:element name="meshes" minOccurs="0" maxOccurs="1">
           <xs:complexType>
-            <xs:simpleContent>
-              <xs:extension base="xs:string">
-                <xs:attribute name="axially_symmetric" type="xs:boolean"/>
-              </xs:extension>
-            </xs:simpleContent>
+            <xs:sequence>
+              <xs:element name="mesh" type="meshType" minOccurs="0" maxOccurs="unbounded"/>
+            </xs:sequence>
           </xs:complexType>
         </xs:element>
-        <xs:element name="geometry" type="xs:string" minOccurs="0"/>
         <xs:element name="processes" minOccurs="0"/> <!--ignore-->
         <xs:element name="time_loop" minOccurs="0"/> <!--ignore-->
         <xs:element name="parameters" minOccurs="0">

Applications/FileIO/XmlIO/Qt/XmlPrjInterface.cpp

@@ -63,6 +63,15 @@ int XmlPrjInterface::readFile(const QString& fileName)
         return 0;
     }
 
+    auto read_single_mesh = [&](QString const& mesh_str) {
+        std::unique_ptr<MeshLib::Mesh> mesh{
+            MeshLib::IO::readMeshFromFile(mesh_str.toStdString())};
+        if (mesh != nullptr)
+        {
+            _project.addMesh(std::move(mesh));
+        }
+    };
+
     QDomNodeList fileList = docElement.childNodes();
 
     for (int i = 0; i < fileList.count(); i++)
@@ -93,11 +102,41 @@ int XmlPrjInterface::readFile(const QString& fileName)
         }
         else if (node_name == "mesh")
         {
-            QString const mesh_str(path + file_name);
-            std::unique_ptr<MeshLib::Mesh> mesh(
-                MeshLib::IO::readMeshFromFile(mesh_str.toStdString()));
-            if (mesh)
-                _project.addMesh(std::move(mesh));
+            read_single_mesh(path + file_name);
+        }
+        else if (node_name == "meshes")
+        {
+            for (QDomNode meshes_node = node.firstChild();
+                 meshes_node != QDomNode();
+                 meshes_node = meshes_node.nextSibling())
+            {
+                if (!meshes_node.isElement())
+                {
+                    ERR("Expected an XML element node.")
+                    return 0;
+                }
+                if (meshes_node.nodeName() != "mesh")
+                {
+                    ERR("Expected an XML element node named 'mesh' got '%s'.",
+                        meshes_node.nodeName().data())
+                    return 0;
+                }
+                if (meshes_node.childNodes().count() != 1)
+                {
+                    ERR("Expected an XML element node named 'mesh' to contain "
+                        "exactly one child node but it has %d children.",
+                        meshes_node.childNodes().count());
+                    return 0;
+                }
+                QDomNode node_text = meshes_node.firstChild();
+                if (!node_text.isText())
+                {
+                    ERR("Expected an XML element node named 'mesh' to contain "
+                        "text.")
+                    return 0;
+                }
+                read_single_mesh(path + node_text.toText().data().trimmed());
+            }
         }
 
         else if (node_name == "parameters")

Documentation/Doxyfile.in

@@ -2023,7 +2023,7 @@ INCLUDE_FILE_PATTERNS  =
 # recursively expanded use the := operator instead of the = operator.
 # This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
 
-PREDEFINED             =
+PREDEFINED             = DOXYGEN_DOCU_ONLY
 
 # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this
 # tag can be used to specify a list of macro names that should be expanded. The

Documentation/ProjectFile/prj/mesh/a_axially_symmetric.md

@@ -1 +1,4 @@
-\ogs_missing_documentation
+Axial symmetry can be used to reduce complexity of a problem by reducing the
+geometrical dimension by one. The symmetry axis is y-axis in 1D and 2D cases.
+
+Only applicable to 1D and 2D meshes.

Documentation/ProjectFile/prj/meshes/i_meshes.md

@@ -0,0 +1 @@
+External mesh file names for the process variables and the boundary conditions.

Documentation/ProjectFile/prj/meshes/mesh/a_axially_symmetric.md

@@ -0,0 +1,4 @@
+Axial symmetry can be used to reduce complexity of a problem by reducing the
+geometrical dimension by one. The symmetry axis is y-axis in 1D and 2D cases.
+
+Only applicable to 1D and 2D meshes.

Documentation/ProjectFile/prj/meshes/mesh/i_mesh.md

@@ -0,0 +1 @@
+Information about input mesh.

Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/t_mesh.md

@@ -0,0 +1,10 @@
+Name of the surface mesh where the boundary condition will be defined.
+
+The surface mesh must contain a nodal integer-valued field (unsigned 64bit)
+named \c bulk_node_ids, and a cell field named \c bulk_element_ids. These fields
+establish the mapping between the nodes of the surface mesh to the nodes in the
+bulk mesh.
+
+\warning It is not checked if the surface mesh and the bulk mesh
+correspond to each other; in particular it is not checked if surface and bulk
+nodes coincide and if surface elements coincide with the faces of bulk elements.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_mesh.md

@@ -0,0 +1,10 @@
+Name of the mesh where the source term will be applied.
+
+The surface mesh must contain a nodal integer-valued field (unsigned 64bit)
+named \c bulk_node_ids, and a cell field named \c bulk_element_ids. These fields
+establish the mapping between the nodes of the surface mesh to the nodes in the
+bulk mesh.
+
+\warning It is not checked if the surface mesh and the bulk mesh
+correspond to each other; in particular it is not checked if surface and bulk
+nodes coincide and if surface elements coincide with the faces of bulk elements.

Documentation/ProjectFile/prj/process_variables/process_variable/t_mesh.md

@@ -0,0 +1 @@
+A mesh name for the process variable's domain of definition.

MeshGeoToolsLib/ConstructMeshesFromGeometries.cpp

@@ -20,6 +20,12 @@
 
 namespace MeshGeoToolsLib
 {
+std::string meshNameFromGeometry(std::string const& geometrical_set_name,
+                                 std::string const& geometry_name)
+{
+    return geometrical_set_name + "_" + geometry_name;
+}
+
 template <typename GeometryVec>
 std::vector<std::unique_ptr<MeshLib::Mesh>>
 constructAdditionalMeshesFromGeometries(
@@ -55,7 +61,7 @@ constructAdditionalMeshesFromGeometries(
                  geometry_name.c_str());
 
             additional_meshes.emplace_back(createMeshFromElementSelection(
-                vec_name + "_" + geometry_name,
+                meshNameFromGeometry(vec_name, geometry_name),
                 MeshLib::cloneElements(
                     boundary_element_searcher.getBoundaryElements(geometry))));
         }

MeshGeoToolsLib/ConstructMeshesFromGeometries.h

@@ -31,4 +31,7 @@ constructAdditionalMeshesFromGeoObjects(GeoLib::GEOObjects const& geo_objects,
                                         MeshLib::Mesh const& mesh,
                                         std::unique_ptr<SearchLength>
                                             search_length_algorithm);
+
+std::string meshNameFromGeometry(std::string const& geometrical_set_name,
+                                 std::string const& geometry_name);
 }  // namespace MeshGeoToolsLib

MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp

@@ -563,13 +563,9 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::newMesh(
     MeshLib::Properties const& properties) const
 {
     return new MeshLib::NodePartitionedMesh(
-        mesh_name + std::to_string(_mpi_comm_size),
-        mesh_nodes, glb_node_ids, mesh_elems,
-        properties,
-        _mesh_info.global_base_nodes,
-        _mesh_info.global_nodes,
-        _mesh_info.base_nodes,
-        _mesh_info.active_base_nodes,
+        mesh_name, mesh_nodes, glb_node_ids, mesh_elems, properties,
+        _mesh_info.global_base_nodes, _mesh_info.global_nodes,
+        _mesh_info.base_nodes, _mesh_info.active_base_nodes,
         _mesh_info.active_nodes);
 }
 

MeshLib/Node.cpp

@@ -49,5 +49,18 @@ void Node::updateCoordinates(double x, double y, double z)
         _elements[i]->computeVolume();
 }
 
+bool isBaseNode(Node const& node)
+{
+    // Check if node is connected.
+    if (node.getNumberOfElements() == 0)
+        return true;
+
+    // In a mesh a node always belongs to at least one element.
+    auto const e = node.getElement(0);
+
+    auto const n_base_nodes = e->getNumberOfBaseNodes();
+    auto const local_index = e->getNodeIDinElement(&node);
+    return local_index < n_base_nodes;
+}
 }
 

MeshLib/Node.h

@@ -102,4 +102,8 @@ class Node final : public MathLib::Point3dWithID
     std::vector<Element*> _elements;
 }; /* class */
 
+/// Returns true if the given node is a base node of a (first) element, or if it
+/// is not connected to any element i.e. an unconnected node.
+bool isBaseNode(Node const& node);
+
 } /* namespace */

NumLib/DOF/LocalToGlobalIndexMap.h

@@ -96,7 +96,7 @@ class LocalToGlobalIndexMap final
     LocalToGlobalIndexMap* deriveBoundaryConstrainedMap(
         int const variable_id,
         std::vector<int> const& component_ids,
-        MeshLib::MeshSubset&& mesh_subset) const;
+        MeshLib::MeshSubset&& new_mesh_subset) const;
 
     /// Derive a LocalToGlobalIndexMap constrained to the mesh subset and mesh
     /// subset's elements. A new mesh component map will be constructed using