patchkernel: add a function to check if factes of surface patch have …

…a consistent orientation
optimad · Oct 12, 2023 · 647114d · 647114d
1 parent 5840fde
commit 647114d
Show file tree

Hide file tree

Showing 3 changed files with 117 additions and 0 deletions.
diff --git a/src/patchkernel/surface_kernel.cpp b/src/patchkernel/surface_kernel.cpp
@@ -886,6 +886,113 @@ void SurfaceKernel::evalVertexNormals(long id, int vertex, std::size_t nVertexNe
     }
 }
 
+/*!
+ * Check if the factes have a consistent orientation.
+ *
+ * \return Returns true if the factes have a consistent orientation, false otherwise.
+ */
+bool SurfaceKernel::isCellOrientationConsistent() const
+{
+    CellConstIterator internalBegin = internalCellConstBegin();
+    CellConstIterator internalEnd   = internalCellConstEnd();
+
+    std::size_t nUnprocessed = getInternalCellCount();
+    PiercedStorage<bool> alreadyProcessed(1, &(getCells()));
+    alreadyProcessed.fill(false);
+
+    bool consistent = true;
+    std::vector<std::size_t> processRawList;
+    for (CellConstIterator seedItr = internalBegin; seedItr != internalEnd; ++seedItr) {
+        // Get seed information
+        long seedRawId = seedItr.getRawIndex();
+        if (alreadyProcessed.rawAt(seedRawId)) {
+            continue;
+        }
+
+        // Compare seed orientation with the orientation of its neighbours
+        processRawList.assign(1, seedRawId);
+        while (!processRawList.empty()) {
+            // Get cell to process
+            std::size_t cellRawId = processRawList.back();
+            processRawList.pop_back();
+
+            // Skip cells already processed
+            if (alreadyProcessed.rawAt(cellRawId)) {
+                continue;
+            }
+
+            // Process cell
+            CellConstIterator cellItr = getCells().rawFind(cellRawId);
+            const Cell &cell = *cellItr;
+            int nCellFaces = cell.getFaceCount();
+            for (int face = 0; face < nCellFaces; ++face) {
+                int nFaceAdjacencies = cell.getAdjacencyCount(face);
+                const long *faceAdjacencies = cell.getAdjacencies(face);
+                for (int i = 0; i < nFaceAdjacencies; ++i) {
+                    // Neighbout information
+                    long neighId = faceAdjacencies[i];
+                    CellConstIterator neighItr = getCells().find(neighId);
+                    const Cell &neigh = getCell(neighId);
+
+                    // Check neighbour consistency
+                    int neighFace = findAdjoinNeighFace(cell, face, neigh);
+                    bool isNeighConsistent = haveSameOrientation(cell, face, neigh, neighFace);
+                    if (!isNeighConsistent) {
+                        consistent = false;
+                        break;
+                    }
+
+                    // Add neighbour to the process list
+                    //
+                    // Only internal cells needs to be processed.
+                    if (neigh.isInterior()) {
+                        std::size_t neighRawId = neighItr.getRawIndex();
+                        if (!alreadyProcessed.rawAt(neighRawId)) {
+                            processRawList.push_back(neighRawId);
+                        }
+                    }
+                }
+
+                // Stop processing the cell is the orientation is not consistent
+                if (!consistent) {
+                    break;
+                }
+            }
+
+            // Cell processing has been completed
+            alreadyProcessed.rawSet(cellRawId, true);
+            --nUnprocessed;
+            if (nUnprocessed == 0) {
+                break;
+            }
+
+            // Stop processing the patch is the orientation is not consistent
+            if (!consistent) {
+                break;
+            }
+        }
+
+        // Stop processing the patch is the orientation is not consistent
+        if (!consistent) {
+            break;
+        }
+
+        // Check if all cells have been processed
+        if (nUnprocessed == 0) {
+            break;
+        }
+    }
+
+#if BITPIT_ENABLE_MPI==1
+    // Comunicate consistency flag among the partitions
+    if (isPartitioned()) {
+        MPI_Allreduce(MPI_IN_PLACE, &consistent, 1, MPI_C_BOOL, MPI_LAND, getCommunicator());
+    }
+#endif
+
+    return consistent;
+}
+
 /*!
  * Adjust the orientation of all facets of the local partition according to the
  * orientation of the first facet stored.

diff --git a/src/patchkernel/surface_kernel.hpp b/src/patchkernel/surface_kernel.hpp
@@ -69,6 +69,7 @@ class SurfaceKernel : public PatchKernel {
                                    std::array<double, 3> *unlimitedNormal, std::array<double, 3> *limitedNormal) const;
     double evalCellSize(long id) const override;
 
+    bool isCellOrientationConsistent() const;
     bool adjustCellOrientation();
     bool adjustCellOrientation(long id, bool invert = false);
     void flipCellOrientation(long id);

diff --git a/test/integration_tests/surfunstructured/test_surfunstructured_parallel_00003.cpp b/test/integration_tests/surfunstructured/test_surfunstructured_parallel_00003.cpp
@@ -212,6 +212,15 @@ if (myRank == 0) {
         log::cout() << "   Error during surface orientation" << endl;
     }
 
+    // Check if orientation is consistent ----------------------------------- //
+    log::cout() << "** Mesh adjust orientation" << endl;
+    bool oriented = mesh.isCellOrientationConsistent();
+    if (oriented) {
+        log::cout() << "   Mesh orientation is consistent" << endl;
+    } else {
+        log::cout() << "   Error while checking if mesh orientation is consisten" << endl;
+    }
+
     // Write mesh ----------------------------------------------------------- //
     log::cout() << "** Writing mesh" << endl;
     mesh.write();