Merge pull request #522 from geodynamics/baagaard/fix-meshiopetsc-labels

Fix boundary labels in MeshIOPetsc

Author: baagaard-usgs

libsrc/pylith/meshio/MeshIOPetsc.cc

@@ -29,10 +29,8 @@
 #include "pylith/utils/error.hh" // USES PYLITH_METHOD_*
 #include "pylith/utils/journals.hh" // USES PYLITH_COMPONENT_*
 
-#include <iomanip> // USES setw(), setiosflags(), resetiosflags()
-#include <strings.h> // USES strcasecmp()
+#include <set> // USES std::set
 #include <cassert> // USES assert()
-#include <fstream> // USES std::ifstream, std::ofstream
 #include <stdexcept> // USES std::runtime_error
 #include <sstream> // USES std::ostringstream
 #include <typeinfo> // USES std::typeid
@@ -50,6 +48,17 @@ namespace pylith {
             static
             void fixMaterialLabel(PetscDM* dmMesh);
 
+            /** Create labels for boundaries based on vertices.
+             *
+             * For each boundary condition label:
+             * 1. Remove all points that are not vertices from label.
+             * 2. Add edges and faces to label based upon vertices.
+             *
+             * @param[inout] dmMesh PETSc DM for mesh.
+             */
+            static
+            void fixBoundaryLabels(PetscDM* dmMesh);
+
         }; // _MeshIOPetsc
     } // meshio
 } // pylith
