Hello PETSc team,
I am trying to extrude some cells from a mesh that comes from a
DMForest. The DMForest is converted to a DMPlex which is then being used
in a modified version of DMPlexExtrude (patch is included). The issue is
the DMForest to DMPlex conversion switches the stratum ordering. It
appears that DMPlexTransform assumes that you have the standard DMPlex
ordering: Cells->Vertices->Edges in 2D. The plex that comes from
DMForest has Cells->Edges->Vertices. This results in an edge being
incorrectly flagged as outside of the range of edges.
I'm including a minimal code showing this. Sample output and the error
message is also shown below.
Thanks,
Dave
Original Plex:
0: [100,221)
1: [221,441)
2: [0,100)
Before Extrude Box: +0.000000e+00, +1.000000e+00 +0.000000e+00,
+1.000000e+00
After Extrude Box: -1.000000e-01, +1.100000e+00 -1.000000e-01,
+1.100000e+00
Converted Plex:
0: [320,441)
1: [100,320)
2: [0,100)
Before Extrude Box: +0.000000e+00, +1.000000e+00 +0.000000e+00,
+1.000000e+00
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: PETSc has generated inconsistent data
[0]PETSC ERROR: Point 100 is not a segment [221, 441)
[0]PETSC ERROR: See
https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!cKv8WrCLXVqU9jE_pHlDNyVtTMk8aT2xxJcPTjksuQP4PiWtA5Kv8E0Ph8jKhOm7MbkM6qFHmK3R5N4y_sQ$
for trouble shooting.
[0]PETSC ERROR: PETSc Development Git Revision: v3.23.3-199-g3b15527e756
Git Date: 2025-06-03 15:20:49 +0000
[0]PETSC ERROR: ./ex with 1 MPI process(es) and PETSC_ARCH
arch-ablate-opt on libuse.eng.buffalo.edu by davidsal Wed Jun 4
10:39:41 2025
[0]PETSC ERROR: Configure options: --with-64-bit-indices=0
--download-metis --download-mpich --download-p4est --download-parmetis
--download-scalapack --download-slepc --download-tetgen
--download-triangle --download-zlib --with-debugging=0 --with-libpng
--with-slepc PETSC_ARCH=arch-ablate-opt --with-fortran-bindings=0
[0]PETSC ERROR: #1 DMPlexTransformGetTargetPoint() at
/home/davidsal/Research/petsc/src/dm/impls/plex/transform/interface/plextransform.c:1003
[0]PETSC ERROR: #2 DMPlexTransformSetConeSizes() at
/home/davidsal/Research/petsc/src/dm/impls/plex/transform/interface/plextransform.c:1400
[0]PETSC ERROR: #3 DMPlexTransformApply() at
/home/davidsal/Research/petsc/src/dm/impls/plex/transform/interface/plextransform.c:2385
[0]PETSC ERROR: #4 DMPlexExtrude() at
/home/davidsal/Research/petsc/src/dm/impls/plex/plexextrude.c:82
[0]PETSC ERROR: #5 extrudeMesh() at ex.c:34
[0]PETSC ERROR: #6 main() at ex.c:86
[0]PETSC ERROR: No PETSc Option Table entries
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to petsc-ma...@mcs.anl.gov----------
diff --git a/include/petscdmplex.h b/include/petscdmplex.h
index 6b52d8b..8ba1f7d 100644
--- a/include/petscdmplex.h
+++ b/include/petscdmplex.h
@@ -187,7 +187,7 @@ PETSC_EXTERN PetscErrorCode DMPlexCreateTPSMesh(MPI_Comm, DMPlexTPSType, const P
PETSC_EXTERN PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm, PetscInt, PetscBool, DM *);
PETSC_EXTERN PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm, const PetscInt[], const PetscReal[], const PetscReal[], const DMBoundaryType[], PetscBool, PetscBool, DM *);
PETSC_EXTERN PetscErrorCode DMPlexCreateHypercubicMesh(MPI_Comm, PetscInt, const PetscInt[], const PetscReal[], const PetscReal[], PetscInt, DM *);
-PETSC_EXTERN PetscErrorCode DMPlexExtrude(DM, PetscInt, PetscReal, PetscBool, PetscBool, PetscBool, const PetscReal[], const PetscReal[], DM *);
+PETSC_EXTERN PetscErrorCode DMPlexExtrude(DM, PetscInt, PetscReal, PetscBool, PetscBool, PetscBool, const PetscReal[], const PetscReal[], DMLabel, DM *);
PETSC_EXTERN PetscErrorCode DMPlexSetIsoperiodicFaceSF(DM, PetscInt, PetscSF *);
PETSC_EXTERN PetscErrorCode DMPlexGetIsoperiodicFaceSF(DM, PetscInt *, const PetscSF **);
diff --git a/src/dm/impls/plex/plexcreate.c b/src/dm/impls/plex/plexcreate.c
index 40352d7..a289c25 100644
--- a/src/dm/impls/plex/plexcreate.c
+++ b/src/dm/impls/plex/plexcreate.c
@@ -2051,7 +2051,7 @@ static PetscErrorCode DMPlexCreateWedgeBoxMesh_Internal(DM dm, const PetscInt fa
PetscCall(DMSetDimension(bdm, 2));
PetscCall(PetscLogEventBegin(DMPLEX_Generate, bdm, 0, 0, 0));
PetscCall(DMPlexCreateBoxMesh_Simplex_Internal(bdm, 2, faces, lower, upper, periodicity, PETSC_TRUE));
- PetscCall(DMPlexExtrude(bdm, faces[2], upper[2] - lower[2], PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, NULL, NULL, &vol));
+ PetscCall(DMPlexExtrude(bdm, faces[2], upper[2] - lower[2], PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, NULL, NULL, NULL, &vol));
PetscCall(PetscLogEventEnd(DMPLEX_Generate, bdm, 0, 0, 0));
PetscCall(DMDestroy(&bdm));
PetscCall(DMPlexReplace_Internal(dm, &vol));
diff --git a/src/dm/impls/plex/plexextrude.c b/src/dm/impls/plex/plexextrude.c
index b8c3fbb..8c2add4 100644
--- a/src/dm/impls/plex/plexextrude.c
+++ b/src/dm/impls/plex/plexextrude.c
@@ -48,7 +48,7 @@
.seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMExtrude()`, `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeSetTensor()`
@*/
-PetscErrorCode DMPlexExtrude(DM dm, PetscInt layers, PetscReal thickness, PetscBool tensor, PetscBool symmetric, PetscBool periodic, const PetscReal normal[], const PetscReal thicknesses[], DM *edm)
+PetscErrorCode DMPlexExtrude(DM dm, PetscInt layers, PetscReal thickness, PetscBool tensor, PetscBool symmetric, PetscBool periodic, const PetscReal normal[], const PetscReal thicknesses[], DMLabel adaptLabel, DM *edm)
{
DMPlexTransform tr;
DM cdm;
@@ -74,6 +74,7 @@ PetscErrorCode DMPlexExtrude(DM dm, PetscInt layers, PetscReal thickness, PetscB
PetscCall(DMPlexTransformExtrudeSetPeriodic(tr, periodic));
if (normal) PetscCall(DMPlexTransformExtrudeSetNormal(tr, normal));
if (thicknesses) PetscCall(DMPlexTransformExtrudeSetThicknesses(tr, layers, thicknesses));
+ if (adaptLabel) PetscCall(DMPlexTransformSetActive(tr, adaptLabel));
PetscCall(DMPlexTransformSetFromOptions(tr));
PetscCall(PetscObjectSetOptions((PetscObject)tr, NULL));
PetscCall(DMPlexTransformSetUp(tr));
@@ -130,7 +131,7 @@ PetscErrorCode DMPlexExtrude(DM dm, PetscInt layers, PetscReal thickness, PetscB
PetscErrorCode DMExtrude_Plex(DM dm, PetscInt layers, DM *edm)
{
PetscFunctionBegin;
- PetscCall(DMPlexExtrude(dm, layers, PETSC_DETERMINE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, NULL, NULL, edm));
+ PetscCall(DMPlexExtrude(dm, layers, PETSC_DETERMINE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, NULL, NULL, NULL, edm));
PetscCall(DMSetMatType(*edm, dm->mattype));
PetscCall(DMViewFromOptions(*edm, NULL, "-check_extrude"));
PetscFunctionReturn(PETSC_SUCCESS);
#include <petscdmplex.h>
#include <petscdmforest.h>
#include <petscdmplextransform.h>
PetscErrorCode extrudeMesh(DM dm, DM *dmExtrude) {
DMLabel extrudeLabel;
PetscInt fStart, fEnd, vStart, vEnd;
PetscInt nLayers = 2;
PetscReal thickness = 0.1;
PetscFunctionBegin;
PetscCall(DMLabelCreate(PETSC_COMM_SELF, "extrudeLabel", &extrudeLabel));
PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); // Faces
PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); // Vertices
for (PetscInt f = fStart; f < fEnd; ++f) {
PetscInt suppSize, numPoints;
PetscInt *points = NULL;
// Only use faces attached to one cell
PetscCall(DMPlexGetSupportSize(dm, f, &suppSize));
if (suppSize > 1) continue;
// Get everything at a lower dimension than the face
PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points));
for (PetscInt cl = 0; cl < numPoints * 2; cl += 2) {
const PetscInt v = points[cl];
if ((v >= vStart) && (v < vEnd)) {
PetscCall(DMLabelSetValue(extrudeLabel, v, DM_ADAPT_REFINE));
}
}
PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points));
}
PetscCall(DMPlexExtrude(dm, nLayers, thickness, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, NULL, NULL, extrudeLabel, dmExtrude));
PetscCall(DMLabelDestroy(&extrudeLabel));
PetscFunctionReturn(PETSC_SUCCESS);
}
int main(int argc,char **argv)
{
PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
DM dmPlex, dmForest, dmExtrude;
PetscInt dim = 2;
PetscInt pStart, pEnd;
PetscInt faces[2] = {10, 10};
PetscReal lower[2] = {0, 0}, upper[2] = {1, 1};
PetscBool interpolate = PETSC_TRUE;
DMBoundaryType periodicity[2] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
// Base plex
PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, faces, lower, upper, periodicity, interpolate, dim, PETSC_TRUE, &dmPlex));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original Plex:\n"));
for (PetscInt d = 0; d <= dim; ++d) {
PetscCall(DMPlexGetDepthStratum(dmPlex, d, &pStart, &pEnd));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t%d: [%d,%d)\n", d, pStart, pEnd));
}
PetscCall(DMGetBoundingBox(dmPlex, lower, upper));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Before Extrude Box: %+e, %+e\t%+e, %+e\n", lower[0], upper[0], lower[1], upper[1]));
PetscCall(extrudeMesh(dmPlex, &dmExtrude));
PetscCall(DMGetBoundingBox(dmExtrude, lower, upper));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, " After Extrude Box: %+e, %+e\t%+e, %+e\n", lower[0], upper[0], lower[1], upper[1]));
PetscCall(DMDestroy(&dmExtrude));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n"));
// forest
PetscCall(DMCreate(PETSC_COMM_WORLD, &dmForest));
PetscCall(DMSetType(dmForest, DMP4EST));
PetscCall(DMForestSetBaseDM(dmForest, dmPlex));
PetscCall(DMSetUp(dmForest));
PetscCall(DMDestroy(&dmPlex));
PetscCall(DMConvert(dmForest, DMPLEX, &dmPlex));
PetscCall(DMDestroy(&dmForest));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Converted Plex:\n"));
for (PetscInt d = 0; d <= dim; ++d) {
PetscCall(DMPlexGetDepthStratum(dmPlex, d, &pStart, &pEnd));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t%d: [%d,%d)\n", d, pStart, pEnd));
}
PetscCall(DMGetBoundingBox(dmPlex, lower, upper));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Before Extrude Box: %+e, %+e\t%+e, %+e\n", lower[0], upper[0], lower[1], upper[1]));
PetscCall(extrudeMesh(dmPlex, &dmExtrude));
PetscCall(DMGetBoundingBox(dmExtrude, lower, upper));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, " After Extrude Box: %+e, %+e\t%+e, %+e\n", lower[0], upper[0], lower[1], upper[1]));
PetscCall(DMDestroy(&dmExtrude));
PetscCall(DMDestroy(&dmPlex));
PetscCall(PetscFinalize());
return PETSC_SUCCESS;
}
