Actual source code: ex71.c

  1: static char help[] = "Tests for submesh creation\n\n";

  3: #include <petscdmplex.h>
  4: #include <petscsf.h>

  6: /* Submesh of a 2 x 2 mesh using 3 processes.

  8: Local numbering on each rank:

 10:      (6)(16)-(7)(17)-(8)       5--14---6--15---7      (10)(20)(11)(21)(12)
 11:       |       |       |        |       |       |        |       |       |
 12:     (18) (1)(19) (2)(20)      16   0  17   1  18      (22) (2)(23) (3)(24)
 13:       |       |       |        |       |       |        |       |       |
 14:      (5)(15)(11)(22)(12)       4--13-(11)(22)(12)      (9)(19)--6--14---7
 15:       |       |       |        |       |       |        |       |       |
 16:      14   0 (23) (3)(24)     (20) (2)(23) (3)(24)     (18) (1) 15   0  16
 17:       |       |       |        |       |       |        |       |       |
 18:       4--13--(9)(21)(10)      (8)(19)-(9)(21)(10)      (8)(17)--4--13---5

 20:            mesh_0                   mesh_1                   mesh_2

 22: where () represents ghost points. We extract the left 2 cells.
 23: With sanitize_submesh = PETSC_FALSE, we get:

 25:      (4)(11)-(5)               3---9---4
 26:       |       |                |       |
 27:     (12) (1)(13)              10   0  11
 28:       |       |                |       |
 29:      (3)(10) (7)               2---8---7
 30:       |       |                |       |
 31:       9   0 (14)             (13) (1) 14
 32:       |       |                |       |
 33:       2---8--(6)              (5)(12)--6

 35: On the other hand, with sanitize_submesh = PETSC_TRUE, we get:

 37:      (4)(11)-(5)               3---9---4
 38:       |       |                |       |
 39:     (12) (1)(13)              10   0  11
 40:       |       |                |       |
 41:      (3)(10) (7)               2---8---7
 42:       |       |                |       |
 43:       9   0  14              (13) (1)(14)
 44:       |       |                |       |
 45:       2---8---6               (5)(12)-(6)

 47:         submesh_0                submesh_1               submesh_2

 49: as points 15 and 4 of mesh_2 are in the closure of a submesh cell owned by rank 0 (point 0 of submesh_0),
 50: and not in the closure of any submesh cell owned by rank 1.

 52: */

 54: typedef struct {
 55:   PetscBool ignoreLabelHalo; /* Ignore filter values in the halo. */
 56:   PetscBool sanitizeSubmesh; /* Sanitize submesh. */
 57: } AppCtx;

 59: PetscErrorCode ProcessOptions(AppCtx *options)
 60: {
 61:   PetscFunctionBegin;
 62:   options->ignoreLabelHalo = PETSC_FALSE;
 63:   options->sanitizeSubmesh = PETSC_FALSE;

 65:   PetscOptionsBegin(PETSC_COMM_SELF, "", "Filtering Problem Options", "DMPLEX");
 66:   PetscCall(PetscOptionsBool("-ignore_label_halo", "Ignore filter values in the halo", "ex80.c", options->ignoreLabelHalo, &options->ignoreLabelHalo, NULL));
 67:   PetscCall(PetscOptionsBool("-sanitize_submesh", "Sanitize submesh", "ex80.c", options->sanitizeSubmesh, &options->sanitizeSubmesh, NULL));
 68:   PetscOptionsEnd();
 69:   PetscFunctionReturn(PETSC_SUCCESS);
 70: }
 71: int main(int argc, char **argv)
 72: {
 73:   DM             dm, subdm;
 74:   PetscSF        ownershipTransferSF;
 75:   DMLabel        filter;
 76:   const PetscInt filterValue = 1;
 77:   MPI_Comm       comm;
 78:   PetscMPIInt    size, rank;
 79:   AppCtx         user;

 81:   PetscFunctionBeginUser;
 82:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 83:   PetscCall(ProcessOptions(&user));
 84:   comm = PETSC_COMM_WORLD;
 85:   PetscCallMPI(MPI_Comm_size(comm, &size));
 86:   if (size != 3) {
 87:     PetscCall(PetscPrintf(comm, "This example is specifically designed for comm size == 3.\n"));
 88:     PetscCall(PetscFinalize());
 89:     return 0;
 90:   }
 91:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
 92:   {
 93:     DM             pdm;
 94:     const PetscInt faces[2] = {2, 2};
 95:     PetscInt       overlap  = 1;

 97:     PetscCall(DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &dm));
 98:     {
 99:       PetscPartitioner part;
100:       PetscInt        *sizes  = NULL;
101:       PetscInt        *points = NULL;

103:       if (rank == 0) {
104:         PetscInt sizes1[3]  = {1, 2, 1};
105:         PetscInt points1[4] = {0, 2, 3, 1};

107:         PetscCall(PetscMalloc2(3, &sizes, 4, &points));
108:         PetscCall(PetscArraycpy(sizes, sizes1, 3));
109:         PetscCall(PetscArraycpy(points, points1, 4));
110:       }
111:       PetscCall(DMPlexGetPartitioner(dm, &part));
112:       PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL));
113:       PetscCall(PetscPartitionerShellSetPartition(part, size, sizes, points));
114:       PetscCall(PetscFree2(sizes, points));
115:     }
116:     PetscCall(DMSetAdjacency(dm, -1, PETSC_FALSE, PETSC_TRUE));
117:     PetscCall(DMPlexDistribute(dm, overlap, NULL, &pdm));
118:     if (pdm) {
119:       PetscCall(DMDestroy(&dm));
120:       dm = pdm;
121:     }
122:   }
123:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "filter", &filter));
124:   switch (rank) {
125:   case 0:
126:     PetscCall(DMLabelSetValue(filter, 0, filterValue));
127:     PetscCall(DMLabelSetValue(filter, 1, filterValue));
128:     break;
129:   case 1:
130:     PetscCall(DMLabelSetValue(filter, 0, filterValue));
131:     PetscCall(DMLabelSetValue(filter, 2, filterValue));
132:     break;
133:   case 2:
134:     break;
135:   }
136:   PetscCall(PetscObjectSetName((PetscObject)dm, "Example_DM"));
137:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
138:   PetscCall(DMPlexFilter(dm, filter, filterValue, user.ignoreLabelHalo, user.sanitizeSubmesh, &ownershipTransferSF, &subdm));
139:   PetscCall(DMLabelDestroy(&filter));
140:   PetscCall(DMDestroy(&dm));
141:   PetscCall(PetscObjectSetName((PetscObject)subdm, "Example_SubDM"));
142:   PetscCall(DMViewFromOptions(subdm, NULL, "-dm_view"));
143:   PetscCall(DMDestroy(&subdm));
144:   PetscCall(PetscObjectSetName((PetscObject)ownershipTransferSF, "Example_Ownership_Transfer_SF"));
145:   PetscCall(PetscSFView(ownershipTransferSF, PETSC_VIEWER_STDOUT_WORLD));
146:   PetscCall(PetscSFDestroy(&ownershipTransferSF));
147:   PetscCall(PetscFinalize());
148:   return 0;
149: }

151: /*TEST

153:   testset:
154:     nsize: 3
155:     args: -dm_view ascii::ascii_info_detail

157:     test:
158:       suffix: 0
159:       args:

161:     test:
162:       suffix: 1
163:       args: -sanitize_submesh

165:     test:
166:       suffix: 2
167:       args: -ignore_label_halo

169: TEST*/