SWE: Add FindPanelEdgeNodes() function to determine which nodes are on the cube edges

The algorithm with bit manipulation was greatly inspired-by and helped-by @matteovigni

Thanks-to @matteovigni

Author: valeriabarra

examples/fluids/shallow-water/qfunctions/setup_geo.h

@@ -171,7 +171,7 @@ CEED_QFUNCTION(SetupGeo)(void *ctx, CeedInt Q,
     qdata[7][i] = dxdXTdxdXinv[0][0];
     qdata[8][i] = dxdXTdxdXinv[1][1];
     qdata[9][i] = dxdXTdxdXinv[0][1];
-    
+
     // Terrain topography, hs
     qdata[10][i] = sin(xx[0]) + cos(xx[1]); // put 0 for constant flat topography
   } // End of Quadrature Point Loop

examples/fluids/shallow-water/shallowwater.c

@@ -56,6 +56,7 @@ int main(int argc, char **argv) {
   Units units;
   problemType problemChoice;
   problemData *problem = NULL;
+  EdgeNode edgenodes;
   StabilizationType stab;
   PetscBool implicit;
   PetscInt degree, qextra, outputfreq, steps, contsteps;
@@ -65,7 +66,7 @@ int main(int argc, char **argv) {
   const CeedInt ncompx = 3;
   PetscInt viz_refine = 0;
   PetscBool read_mesh, simplex, test;
-  PetscInt topodim = 2, ncompq = 3, lnodes;
+  PetscInt topodim = 2, ncompq = 3, lnodes, nedgenodes;
   // libCEED context
   char ceedresource[PETSC_MAX_PATH_LEN] = "/cpu/self",
                                           filename[PETSC_MAX_PATH_LEN];
@@ -374,6 +375,9 @@ int main(int argc, char **argv) {
   // Set up global mass vector
   ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr);
 
+  // Setup global lat-long vector for different panels (charts) of the cube
+  ierr = FindPanelEdgeNodes(dm, ncompq, &nedgenodes, &edgenodes); CHKERRQ(ierr);
+
   // Setup libCEED's objects
   ierr = PetscMalloc1(1, &ceeddata); CHKERRQ(ierr);
   ierr = SetupLibceed(dm, ceed, degree, qextra, ncompx, ncompq, user, ceeddata,

examples/fluids/shallow-water/src/setup.c

@@ -20,6 +20,7 @@
 #include <stdbool.h>
 #include <string.h>
 #include <petsc.h>
+#include <petscsys.h>
 #include <petscdmplex.h>
 #include <petscfe.h>
 #include <ceed.h>
@@ -66,6 +67,78 @@ problemData problemOptions[] = {
   }
 };
 
+// -----------------------------------------------------------------------------
+// Auxiliary function to determine if nodes belong to cube faces (panels)
+// -----------------------------------------------------------------------------
+
+PetscErrorCode FindPanelEdgeNodes(DM dm, PetscInt ncomp, PetscInt *nedgenodes,
+                                  EdgeNode *edgenodes) {
+
+  PetscInt ierr;
+  PetscInt cstart, cend, nstart, nend, nnodes, depth, edgenode = 0;
+  PetscSection section;
+  PetscFunctionBeginUser;
+  // Get Nelem
+  ierr = DMGetSection(dm, &section); CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 0, &cstart, &cend); CHKERRQ(ierr);
+  ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, depth, &nstart, &nend); CHKERRQ(ierr);
+  nnodes = nend - nstart;
+  unsigned int bitmap[nnodes];
+  ierr = PetscMemzero(bitmap, sizeof(bitmap)); CHKERRQ(ierr);
+
+  // Get indices
+  for (PetscInt c = cstart; c < cend; c++) { // Traverse elements
+    PetscInt numindices, *indices, n, panel;
+    // Query element panel
+    ierr = DMGetLabelValue(dm, "panel", c, &panel); CHKERRQ(ierr);
+    ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE,
+                                   &numindices, &indices, NULL, NULL);
+    CHKERRQ(ierr);
+    for (n = 0; n < numindices; n += ncomp) { // Traverse nodes per element
+      PetscInt  bitmapidx = indices[n] / ncomp;
+      bitmap[bitmapidx] |= (1 << panel);
+    }
+    ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE,
+                                       &numindices, &indices, NULL, NULL);
+    CHKERRQ(ierr);
+  }
+  // Read the 1's in the resulting bitmap and extract edge nodes only
+  ierr = PetscMalloc1(nnodes + 24, edgenodes); CHKERRQ(ierr);
+  for (PetscInt i = 0; i < nnodes; i++) {
+    PetscInt ones = 0, panels[3];
+    for (PetscInt p = 0; p < 6; p++) {
+      if (bitmap[i] & 1)
+        panels[ones++] = p;
+      bitmap[i] >>= 1;
+    }
+    if (ones == 2) {
+      (*edgenodes)[edgenode].idx = i;
+      (*edgenodes)[edgenode].panelA = panels[0];
+      (*edgenodes)[edgenode].panelB = panels[1];
+      edgenode++;
+    }
+    else if (ones == 3) {
+      (*edgenodes)[edgenode].idx = i;
+      (*edgenodes)[edgenode].panelA = panels[0];
+      (*edgenodes)[edgenode].panelB = panels[1];
+      edgenode++;
+      (*edgenodes)[edgenode].idx = i;
+      (*edgenodes)[edgenode].panelA = panels[0];
+      (*edgenodes)[edgenode].panelB = panels[2];
+      edgenode++;
+      (*edgenodes)[edgenode].idx = i;
+      (*edgenodes)[edgenode].panelA = panels[1];
+      (*edgenodes)[edgenode].panelB = panels[2];
+      edgenode++;
+    }
+  }
+  ierr = PetscRealloc(edgenode, edgenodes); CHKERRQ(ierr);
+  *nedgenodes = edgenode;
+
+  PetscFunctionReturn(0);
+}
+
 // -----------------------------------------------------------------------------
 // Auxiliary function to create PETSc FE space for a given degree
 // -----------------------------------------------------------------------------
