static char help[] = "Tests for parallel mesh loading and parallel topological interpolation\n\n"; #include /* List of test meshes Network ------- Test 0 (2 ranks): network=0: --------- cell 0 cell 1 cell 2 nCells-1 (edge) 0 ------ 1 ------ 2 ------ 3 -- -- v -- -- nCells (vertex) vertex distribution: rank 0: 0 1 rank 1: 2 3 ... nCells cell(edge) distribution: rank 0: 0 1 rank 1: 2 ... nCells-1 network=1: --------- v2 ^ | cell 2 | v0 --cell 0--> v3--cell 1--> v1 vertex distribution: rank 0: 0 1 3 rank 1: 2 cell(edge) distribution: rank 0: 0 1 rank 1: 2 example: mpiexec -n 2 ./ex18 -distribute 1 -dim 1 -orig_dm_view -dist_dm_view -dist_dm_view -petscpartitioner_type parmetis -ncells 50 Triangle -------- Test 0 (2 ranks): Two triangles sharing a face 2 / | \ / | \ / | \ 0 0 | 1 3 \ | / \ | / \ | / 1 vertex distribution: rank 0: 0 1 rank 1: 2 3 cell distribution: rank 0: 0 rank 1: 1 Test 1 (3 ranks): Four triangles partitioned across 3 ranks 0 _______ 3 | \ / | | \ 1 / | | \ / | | 0 2 2 | | / \ | | / 3 \ | | / \ | 1 ------- 4 vertex distribution: rank 0: 0 1 rank 1: 2 3 rank 2: 4 cell distribution: rank 0: 0 rank 1: 1 rank 2: 2 3 Test 2 (3 ranks): Four triangles partitioned across 3 ranks 1 _______ 3 | \ / | | \ 1 / | | \ / | | 0 0 2 | | / \ | | / 3 \ | | / \ | 2 ------- 4 vertex distribution: rank 0: 0 1 rank 1: 2 3 rank 2: 4 cell distribution: rank 0: 0 rank 1: 1 rank 2: 2 3 Tetrahedron ----------- Test 0: Two tets sharing a face cell 3 _______ cell 0 / | \ \ 1 / | \ \ / | \ \ 0----|----4-----2 \ | / / \ | / / \ | / / 1------- y | x |/ *----z vertex distribution: rank 0: 0 1 rank 1: 2 3 4 cell distribution: rank 0: 0 rank 1: 1 Quadrilateral ------------- Test 0 (2 ranks): Two quads sharing a face 3-------2-------5 | | | | 0 | 1 | | | | 0-------1-------4 vertex distribution: rank 0: 0 1 2 rank 1: 3 4 5 cell distribution: rank 0: 0 rank 1: 1 TODO Test 1: A quad and a triangle sharing a face 5-------4 | | \ | 0 | \ | | 1 \ 2-------3----6 Hexahedron ---------- Test 0 (2 ranks): Two hexes sharing a face cell 7-------------6-------------11 cell 0 /| /| /| 1 / | F1 / | F7 / | / | / | / | 4-------------5-------------10 | | | F4 | | F10 | | | | | | | | |F5 | |F3 | |F9 | | | F2 | | F8 | | | 3---------|---2---------|---9 | / | / | / | / F0 | / F6 | / |/ |/ |/ 0-------------1-------------8 vertex distribution: rank 0: 0 1 2 3 4 5 rank 1: 6 7 8 9 10 11 cell distribution: rank 0: 0 rank 1: 1 */ typedef enum {NONE, CREATE, AFTER_CREATE, AFTER_DISTRIBUTE} InterpType; typedef struct { PetscInt debug; /* The debugging level */ PetscInt testNum; /* Indicates the mesh to create */ PetscInt dim; /* The topological mesh dimension */ PetscBool cellSimplex; /* Use simplices or hexes */ PetscBool distribute; /* Distribute the mesh */ InterpType interpolate; /* Interpolate the mesh before or after DMPlexDistribute() */ PetscBool useGenerator; /* Construct mesh with a mesh generator */ PetscBool testOrientIF; /* Test for different original interface orientations */ PetscBool testHeavy; /* Run the heavy PointSF test */ PetscBool customView; /* Show results of DMPlexIsInterpolated() etc. */ PetscInt ornt[2]; /* Orientation of interface on rank 0 and rank 1 */ PetscInt faces[3]; /* Number of faces per dimension for generator */ PetscScalar coords[128]; PetscReal coordsTol; PetscInt ncoords; PetscInt pointsToExpand[128]; PetscInt nPointsToExpand; PetscBool testExpandPointsEmpty; char filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */ } AppCtx; struct _n_PortableBoundary { Vec coordinates; PetscInt depth; PetscSection *sections; }; typedef struct _n_PortableBoundary * PortableBoundary; #if defined(PETSC_USE_LOG) static PetscLogStage stage[3]; #endif static PetscErrorCode DMPlexCheckPointSFHeavy(DM, PortableBoundary); static PetscErrorCode DMPlexSetOrientInterface_Private(DM,PetscBool); static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM, PortableBoundary *); static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM, IS, PetscSection, IS *); static PetscErrorCode PortableBoundaryDestroy(PortableBoundary *bnd) { PetscInt d; PetscFunctionBegin; if (!*bnd) PetscFunctionReturn(0); PetscCall(VecDestroy(&(*bnd)->coordinates)); for (d=0; d < (*bnd)->depth; d++) { PetscCall(PetscSectionDestroy(&(*bnd)->sections[d])); } PetscCall(PetscFree((*bnd)->sections)); PetscCall(PetscFree(*bnd)); PetscFunctionReturn(0); } static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) { const char *interpTypes[4] = {"none", "create", "after_create", "after_distribute"}; PetscInt interp=NONE, dim; PetscBool flg1, flg2; PetscErrorCode ierr; PetscFunctionBegin; options->debug = 0; options->testNum = 0; options->dim = 2; options->cellSimplex = PETSC_TRUE; options->distribute = PETSC_FALSE; options->interpolate = NONE; options->useGenerator = PETSC_FALSE; options->testOrientIF = PETSC_FALSE; options->testHeavy = PETSC_TRUE; options->customView = PETSC_FALSE; options->testExpandPointsEmpty = PETSC_FALSE; options->ornt[0] = 0; options->ornt[1] = 0; options->faces[0] = 2; options->faces[1] = 2; options->faces[2] = 2; options->filename[0] = '\0'; options->coordsTol = PETSC_DEFAULT; ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");PetscCall(ierr); PetscCall(PetscOptionsBoundedInt("-debug", "The debugging level", "ex18.c", options->debug, &options->debug, NULL,0)); PetscCall(PetscOptionsBoundedInt("-testnum", "The mesh to create", "ex18.c", options->testNum, &options->testNum, NULL,0)); PetscCall(PetscOptionsBool("-cell_simplex", "Generate simplices if true, otherwise hexes", "ex18.c", options->cellSimplex, &options->cellSimplex, NULL)); PetscCall(PetscOptionsBool("-distribute", "Distribute the mesh", "ex18.c", options->distribute, &options->distribute, NULL)); PetscCall(PetscOptionsEList("-interpolate", "Type of mesh interpolation (none, create, after_create, after_distribute)", "ex18.c", interpTypes, 4, interpTypes[options->interpolate], &interp, NULL)); options->interpolate = (InterpType) interp; PetscCheckFalse(!options->distribute && options->interpolate == AFTER_DISTRIBUTE,comm, PETSC_ERR_SUP, "-interpolate after_distribute needs -distribute 1"); PetscCall(PetscOptionsBool("-use_generator", "Use a mesh generator to build the mesh", "ex18.c", options->useGenerator, &options->useGenerator, NULL)); options->ncoords = 128; PetscCall(PetscOptionsScalarArray("-view_vertices_from_coords", "Print DAG points corresponding to vertices with given coordinates", "ex18.c", options->coords, &options->ncoords, NULL)); PetscCall(PetscOptionsReal("-view_vertices_from_coords_tol", "Tolerance for -view_vertices_from_coords", "ex18.c", options->coordsTol, &options->coordsTol, NULL)); options->nPointsToExpand = 128; PetscCall(PetscOptionsIntArray("-test_expand_points", "Expand given array of DAG point using DMPlexGetConeRecursive() and print results", "ex18.c", options->pointsToExpand, &options->nPointsToExpand, NULL)); if (options->nPointsToExpand) { PetscCall(PetscOptionsBool("-test_expand_points_empty", "For -test_expand_points, rank 0 will have empty input array", "ex18.c", options->testExpandPointsEmpty, &options->testExpandPointsEmpty, NULL)); } PetscCall(PetscOptionsBool("-test_heavy", "Run the heavy PointSF test", "ex18.c", options->testHeavy, &options->testHeavy, NULL)); PetscCall(PetscOptionsBool("-custom_view", "Custom DMPlex view", "ex18.c", options->customView, &options->customView, NULL)); PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex18.c", options->dim, &options->dim, &flg1,1,3)); dim = 3; PetscCall(PetscOptionsIntArray("-faces", "Number of faces per dimension", "ex18.c", options->faces, &dim, &flg2)); if (flg2) { PetscCheckFalse(flg1 && dim != options->dim,comm, PETSC_ERR_ARG_OUTOFRANGE, "specified -dim %D is not equal to length %D of -faces (note that -dim can be omitted)", options->dim, dim); options->dim = dim; } PetscCall(PetscOptionsString("-filename", "The mesh file", "ex18.c", options->filename, options->filename, sizeof(options->filename), NULL)); PetscCall(PetscOptionsBoundedInt("-rotate_interface_0", "Rotation (relative orientation) of interface on rank 0; implies -interpolate create -distribute 0", "ex18.c", options->ornt[0], &options->ornt[0], &options->testOrientIF,0)); PetscCall(PetscOptionsBoundedInt("-rotate_interface_1", "Rotation (relative orientation) of interface on rank 1; implies -interpolate create -distribute 0", "ex18.c", options->ornt[1], &options->ornt[1], &flg2,0)); PetscCheck(flg2 == options->testOrientIF,comm, PETSC_ERR_ARG_OUTOFRANGE, "neither or both -rotate_interface_0 -rotate_interface_1 must be set"); if (options->testOrientIF) { PetscInt i; for (i=0; i<2; i++) { if (options->ornt[i] >= 10) options->ornt[i] = -(options->ornt[i]-10); /* 11 12 13 become -1 -2 -3 */ } options->filename[0] = 0; options->useGenerator = PETSC_FALSE; options->dim = 3; options->cellSimplex = PETSC_TRUE; options->interpolate = CREATE; options->distribute = PETSC_FALSE; } ierr = PetscOptionsEnd();PetscCall(ierr); PetscFunctionReturn(0); } static PetscErrorCode CreateMesh_1D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) { PetscInt testNum = user->testNum; PetscMPIInt rank,size; PetscInt numCorners=2,i; PetscInt numCells,numVertices,network; PetscInt *cells; PetscReal *coords; PetscFunctionBegin; PetscCallMPI(MPI_Comm_rank(comm, &rank)); PetscCallMPI(MPI_Comm_size(comm, &size)); PetscCheck(size <= 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for <=2 processes",testNum); numCells = 3; PetscCall(PetscOptionsGetInt(NULL, NULL, "-ncells", &numCells, NULL)); PetscCheck(numCells >= 3,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test ncells must >=3",numCells); if (size == 1) { numVertices = numCells + 1; PetscCall(PetscMalloc2(2*numCells,&cells,2*numVertices,&coords)); for (i=0; idim, numCells, numVertices, numCorners, PETSC_FALSE, cells, user->dim, coords, dm)); PetscCall(PetscFree2(cells,coords)); PetscFunctionReturn(0); } network = 0; PetscCall(PetscOptionsGetInt(NULL, NULL, "-network_case", &network, NULL)); if (network == 0) { switch (rank) { case 0: { numCells = 2; numVertices = numCells; PetscCall(PetscMalloc2(2*numCells,&cells,2*numCells,&coords)); cells[0] = 0; cells[1] = 1; cells[2] = 1; cells[3] = 2; coords[0] = 0.; coords[1] = 1.; coords[2] = 1.; coords[3] = 2.; } break; case 1: { numCells -= 2; numVertices = numCells + 1; PetscCall(PetscMalloc2(2*numCells,&cells,2*numCells,&coords)); for (i=0; idim, numCells, numVertices, PETSC_DECIDE, numCorners, PETSC_FALSE, cells, user->dim, coords, NULL, NULL, dm)); PetscCall(PetscFree2(cells,coords)); PetscFunctionReturn(0); } static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) { PetscInt testNum = user->testNum, p; PetscMPIInt rank, size; PetscFunctionBegin; PetscCallMPI(MPI_Comm_rank(comm, &rank)); PetscCallMPI(MPI_Comm_size(comm, &size)); switch (testNum) { case 0: PetscCheck(size == 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); switch (rank) { case 0: { const PetscInt numCells = 1, numVertices = 2, numCorners = 3; const PetscInt cells[3] = {0, 1, 2}; PetscReal coords[4] = {-0.5, 0.5, 0.0, 0.0}; PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; case 1: { const PetscInt numCells = 1, numVertices = 2, numCorners = 3; const PetscInt cells[3] = {1, 3, 2}; PetscReal coords[4] = {0.0, 1.0, 0.5, 0.5}; PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); } break; case 1: PetscCheck(size == 3,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 3 processes", testNum); switch (rank) { case 0: { const PetscInt numCells = 1, numVertices = 2, numCorners = 3; const PetscInt cells[3] = {0, 1, 2}; PetscReal coords[4] = {0.0, 1.0, 0.0, 0.0}; PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; case 1: { const PetscInt numCells = 1, numVertices = 2, numCorners = 3; const PetscInt cells[3] = {0, 2, 3}; PetscReal coords[4] = {0.5, 0.5, 1.0, 1.0}; PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; case 2: { const PetscInt numCells = 2, numVertices = 1, numCorners = 3; const PetscInt cells[6] = {2, 4, 3, 2, 1, 4}; PetscReal coords[2] = {1.0, 0.0}; PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); } break; case 2: PetscCheck(size == 3,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 3 processes", testNum); switch (rank) { case 0: { const PetscInt numCells = 1, numVertices = 2, numCorners = 3; const PetscInt cells[3] = {1, 2, 0}; PetscReal coords[4] = {0.5, 0.5, 0.0, 1.0}; PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; case 1: { const PetscInt numCells = 1, numVertices = 2, numCorners = 3; const PetscInt cells[3] = {1, 0, 3}; PetscReal coords[4] = {0.0, 0.0, 1.0, 1.0}; PetscInt markerPoints[6] = {1, 1, 2, 1, 3, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; case 2: { const PetscInt numCells = 2, numVertices = 1, numCorners = 3; const PetscInt cells[6] = {0, 4, 3, 0, 2, 4}; PetscReal coords[2] = {1.0, 0.0}; PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); } break; default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); } PetscFunctionReturn(0); } static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) { PetscInt testNum = user->testNum, p; PetscMPIInt rank, size; PetscFunctionBegin; PetscCallMPI(MPI_Comm_rank(comm, &rank)); PetscCallMPI(MPI_Comm_size(comm, &size)); switch (testNum) { case 0: PetscCheck(size == 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); switch (rank) { case 0: { const PetscInt numCells = 1, numVertices = 2, numCorners = 4; const PetscInt cells[4] = {0, 2, 1, 3}; PetscReal coords[6] = {0.0, 0.0, -0.5, 0.0, -0.5, 0.0}; PetscInt markerPoints[8] = {1, 1, 2, 1, 3, 1, 4, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; case 1: { const PetscInt numCells = 1, numVertices = 3, numCorners = 4; const PetscInt cells[4] = {1, 2, 4, 3}; PetscReal coords[9] = {1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5}; PetscInt markerPoints[8] = {1, 1, 2, 1, 3, 1, 4, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); } break; default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); } if (user->testOrientIF) { PetscInt ifp[] = {8, 6}; PetscCall(PetscObjectSetName((PetscObject) *dm, "Mesh before orientation")); PetscCall(DMViewFromOptions(*dm, NULL, "-before_orientation_dm_view")); /* rotate interface face ifp[rank] by given orientation ornt[rank] */ PetscCall(DMPlexOrientPoint(*dm, ifp[rank], user->ornt[rank])); PetscCall(DMViewFromOptions(*dm, NULL, "-before_orientation_dm_view")); PetscCall(DMPlexCheckFaces(*dm, 0)); PetscCall(DMPlexOrientInterface_Internal(*dm)); PetscCall(PetscPrintf(comm, "Orientation test PASSED\n")); } PetscFunctionReturn(0); } static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) { PetscInt testNum = user->testNum, p; PetscMPIInt rank, size; PetscFunctionBegin; PetscCallMPI(MPI_Comm_rank(comm, &rank)); PetscCallMPI(MPI_Comm_size(comm, &size)); switch (testNum) { case 0: PetscCheck(size == 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); switch (rank) { case 0: { const PetscInt numCells = 1, numVertices = 3, numCorners = 4; const PetscInt cells[4] = {0, 1, 2, 3}; PetscReal coords[6] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0}; PetscInt markerPoints[4*2] = {1, 1, 2, 1, 3, 1, 4, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; case 1: { const PetscInt numCells = 1, numVertices = 3, numCorners = 4; const PetscInt cells[4] = {1, 4, 5, 2}; PetscReal coords[6] = {-0.5, 1.0, 0.5, 0.0, 0.5, 1.0}; PetscInt markerPoints[4*2] = {1, 1, 2, 1, 3, 1, 4, 1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); } break; default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); } PetscFunctionReturn(0); } static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscBool interpolate, AppCtx *user, DM *dm) { PetscInt testNum = user->testNum, p; PetscMPIInt rank, size; PetscFunctionBegin; PetscCallMPI(MPI_Comm_rank(comm, &rank)); PetscCallMPI(MPI_Comm_size(comm, &size)); switch (testNum) { case 0: PetscCheck(size == 2,comm, PETSC_ERR_ARG_OUTOFRANGE, "Test mesh %d only for 2 processes", testNum); switch (rank) { case 0: { const PetscInt numCells = 1, numVertices = 6, numCorners = 8; const PetscInt cells[8] = {0, 3, 2, 1, 4, 5, 6, 7}; PetscReal coords[6*3] = {-0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0, -0.5,0.0,1.0, 0.0,0.0,1.0}; PetscInt markerPoints[8*2] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; case 1: { const PetscInt numCells = 1, numVertices = 6, numCorners = 8; const PetscInt cells[8] = {1, 2, 9, 8, 5, 10, 11, 6}; PetscReal coords[6*3] = {0.0,1.0,1.0, -0.5,1.0,1.0, 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0, 0.5,1.0,1.0}; PetscInt markerPoints[8*2] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, user->dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, cells, user->dim, coords, NULL, NULL, dm)); for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); } break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh for rank %d", rank); } break; default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %D", testNum); } PetscFunctionReturn(0); } static PetscErrorCode CustomView(DM dm, PetscViewer v) { DMPlexInterpolatedFlag interpolated; PetscBool distributed; PetscFunctionBegin; PetscCall(DMPlexIsDistributed(dm, &distributed)); PetscCall(DMPlexIsInterpolatedCollective(dm, &interpolated)); PetscCall(PetscViewerASCIIPrintf(v, "DMPlexIsDistributed: %s\n", PetscBools[distributed])); PetscCall(PetscViewerASCIIPrintf(v, "DMPlexIsInterpolatedCollective: %s\n", DMPlexInterpolatedFlags[interpolated])); PetscFunctionReturn(0); } static PetscErrorCode CreateMeshFromFile(MPI_Comm comm, AppCtx *user, DM *dm, DM *serialDM) { const char *filename = user->filename; PetscBool testHeavy = user->testHeavy; PetscBool interpCreate = user->interpolate == CREATE ? PETSC_TRUE : PETSC_FALSE; PetscBool distributed = PETSC_FALSE; PetscFunctionBegin; *serialDM = NULL; if (testHeavy && interpCreate) PetscCall(DMPlexSetOrientInterface_Private(NULL, PETSC_FALSE)); PetscCall(PetscLogStagePush(stage[0])); PetscCall(DMPlexCreateFromFile(comm, filename, "ex18_plex", interpCreate, dm)); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */ PetscCall(PetscLogStagePop()); if (testHeavy && interpCreate) PetscCall(DMPlexSetOrientInterface_Private(NULL, PETSC_TRUE)); PetscCall(DMPlexIsDistributed(*dm, &distributed)); PetscCall(PetscPrintf(comm, "DMPlexCreateFromFile produced %s mesh.\n", distributed ? "distributed" : "serial")); if (testHeavy && distributed) { PetscCall(PetscOptionsSetValue(NULL, "-dm_plex_hdf5_force_sequential", NULL)); PetscCall(DMPlexCreateFromFile(comm, filename, "ex18_plex", interpCreate, serialDM)); PetscCall(DMPlexIsDistributed(*serialDM, &distributed)); PetscCheck(!distributed,comm, PETSC_ERR_PLIB, "unable to create a serial DM from file"); } PetscCall(DMGetDimension(*dm, &user->dim)); PetscFunctionReturn(0); } static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) { PetscPartitioner part; PortableBoundary boundary = NULL; DM serialDM = NULL; PetscBool cellSimplex = user->cellSimplex; PetscBool useGenerator = user->useGenerator; PetscBool interpCreate = user->interpolate == CREATE ? PETSC_TRUE : PETSC_FALSE; PetscBool interpSerial = user->interpolate == AFTER_CREATE ? PETSC_TRUE : PETSC_FALSE; PetscBool interpParallel = user->interpolate == AFTER_DISTRIBUTE ? PETSC_TRUE : PETSC_FALSE; PetscBool testHeavy = user->testHeavy; PetscMPIInt rank; PetscFunctionBegin; PetscCallMPI(MPI_Comm_rank(comm, &rank)); if (user->filename[0]) { PetscCall(CreateMeshFromFile(comm, user, dm, &serialDM)); } else if (useGenerator) { PetscCall(PetscLogStagePush(stage[0])); PetscCall(DMPlexCreateBoxMesh(comm, user->dim, cellSimplex, user->faces, NULL, NULL, NULL, interpCreate, dm)); PetscCall(PetscLogStagePop()); } else { PetscCall(PetscLogStagePush(stage[0])); switch (user->dim) { case 1: PetscCall(CreateMesh_1D(comm, interpCreate, user, dm)); break; case 2: if (cellSimplex) { PetscCall(CreateSimplex_2D(comm, interpCreate, user, dm)); } else { PetscCall(CreateQuad_2D(comm, interpCreate, user, dm)); } break; case 3: if (cellSimplex) { PetscCall(CreateSimplex_3D(comm, interpCreate, user, dm)); } else { PetscCall(CreateHex_3D(comm, interpCreate, user, dm)); } break; default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", user->dim); } PetscCall(PetscLogStagePop()); } PetscCheck(user->ncoords % user->dim == 0,comm, PETSC_ERR_ARG_OUTOFRANGE, "length of coordinates array %D must be divisable by spatial dimension %D", user->ncoords, user->dim); PetscCall(PetscObjectSetName((PetscObject) *dm, "Original Mesh")); PetscCall(DMViewFromOptions(*dm, NULL, "-orig_dm_view")); if (interpSerial) { DM idm; if (testHeavy) PetscCall(DMPlexSetOrientInterface_Private(*dm, PETSC_FALSE)); PetscCall(PetscLogStagePush(stage[2])); PetscCall(DMPlexInterpolate(*dm, &idm)); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */ PetscCall(PetscLogStagePop()); if (testHeavy) PetscCall(DMPlexSetOrientInterface_Private(*dm, PETSC_TRUE)); PetscCall(DMDestroy(dm)); *dm = idm; PetscCall(PetscObjectSetName((PetscObject) *dm, "Interpolated Mesh")); PetscCall(DMViewFromOptions(*dm, NULL, "-intp_dm_view")); } /* Set partitioner options */ PetscCall(DMPlexGetPartitioner(*dm, &part)); if (part) { PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSIMPLE)); PetscCall(PetscPartitionerSetFromOptions(part)); } if (user->customView) PetscCall(CustomView(*dm, PETSC_VIEWER_STDOUT_(comm))); if (testHeavy) { PetscBool distributed; PetscCall(DMPlexIsDistributed(*dm, &distributed)); if (!serialDM && !distributed) { serialDM = *dm; PetscCall(PetscObjectReference((PetscObject)*dm)); } if (serialDM) { PetscCall(DMPlexGetExpandedBoundary_Private(serialDM, &boundary)); } if (boundary) { /* check DM which has been created in parallel and already interpolated */ PetscCall(DMPlexCheckPointSFHeavy(*dm, boundary)); } /* Orient interface because it could be deliberately skipped above. It is idempotent. */ PetscCall(DMPlexOrientInterface_Internal(*dm)); } if (user->distribute) { DM pdm = NULL; /* Redistribute mesh over processes using that partitioner */ PetscCall(PetscLogStagePush(stage[1])); PetscCall(DMPlexDistribute(*dm, 0, NULL, &pdm)); PetscCall(PetscLogStagePop()); if (pdm) { PetscCall(DMDestroy(dm)); *dm = pdm; PetscCall(PetscObjectSetName((PetscObject) *dm, "Redistributed Mesh")); PetscCall(DMViewFromOptions(*dm, NULL, "-dist_dm_view")); } if (interpParallel) { DM idm; if (testHeavy) PetscCall(DMPlexSetOrientInterface_Private(*dm, PETSC_FALSE)); PetscCall(PetscLogStagePush(stage[2])); PetscCall(DMPlexInterpolate(*dm, &idm)); /* with DMPlexOrientInterface_Internal() call skipped so that PointSF issues are left to DMPlexCheckPointSFHeavy() */ PetscCall(PetscLogStagePop()); if (testHeavy) PetscCall(DMPlexSetOrientInterface_Private(*dm, PETSC_TRUE)); PetscCall(DMDestroy(dm)); *dm = idm; PetscCall(PetscObjectSetName((PetscObject) *dm, "Interpolated Redistributed Mesh")); PetscCall(DMViewFromOptions(*dm, NULL, "-intp_dm_view")); } } if (testHeavy) { if (boundary) { PetscCall(DMPlexCheckPointSFHeavy(*dm, boundary)); } /* Orient interface because it could be deliberately skipped above. It is idempotent. */ PetscCall(DMPlexOrientInterface_Internal(*dm)); } PetscCall(PetscObjectSetName((PetscObject) *dm, "Parallel Mesh")); PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); PetscCall(DMSetFromOptions(*dm)); PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); if (user->customView) PetscCall(CustomView(*dm, PETSC_VIEWER_STDOUT_(comm))); PetscCall(DMDestroy(&serialDM)); PetscCall(PortableBoundaryDestroy(&boundary)); PetscFunctionReturn(0); } #define ps2d(number) ((double) PetscRealPart(number)) static inline PetscErrorCode coord2str(char buf[], size_t len, PetscInt dim, const PetscScalar coords[], PetscReal tol) { PetscFunctionBegin; PetscCheck(dim <= 3,PETSC_COMM_SELF, PETSC_ERR_SUP, "dim must be less than or equal 3"); if (tol >= 1e-3) { switch (dim) { case 1: PetscCall(PetscSNPrintf(buf,len,"(%12.3f)",ps2d(coords[0]))); case 2: PetscCall(PetscSNPrintf(buf,len,"(%12.3f, %12.3f)",ps2d(coords[0]),ps2d(coords[1]))); case 3: PetscCall(PetscSNPrintf(buf,len,"(%12.3f, %12.3f, %12.3f)",ps2d(coords[0]),ps2d(coords[1]),ps2d(coords[2]))); } } else { switch (dim) { case 1: PetscCall(PetscSNPrintf(buf,len,"(%12.6f)",ps2d(coords[0]))); case 2: PetscCall(PetscSNPrintf(buf,len,"(%12.6f, %12.6f)",ps2d(coords[0]),ps2d(coords[1]))); case 3: PetscCall(PetscSNPrintf(buf,len,"(%12.6f, %12.6f, %12.6f)",ps2d(coords[0]),ps2d(coords[1]),ps2d(coords[2]))); } } PetscFunctionReturn(0); } static PetscErrorCode ViewVerticesFromCoords(DM dm, Vec coordsVec, PetscReal tol, PetscViewer viewer) { PetscInt dim, i, npoints; IS pointsIS; const PetscInt *points; const PetscScalar *coords; char coordstr[128]; MPI_Comm comm; PetscMPIInt rank; PetscFunctionBegin; PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); PetscCallMPI(MPI_Comm_rank(comm, &rank)); PetscCall(DMGetDimension(dm, &dim)); PetscCall(PetscViewerASCIIPushSynchronized(viewer)); PetscCall(DMPlexFindVertices(dm, coordsVec, tol, &pointsIS)); PetscCall(ISGetIndices(pointsIS, &points)); PetscCall(ISGetLocalSize(pointsIS, &npoints)); PetscCall(VecGetArrayRead(coordsVec, &coords)); for (i=0; i < npoints; i++) { PetscCall(coord2str(coordstr, sizeof(coordstr), dim, &coords[i*dim], tol)); if (rank == 0 && i) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "-----\n")); PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %s --> points[%D] = %D\n", rank, coordstr, i, points[i])); PetscCall(PetscViewerFlush(viewer)); } PetscCall(PetscViewerASCIIPopSynchronized(viewer)); PetscCall(VecRestoreArrayRead(coordsVec, &coords)); PetscCall(ISRestoreIndices(pointsIS, &points)); PetscCall(ISDestroy(&pointsIS)); PetscFunctionReturn(0); } static PetscErrorCode TestExpandPoints(DM dm, AppCtx *user) { IS is; PetscSection *sects; IS *iss; PetscInt d,depth; PetscMPIInt rank; PetscViewer viewer=PETSC_VIEWER_STDOUT_WORLD, sviewer; PetscFunctionBegin; PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank)); if (user->testExpandPointsEmpty && rank == 0) { PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &is)); } else { PetscCall(ISCreateGeneral(PETSC_COMM_SELF, user->nPointsToExpand, user->pointsToExpand, PETSC_USE_POINTER, &is)); } PetscCall(DMPlexGetConeRecursive(dm, is, &depth, &iss, §s)); PetscCall(PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] ==========================\n",rank)); for (d=depth-1; d>=0; d--) { IS checkIS; PetscBool flg; PetscCall(PetscViewerASCIIPrintf(sviewer, "depth %D ---------------\n",d)); PetscCall(PetscSectionView(sects[d], sviewer)); PetscCall(ISView(iss[d], sviewer)); /* check reverse operation */ if (d < depth-1) { PetscCall(DMPlexExpandedConesToFaces_Private(dm, iss[d], sects[d], &checkIS)); PetscCall(ISEqualUnsorted(checkIS, iss[d+1], &flg)); PetscCheck(flg,PetscObjectComm((PetscObject) checkIS), PETSC_ERR_PLIB, "DMPlexExpandedConesToFaces_Private produced wrong IS"); PetscCall(ISDestroy(&checkIS)); } } PetscCall(PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); PetscCall(PetscViewerFlush(viewer)); PetscCall(DMPlexRestoreConeRecursive(dm, is, &depth, &iss, §s)); PetscCall(ISDestroy(&is)); PetscFunctionReturn(0); } static PetscErrorCode DMPlexExpandedConesToFaces_Private(DM dm, IS is, PetscSection section, IS *newis) { PetscInt n,n1,ncone,numCoveredPoints,o,p,q,start,end; const PetscInt *coveredPoints; const PetscInt *arr, *cone; PetscInt *newarr; PetscFunctionBegin; PetscCall(ISGetLocalSize(is, &n)); PetscCall(PetscSectionGetStorageSize(section, &n1)); PetscCall(PetscSectionGetChart(section, &start, &end)); PetscCheck(n == n1,PETSC_COMM_SELF, PETSC_ERR_PLIB, "IS size = %D != %D = section storage size", n, n1); PetscCall(ISGetIndices(is, &arr)); PetscCall(PetscMalloc1(end-start, &newarr)); for (q=start; q= 0) { PetscCall(DMPlexGetJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints)); PetscCheck(numCoveredPoints <= 1,PETSC_COMM_SELF, PETSC_ERR_PLIB, "more than one covered points for section point %D",q); if (numCoveredPoints) p = coveredPoints[0]; else p = -2; PetscCall(DMPlexRestoreJoin(dm, ncone, cone, &numCoveredPoints, &coveredPoints)); } } newarr[q-start] = p; } PetscCall(ISRestoreIndices(is, &arr)); PetscCall(ISCreateGeneral(PETSC_COMM_SELF, end-start, newarr, PETSC_OWN_POINTER, newis)); PetscFunctionReturn(0); } static PetscErrorCode DMPlexExpandedVerticesToFaces_Private(DM dm, IS boundary_expanded_is, PetscInt depth, PetscSection sections[], IS *boundary_is) { PetscInt d; IS is,newis; PetscFunctionBegin; is = boundary_expanded_is; PetscCall(PetscObjectReference((PetscObject)is)); for (d = 0; d < depth-1; ++d) { PetscCall(DMPlexExpandedConesToFaces_Private(dm, is, sections[d], &newis)); PetscCall(ISDestroy(&is)); is = newis; } *boundary_is = is; PetscFunctionReturn(0); } #define CHKERRQI(incall,ierr) if (ierr) {incall = PETSC_FALSE; } static PetscErrorCode DMLabelViewFromOptionsOnComm_Private(DMLabel label, const char optionname[], MPI_Comm comm) { PetscViewer viewer; PetscBool flg; static PetscBool incall = PETSC_FALSE; PetscViewerFormat format; PetscFunctionBegin; if (incall) PetscFunctionReturn(0); incall = PETSC_TRUE; CHKERRQI(incall,PetscOptionsGetViewer(comm,((PetscObject)label)->options,((PetscObject)label)->prefix,optionname,&viewer,&format,&flg)); if (flg) { CHKERRQI(incall,PetscViewerPushFormat(viewer,format)); CHKERRQI(incall,DMLabelView(label, viewer)); CHKERRQI(incall,PetscViewerPopFormat(viewer)); CHKERRQI(incall,PetscViewerDestroy(&viewer)); } incall = PETSC_FALSE; PetscFunctionReturn(0); } /* TODO: this is hotfixing DMLabelGetStratumIS() - it should be fixed systematically instead */ static inline PetscErrorCode DMLabelGetStratumISOnComm_Private(DMLabel label, PetscInt value, MPI_Comm comm, IS *is) { IS tmpis; PetscFunctionBegin; PetscCall(DMLabelGetStratumIS(label, value, &tmpis)); if (!tmpis) PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_USE_POINTER, &tmpis)); PetscCall(ISOnComm(tmpis, comm, PETSC_COPY_VALUES, is)); PetscCall(ISDestroy(&tmpis)); PetscFunctionReturn(0); } /* currently only for simple PetscSection without fields or constraints */ static PetscErrorCode PetscSectionReplicate_Private(MPI_Comm comm, PetscMPIInt rootrank, PetscSection sec0, PetscSection *secout) { PetscSection sec; PetscInt chart[2], p; PetscInt *dofarr; PetscMPIInt rank; PetscFunctionBegin; PetscCallMPI(MPI_Comm_rank(comm, &rank)); if (rank == rootrank) { PetscCall(PetscSectionGetChart(sec0, &chart[0], &chart[1])); } PetscCallMPI(MPI_Bcast(chart, 2, MPIU_INT, rootrank, comm)); PetscCall(PetscMalloc1(chart[1]-chart[0], &dofarr)); if (rank == rootrank) { for (p = chart[0]; p < chart[1]; p++) { PetscCall(PetscSectionGetDof(sec0, p, &dofarr[p-chart[0]])); } } PetscCallMPI(MPI_Bcast(dofarr, chart[1]-chart[0], MPIU_INT, rootrank, comm)); PetscCall(PetscSectionCreate(comm, &sec)); PetscCall(PetscSectionSetChart(sec, chart[0], chart[1])); for (p = chart[0]; p < chart[1]; p++) { PetscCall(PetscSectionSetDof(sec, p, dofarr[p-chart[0]])); } PetscCall(PetscSectionSetUp(sec)); PetscCall(PetscFree(dofarr)); *secout = sec; PetscFunctionReturn(0); } static PetscErrorCode DMPlexExpandedVerticesCoordinatesToFaces_Private(DM ipdm, PortableBoundary bnd, IS *face_is) { IS faces_expanded_is; PetscFunctionBegin; PetscCall(DMPlexFindVertices(ipdm, bnd->coordinates, 0.0, &faces_expanded_is)); PetscCall(DMPlexExpandedVerticesToFaces_Private(ipdm, faces_expanded_is, bnd->depth, bnd->sections, face_is)); PetscCall(ISDestroy(&faces_expanded_is)); PetscFunctionReturn(0); } /* hack disabling DMPlexOrientInterface() call in DMPlexInterpolate() via -dm_plex_interpolate_orient_interfaces option */ static PetscErrorCode DMPlexSetOrientInterface_Private(DM dm, PetscBool enable) { PetscOptions options = NULL; const char *prefix = NULL; const char opt[] = "-dm_plex_interpolate_orient_interfaces"; char prefix_opt[512]; PetscBool flg, set; static PetscBool wasSetTrue = PETSC_FALSE; PetscFunctionBegin; if (dm) { PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); options = ((PetscObject)dm)->options; } PetscCall(PetscStrcpy(prefix_opt, "-")); PetscCall(PetscStrlcat(prefix_opt, prefix, sizeof(prefix_opt))); PetscCall(PetscStrlcat(prefix_opt, &opt[1], sizeof(prefix_opt))); PetscCall(PetscOptionsGetBool(options, prefix, opt, &flg, &set)); if (!enable) { if (set && flg) wasSetTrue = PETSC_TRUE; PetscCall(PetscOptionsSetValue(options, prefix_opt, "0")); } else if (set && !flg) { if (wasSetTrue) { PetscCall(PetscOptionsSetValue(options, prefix_opt, "1")); } else { /* default is PETSC_TRUE */ PetscCall(PetscOptionsClearValue(options, prefix_opt)); } wasSetTrue = PETSC_FALSE; } if (PetscDefined(USE_DEBUG)) { PetscCall(PetscOptionsGetBool(options, prefix, opt, &flg, &set)); PetscCheckFalse(set && flg != enable,PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "PetscOptionsSetValue did not have the desired effect"); } PetscFunctionReturn(0); } /* get coordinate description of the whole-domain boundary */ static PetscErrorCode DMPlexGetExpandedBoundary_Private(DM dm, PortableBoundary *boundary) { PortableBoundary bnd0, bnd; MPI_Comm comm; DM idm; DMLabel label; PetscInt d; const char boundaryName[] = "DMPlexDistributeInterpolateMarkInterface_boundary"; IS boundary_is; IS *boundary_expanded_iss; PetscMPIInt rootrank = 0; PetscMPIInt rank, size; PetscInt value = 1; DMPlexInterpolatedFlag intp; PetscBool flg; PetscFunctionBegin; PetscCall(PetscNew(&bnd)); PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); PetscCallMPI(MPI_Comm_rank(comm, &rank)); PetscCallMPI(MPI_Comm_size(comm, &size)); PetscCall(DMPlexIsDistributed(dm, &flg)); PetscCheck(!flg,comm, PETSC_ERR_ARG_WRONG, "serial DM (all points on one rank) needed"); /* interpolate serial DM if not yet interpolated */ PetscCall(DMPlexIsInterpolatedCollective(dm, &intp)); if (intp == DMPLEX_INTERPOLATED_FULL) { idm = dm; PetscCall(PetscObjectReference((PetscObject)dm)); } else { PetscCall(DMPlexInterpolate(dm, &idm)); PetscCall(DMViewFromOptions(idm, NULL, "-idm_view")); } /* mark whole-domain boundary of the serial DM */ PetscCall(DMLabelCreate(PETSC_COMM_SELF, boundaryName, &label)); PetscCall(DMAddLabel(idm, label)); PetscCall(DMPlexMarkBoundaryFaces(idm, value, label)); PetscCall(DMLabelViewFromOptionsOnComm_Private(label, "-idm_boundary_view", comm)); PetscCall(DMLabelGetStratumIS(label, value, &boundary_is)); /* translate to coordinates */ PetscCall(PetscNew(&bnd0)); PetscCall(DMGetCoordinatesLocalSetUp(idm)); if (rank == rootrank) { PetscCall(DMPlexGetConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections)); PetscCall(DMGetCoordinatesLocalTuple(dm, boundary_expanded_iss[0], NULL, &bnd0->coordinates)); /* self-check */ { IS is0; PetscCall(DMPlexExpandedVerticesCoordinatesToFaces_Private(idm, bnd0, &is0)); PetscCall(ISEqual(is0, boundary_is, &flg)); PetscCheck(flg,PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMPlexExpandedVerticesCoordinatesToFaces_Private produced a wrong IS"); PetscCall(ISDestroy(&is0)); } } else { PetscCall(VecCreateSeq(PETSC_COMM_SELF, 0, &bnd0->coordinates)); } { Vec tmp; VecScatter sc; IS xis; PetscInt n; /* just convert seq vectors to mpi vector */ PetscCall(VecGetLocalSize(bnd0->coordinates, &n)); PetscCallMPI(MPI_Bcast(&n, 1, MPIU_INT, rootrank, comm)); if (rank == rootrank) { PetscCall(VecCreateMPI(comm, n, n, &tmp)); } else { PetscCall(VecCreateMPI(comm, 0, n, &tmp)); } PetscCall(VecCopy(bnd0->coordinates, tmp)); PetscCall(VecDestroy(&bnd0->coordinates)); bnd0->coordinates = tmp; /* replicate coordinates from root rank to all ranks */ PetscCall(VecCreateMPI(comm, n, n*size, &bnd->coordinates)); PetscCall(ISCreateStride(comm, n, 0, 1, &xis)); PetscCall(VecScatterCreate(bnd0->coordinates, xis, bnd->coordinates, NULL, &sc)); PetscCall(VecScatterBegin(sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD)); PetscCall(VecScatterEnd( sc, bnd0->coordinates, bnd->coordinates, INSERT_VALUES, SCATTER_FORWARD)); PetscCall(VecScatterDestroy(&sc)); PetscCall(ISDestroy(&xis)); } bnd->depth = bnd0->depth; PetscCallMPI(MPI_Bcast(&bnd->depth, 1, MPIU_INT, rootrank, comm)); PetscCall(PetscMalloc1(bnd->depth, &bnd->sections)); for (d=0; ddepth; d++) { PetscCall(PetscSectionReplicate_Private(comm, rootrank, (rank == rootrank) ? bnd0->sections[d] : NULL, &bnd->sections[d])); } if (rank == rootrank) { PetscCall(DMPlexRestoreConeRecursive(idm, boundary_is, &bnd0->depth, &boundary_expanded_iss, &bnd0->sections)); } PetscCall(PortableBoundaryDestroy(&bnd0)); PetscCall(DMRemoveLabelBySelf(idm, &label, PETSC_TRUE)); PetscCall(DMLabelDestroy(&label)); PetscCall(ISDestroy(&boundary_is)); PetscCall(DMDestroy(&idm)); *boundary = bnd; PetscFunctionReturn(0); } /* get faces of inter-partition interface */ static PetscErrorCode DMPlexGetInterfaceFaces_Private(DM ipdm, IS boundary_faces_is, IS *interface_faces_is) { MPI_Comm comm; DMLabel label; IS part_boundary_faces_is; const char partBoundaryName[] = "DMPlexDistributeInterpolateMarkInterface_partBoundary"; PetscInt value = 1; DMPlexInterpolatedFlag intp; PetscFunctionBegin; PetscCall(PetscObjectGetComm((PetscObject)ipdm, &comm)); PetscCall(DMPlexIsInterpolatedCollective(ipdm, &intp)); PetscCheck(intp == DMPLEX_INTERPOLATED_FULL,comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex"); /* get ipdm partition boundary (partBoundary) */ PetscCall(DMLabelCreate(PETSC_COMM_SELF, partBoundaryName, &label)); PetscCall(DMAddLabel(ipdm, label)); PetscCall(DMPlexMarkBoundaryFaces(ipdm, value, label)); PetscCall(DMLabelViewFromOptionsOnComm_Private(label, "-ipdm_part_boundary_view", comm)); PetscCall(DMLabelGetStratumISOnComm_Private(label, value, comm, &part_boundary_faces_is)); PetscCall(DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE)); PetscCall(DMLabelDestroy(&label)); /* remove ipdm whole-domain boundary (boundary_faces_is) from ipdm partition boundary (part_boundary_faces_is), resulting just in inter-partition interface */ PetscCall(ISDifference(part_boundary_faces_is,boundary_faces_is,interface_faces_is)); PetscCall(ISDestroy(&part_boundary_faces_is)); PetscFunctionReturn(0); } /* compute inter-partition interface including edges and vertices */ static PetscErrorCode DMPlexComputeCompleteInterface_Private(DM ipdm, IS interface_faces_is, IS *interface_is) { DMLabel label; PetscInt value = 1; const char interfaceName[] = "DMPlexDistributeInterpolateMarkInterface_interface"; DMPlexInterpolatedFlag intp; MPI_Comm comm; PetscFunctionBegin; PetscCall(PetscObjectGetComm((PetscObject)ipdm, &comm)); PetscCall(DMPlexIsInterpolatedCollective(ipdm, &intp)); PetscCheck(intp == DMPLEX_INTERPOLATED_FULL,comm, PETSC_ERR_ARG_WRONG, "only for fully interpolated DMPlex"); PetscCall(DMLabelCreate(PETSC_COMM_SELF, interfaceName, &label)); PetscCall(DMAddLabel(ipdm, label)); PetscCall(DMLabelSetStratumIS(label, value, interface_faces_is)); PetscCall(DMLabelViewFromOptionsOnComm_Private(label, "-interface_faces_view", comm)); PetscCall(DMPlexLabelComplete(ipdm, label)); PetscCall(DMLabelViewFromOptionsOnComm_Private(label, "-interface_view", comm)); PetscCall(DMLabelGetStratumISOnComm_Private(label, value, comm, interface_is)); PetscCall(PetscObjectSetName((PetscObject)*interface_is, "interface_is")); PetscCall(ISViewFromOptions(*interface_is, NULL, "-interface_is_view")); PetscCall(DMRemoveLabelBySelf(ipdm, &label, PETSC_TRUE)); PetscCall(DMLabelDestroy(&label)); PetscFunctionReturn(0); } static PetscErrorCode PointSFGetOutwardInterfacePoints(PetscSF sf, IS *is) { PetscInt n; const PetscInt *arr; PetscFunctionBegin; PetscCall(PetscSFGetGraph(sf, NULL, &n, &arr, NULL)); PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_USE_POINTER, is)); PetscFunctionReturn(0); } static PetscErrorCode PointSFGetInwardInterfacePoints(PetscSF sf, IS *is) { PetscInt n; const PetscInt *rootdegree; PetscInt *arr; PetscFunctionBegin; PetscCall(PetscSFSetUp(sf)); PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); PetscCall(PetscSFComputeMultiRootOriginalNumbering(sf, rootdegree, &n, &arr)); PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)sf), n, arr, PETSC_OWN_POINTER, is)); PetscFunctionReturn(0); } static PetscErrorCode PointSFGetInterfacePoints_Private(PetscSF pointSF, IS *is) { IS pointSF_out_is, pointSF_in_is; PetscFunctionBegin; PetscCall(PointSFGetOutwardInterfacePoints(pointSF, &pointSF_out_is)); PetscCall(PointSFGetInwardInterfacePoints(pointSF, &pointSF_in_is)); PetscCall(ISExpand(pointSF_out_is, pointSF_in_is, is)); PetscCall(ISDestroy(&pointSF_out_is)); PetscCall(ISDestroy(&pointSF_in_is)); PetscFunctionReturn(0); } #define CHKERRMY(ierr) PetscCheck(!ierr,PETSC_COMM_SELF, PETSC_ERR_PLIB, "PointSF is wrong. Unable to show details!") static PetscErrorCode ViewPointsWithType_Internal(DM dm, IS pointsIS, PetscViewer v) { DMLabel label; PetscSection coordsSection; Vec coordsVec; PetscScalar *coordsScalar; PetscInt coneSize, depth, dim, i, p, npoints; const PetscInt *points; PetscFunctionBegin; PetscCall(DMGetDimension(dm, &dim)); PetscCall(DMGetCoordinateSection(dm, &coordsSection)); PetscCall(DMGetCoordinatesLocal(dm, &coordsVec)); PetscCall(VecGetArray(coordsVec, &coordsScalar)); PetscCall(ISGetLocalSize(pointsIS, &npoints)); PetscCall(ISGetIndices(pointsIS, &points)); PetscCall(DMPlexGetDepthLabel(dm, &label)); PetscCall(PetscViewerASCIIPushTab(v)); for (i=0; i