From 47448642026bee0c677756830501293f80d2d261 Mon Sep 17 00:00:00 2001 From: valeriabarra Date: Thu, 16 Jul 2020 20:18:59 -0600 Subject: [PATCH] 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 --- .../shallow-water/qfunctions/setup_geo.h | 2 +- examples/fluids/shallow-water/shallowwater.c | 6 +- examples/fluids/shallow-water/src/setup.c | 83 +++++++++++++++++-- examples/fluids/shallow-water/sw_headers.h | 14 ++++ 4 files changed, 98 insertions(+), 7 deletions(-) diff --git a/examples/fluids/shallow-water/qfunctions/setup_geo.h b/examples/fluids/shallow-water/qfunctions/setup_geo.h index 3949c10e50..b79a85083b 100644 --- a/examples/fluids/shallow-water/qfunctions/setup_geo.h +++ b/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 diff --git a/examples/fluids/shallow-water/shallowwater.c b/examples/fluids/shallow-water/shallowwater.c index 683f044fef..4e769f9fa2 100644 --- a/examples/fluids/shallow-water/shallowwater.c +++ b/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, diff --git a/examples/fluids/shallow-water/src/setup.c b/examples/fluids/shallow-water/src/setup.c index 61dbe6b2c2..8fa6f82db7 100644 --- a/examples/fluids/shallow-water/src/setup.c +++ b/examples/fluids/shallow-water/src/setup.c @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -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, §ion); 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