@@ -178,7 +251,7 @@ PetscErrorCode SetupDM(DM dm, PetscInt degree, PetscInt ncompq,
 // PETSc sphere auxiliary function
 // -----------------------------------------------------------------------------
 
-// Utility function taken from petsc/src/dm/impls/plex/examples/tutorials/ex7.c
+// Utility function taken from petsc/src/dm/impls/plex/tutorials/ex7.c
 PetscErrorCode ProjectToUnitSphere(DM dm) {
   Vec            coordinates;
   PetscScalar   *coords;
@@ -223,13 +296,13 @@ PetscErrorCode CreateRestrictionPlex(Ceed ceed, DM dm, CeedInt P,
 
   // Get indices
   ierr = PetscMalloc1(nelem*P*P, &erestrict); CHKERRQ(ierr);
-  for (c=cStart, eoffset = 0; c<cEnd; c++) {
+  for (c = cStart, eoffset = 0; c < cEnd; c++) {
     PetscInt numindices, *indices, i;
     ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE,
                                    &numindices, &indices, NULL, NULL);
     CHKERRQ(ierr);
-    for (i=0; i<numindices; i+=ncomp) {
-      for (PetscInt j=0; j<ncomp; j++) {
+    for (i = 0; i < numindices; i += ncomp) {
+      for (PetscInt j = 0; j < ncomp; j++) {
         if (indices[i+j] != indices[i] + (PetscInt)(copysign(j, indices[i])))
           SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
                    "Cell %D closure indices not interlaced", c);
@@ -244,7 +317,7 @@ PetscErrorCode CreateRestrictionPlex(Ceed ceed, DM dm, CeedInt P,
   }
   if (eoffset != nelem*P*P) SETERRQ3(PETSC_COMM_SELF,
         PETSC_ERR_LIB, "ElemRestriction of size (%D,%D) initialized %D nodes",
-        nelem, P*P,eoffset);
+        nelem, P*P, eoffset);
 
   // Setup CEED restriction
   ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr);

examples/fluids/shallow-water/sw_headers.h

@@ -93,6 +93,7 @@ static const char *const StabilizationTypes[] = {
 // PETSc user data
 typedef struct User_ *User;
 typedef struct Units_ *Units;
+typedef struct EdgeNode_ *EdgeNode;
 
 struct User_ {
   MPI_Comm comm;
@@ -117,6 +118,11 @@ struct Units_ {
   PetscScalar mpersquareds;
 };
 
+struct EdgeNode_ {
+  PetscInt idx; // Node index
+  PetscInt panelA, panelB; // Indices of panels sharing the edge node
+};
+
 // libCEED data struct
 typedef struct CeedData_ *CeedData;
 struct CeedData_ {
@@ -132,6 +138,14 @@ struct CeedData_ {
 // External variables
 extern problemData problemOptions[];
 
+// -----------------------------------------------------------------------------
+// Auxiliary functions for cube face (panel) charts
+// -----------------------------------------------------------------------------
+
+// Auxiliary function to determine if nodes belong to cube faces (panels)
+PetscErrorCode FindPanelEdgeNodes(DM dm, PetscInt ncomp, PetscInt *nedgenodes,
+                                  EdgeNode *edgenodes);
+
 // -----------------------------------------------------------------------------
 // Setup DM functions
 // -----------------------------------------------------------------------------