@@ -113,6 +122,7 @@ pylith::meshio::MeshIOPetsc::_read(void) {
     err = DMPlexDistributeSetDefault(dmMesh, PETSC_FALSE);PYLITH_CHECK_ERROR(err);
     err = DMSetFromOptions(dmMesh);PYLITH_CHECK_ERROR(err);
     _MeshIOPetsc::fixMaterialLabel(&dmMesh);
+    _MeshIOPetsc::fixBoundaryLabels(&dmMesh);
     _mesh->setDM(dmMesh);
 
     PYLITH_METHOD_END;
@@ -126,7 +136,7 @@ pylith::meshio::MeshIOPetsc::_write(void) const {  }
 
 
 // ------------------------------------------------------------------------------------------------
-// Fix material label;
+// Remove everything but cells from label for materials.
 void
 pylith::meshio::_MeshIOPetsc::fixMaterialLabel(PetscDM* dmMesh) {
     PYLITH_METHOD_BEGIN;
@@ -161,3 +171,94 @@ pylith::meshio::_MeshIOPetsc::fixMaterialLabel(PetscDM* dmMesh) {
 
     PYLITH_METHOD_END;
 } // fixMaterialLabel
+
+
+// ------------------------------------------------------------------------------------------------
+// Fix boundary condition labels. Remove points other than vertices from all labels other
+// than the material label.
+void
+pylith::meshio::_MeshIOPetsc::fixBoundaryLabels(PetscDM* dmMesh) {
+    PYLITH_METHOD_BEGIN;
+    assert(dmMesh);
+    PetscErrorCode err = 0;
+
+    // Create set with labels to ignore.
+    std::set<std::string> labelsIgnore;
+    labelsIgnore.insert(std::string(pylith::topology::Mesh::cells_label_name));
+    labelsIgnore.insert("celltype");
+    labelsIgnore.insert("depth");
+
+    const PetscInt vertexDepth = 0;
+    PetscInt vStart = -1, vEnd = -1;
+    err = DMPlexGetDepthStratum(*dmMesh, vertexDepth, &vStart, &vEnd);PYLITH_CHECK_ERROR(err);
+
+    PetscInt numLabels = 0;
+    err = DMGetNumLabels(*dmMesh, &numLabels);PYLITH_CHECK_ERROR(err);
+    for (PetscInt iLabel = 0; iLabel < numLabels; ++iLabel) {
+        const char* labelName = NULL;
+        err = DMGetLabelName(*dmMesh, iLabel, &labelName);PYLITH_CHECK_ERROR(err);
+        if (labelsIgnore.count(std::string(labelName)) > 0) { continue; }
+
+        PetscDMLabel dmLabel = NULL;
+        err = DMGetLabelByNum(*dmMesh, iLabel, &dmLabel);PYLITH_CHECK_ERROR(err);
+        PetscInt pStart = -1, pEnd = -1;
+        err = DMLabelGetBounds(dmLabel, &pStart, &pEnd);PYLITH_CHECK_ERROR(err);
+
+        for (PetscInt point = pStart; point < pEnd; ++point) {
+            if ((point >= vStart) && (point < vEnd)) { continue; }
+            PetscBool hasLabel = PETSC_FALSE;
+            err = DMLabelHasPoint(dmLabel, point, &hasLabel);PYLITH_CHECK_ERROR(err);
+            if (hasLabel) {
+                PetscInt labelValue;
+                err = DMLabelGetValue(dmLabel, point, &labelValue);PYLITH_CHECK_ERROR(err);
+                err = DMLabelClearValue(dmLabel, point, labelValue);PYLITH_CHECK_ERROR(err);
+            } // if
+        } // for
+        err = DMLabelDestroyIndex(dmLabel);PYLITH_CHECK_ERROR(err);
+
+        // Readd edges and faces based upon vertices.
+        // Add any non-cells which have all vertices in label.
+        PetscInt cStart = -1, cEnd = -1;
+        err = DMPlexGetHeightStratum(*dmMesh, 0, &cStart, &cEnd);PYLITH_CHECK_ERROR(err);
+
+        err = DMLabelComputeIndex(dmLabel);PYLITH_CHECK_ERROR(err);
+        err = DMLabelGetBounds(dmLabel, &pStart, &pEnd);PYLITH_CHECK_ERROR(err);
+        err = DMLabelDestroyIndex(dmLabel);PYLITH_CHECK_ERROR(err);
+        for (PetscInt vertex = pStart; vertex < pEnd; ++vertex) {
+            PetscInt labelValue = -1;
+            err = DMLabelGetValue(dmLabel, vertex, &labelValue);PYLITH_CHECK_ERROR(err);
+            if (labelValue < 1) { continue; }
+
+            PetscInt *star = NULL, starSize;
+            err = DMPlexGetTransitiveClosure(*dmMesh, vertex, PETSC_FALSE, &starSize, &star);PYLITH_CHECK_ERROR(err);
+            for (PetscInt s = 0; s < starSize*2; s += 2) {
+                const PetscInt point = star[s];
+
+                if ((point >= cStart) && (point < cEnd)) { continue;}
+
+                // All vertices in closure must be in label to add point to label.
+                PetscInt *closure = NULL, closureSize, value;
+                PetscBool mark = PETSC_TRUE;
+                err = DMPlexGetTransitiveClosure(*dmMesh, point, PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
+                for (PetscInt c = 0; c < closureSize*2; c += 2) {
+                    if ((closure[c] >= vStart) && (closure[c] < vEnd)) {
+                        err = DMLabelGetValue(dmLabel, closure[c], &value);PYLITH_CHECK_ERROR(err);
+                        if (value != labelValue) {
+                            mark = PETSC_FALSE;
+                            break;
+                        } // if
+                    } // if
+                } // for
+                err = DMPlexRestoreTransitiveClosure(*dmMesh, point, PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
+                if (mark) {
+                    err = DMLabelSetValue(dmLabel, point, labelValue);PYLITH_CHECK_ERROR(err);
+                } // if
+            } // for
+            err = DMPlexRestoreTransitiveClosure(*dmMesh, vertex, PETSC_FALSE, &starSize, &star);PYLITH_CHECK_ERROR(err);
+        } // for
+        err = DMLabelComputeIndex(dmLabel);PYLITH_CHECK_ERROR(err);
+
+    } // for
+
+    PYLITH_METHOD_END;
+} // fixBoundaryLabels

pylith/meshio/gmsh_utils.py

@@ -189,4 +189,5 @@ def _add_arguments(parser):
         pass
 
 
-# End of file
+# End of file
+

tests/mmstests/linearelasticity/data/generate_gmsh_2d.py

@@ -22,6 +22,7 @@ def __init__(self):
             "required": True,
             "choices": ["tri", "quad"],
             }
+        self.filename = "tri.msh"
 
     def create_geometry(self):
         """Create geometry.

tests/mmstests/linearelasticity/data/generate_gmsh_3d.py

@@ -8,13 +8,14 @@ class App(GenerateMesh):
     Block is DOMAIN_X by DOMAIN_Y by DOMAIN_Z with discretization size DX.
     """
     DOMAIN_X = DOMAIN_Y = DOMAIN_Z = 8.0e+3
-    DX = 4.0e+3
+    DX = 2.0e+3
 
     def __init__(self):
         self.cell_choices = {
             "required": True,
             "choices": ["tet", "hex"],
             }
+        self.filename = "tet.msh"
 
     def create_geometry(self):
         """Create geometry.