Actual source code: ex11.c

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

  3: #include <petscdmplex.h>
  4: #include <petsc/private/dmimpl.h>

  6: static PetscErrorCode TestInsertion()
  7: {
  8:   DMLabel        label, label2;
  9:   const PetscInt values[5] = {0, 3, 4, -1, 176}, N = 10000;
 10:   PetscInt       i, v;

 12:   PetscFunctionBegin;
 13:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Test Label", &label));
 14:   PetscCall(DMLabelSetDefaultValue(label, -100));
 15:   for (i = 0; i < N; ++i) PetscCall(DMLabelSetValue(label, i, values[i % 5]));
 16:   /* Test get in hash mode */
 17:   for (i = 0; i < N; ++i) {
 18:     PetscInt val;

 20:     PetscCall(DMLabelGetValue(label, i, &val));
 21:     PetscCheck(val == values[i % 5], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Value %" PetscInt_FMT " for point %" PetscInt_FMT " should be %" PetscInt_FMT, val, i, values[i % 5]);
 22:   }
 23:   /* Test stratum */
 24:   for (v = 0; v < 5; ++v) {
 25:     IS              stratum;
 26:     const PetscInt *points;
 27:     PetscInt        n;

 29:     PetscCall(DMLabelGetStratumIS(label, values[v], &stratum));
 30:     PetscCheck(stratum, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Stratum %" PetscInt_FMT " is empty!", v);
 31:     PetscCall(ISGetIndices(stratum, &points));
 32:     PetscCall(ISGetLocalSize(stratum, &n));
 33:     for (i = 0; i < n; ++i) PetscCheck(points[i] == i * 5 + v, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " should be %" PetscInt_FMT, points[i], i * 5 + v);
 34:     PetscCall(ISRestoreIndices(stratum, &points));
 35:     PetscCall(ISDestroy(&stratum));
 36:   }
 37:   /* Test get in array mode */
 38:   for (i = 0; i < N; ++i) {
 39:     PetscInt val;

 41:     PetscCall(DMLabelGetValue(label, i, &val));
 42:     PetscCheck(val == values[i % 5], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Value %" PetscInt_FMT " should be %" PetscInt_FMT, val, values[i % 5]);
 43:   }
 44:   /* Test Duplicate */
 45:   PetscCall(DMLabelDuplicate(label, &label2));
 46:   for (i = 0; i < N; ++i) {
 47:     PetscInt val;

 49:     PetscCall(DMLabelGetValue(label2, i, &val));
 50:     PetscCheck(val == values[i % 5], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Value %" PetscInt_FMT " should be %" PetscInt_FMT, val, values[i % 5]);
 51:   }
 52:   PetscCall(DMLabelDestroy(&label2));
 53:   PetscCall(DMLabelDestroy(&label));
 54:   PetscFunctionReturn(PETSC_SUCCESS);
 55: }

 57: static PetscErrorCode TestEmptyStrata(MPI_Comm comm)
 58: {
 59:   DM               dm, dmDist;
 60:   PetscPartitioner part;
 61:   PetscInt         c0[6]  = {2, 3, 6, 7, 9, 11};
 62:   PetscInt         c1[6]  = {4, 5, 7, 8, 10, 12};
 63:   PetscInt         c2[4]  = {13, 15, 19, 21};
 64:   PetscInt         c3[4]  = {14, 16, 20, 22};
 65:   PetscInt         c4[4]  = {15, 17, 21, 23};
 66:   PetscInt         c5[4]  = {16, 18, 22, 24};
 67:   PetscInt         c6[4]  = {13, 14, 19, 20};
 68:   PetscInt         c7[4]  = {15, 16, 21, 22};
 69:   PetscInt         c8[4]  = {17, 18, 23, 24};
 70:   PetscInt         c9[4]  = {13, 14, 15, 16};
 71:   PetscInt         c10[4] = {15, 16, 17, 18};
 72:   PetscInt         c11[4] = {19, 20, 21, 22};
 73:   PetscInt         c12[4] = {21, 22, 23, 24};
 74:   PetscInt         dim    = 3;
 75:   PetscMPIInt      rank;

 77:   PetscFunctionBegin;
 78:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
 79:   /* A 3D box with two adjacent cells, sharing one face and four vertices */
 80:   PetscCall(DMCreate(comm, &dm));
 81:   PetscCall(DMSetType(dm, DMPLEX));
 82:   PetscCall(DMSetDimension(dm, dim));
 83:   if (rank == 0) {
 84:     PetscCall(DMPlexSetChart(dm, 0, 25));
 85:     PetscCall(DMPlexSetConeSize(dm, 0, 6));
 86:     PetscCall(DMPlexSetConeSize(dm, 1, 6));
 87:     PetscCall(DMPlexSetConeSize(dm, 2, 4));
 88:     PetscCall(DMPlexSetConeSize(dm, 3, 4));
 89:     PetscCall(DMPlexSetConeSize(dm, 4, 4));
 90:     PetscCall(DMPlexSetConeSize(dm, 5, 4));
 91:     PetscCall(DMPlexSetConeSize(dm, 6, 4));
 92:     PetscCall(DMPlexSetConeSize(dm, 7, 4));
 93:     PetscCall(DMPlexSetConeSize(dm, 8, 4));
 94:     PetscCall(DMPlexSetConeSize(dm, 9, 4));
 95:     PetscCall(DMPlexSetConeSize(dm, 10, 4));
 96:     PetscCall(DMPlexSetConeSize(dm, 11, 4));
 97:     PetscCall(DMPlexSetConeSize(dm, 12, 4));
 98:   }
 99:   PetscCall(DMSetUp(dm));
100:   if (rank == 0) {
101:     PetscCall(DMPlexSetCone(dm, 0, c0));
102:     PetscCall(DMPlexSetCone(dm, 1, c1));
103:     PetscCall(DMPlexSetCone(dm, 2, c2));
104:     PetscCall(DMPlexSetCone(dm, 3, c3));
105:     PetscCall(DMPlexSetCone(dm, 4, c4));
106:     PetscCall(DMPlexSetCone(dm, 5, c5));
107:     PetscCall(DMPlexSetCone(dm, 6, c6));
108:     PetscCall(DMPlexSetCone(dm, 7, c7));
109:     PetscCall(DMPlexSetCone(dm, 8, c8));
110:     PetscCall(DMPlexSetCone(dm, 9, c9));
111:     PetscCall(DMPlexSetCone(dm, 10, c10));
112:     PetscCall(DMPlexSetCone(dm, 11, c11));
113:     PetscCall(DMPlexSetCone(dm, 12, c12));
114:   }
115:   PetscCall(DMPlexSymmetrize(dm));
116:   /* Create a user managed depth label, so that we can leave out edges */
117:   {
118:     DMLabel  label;
119:     PetscInt numValues, maxValues = 0, v;

121:     PetscCall(DMCreateLabel(dm, "depth"));
122:     PetscCall(DMPlexGetDepthLabel(dm, &label));
123:     if (rank == 0) {
124:       PetscInt i;

126:       for (i = 0; i < 25; ++i) {
127:         if (i < 2) PetscCall(DMLabelSetValue(label, i, 3));
128:         else if (i < 13) PetscCall(DMLabelSetValue(label, i, 2));
129:         else {
130:           if (i == 13) PetscCall(DMLabelAddStratum(label, 1));
131:           PetscCall(DMLabelSetValue(label, i, 0));
132:         }
133:       }
134:     }
135:     PetscCall(DMLabelGetNumValues(label, &numValues));
136:     PetscCallMPI(MPIU_Allreduce(&numValues, &maxValues, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
137:     for (v = numValues; v < maxValues; ++v) PetscCall(DMLabelAddStratum(label, v));
138:   }
139:   {
140:     DMLabel label;
141:     PetscCall(DMPlexGetDepthLabel(dm, &label));
142:     PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_(comm)));
143:   }
144:   PetscCall(DMPlexGetPartitioner(dm, &part));
145:   PetscCall(PetscPartitionerSetFromOptions(part));
146:   PetscCall(DMPlexDistribute(dm, 1, NULL, &dmDist));
147:   if (dmDist) {
148:     PetscCall(DMDestroy(&dm));
149:     dm = dmDist;
150:   }
151:   {
152:     DMLabel label;
153:     PetscCall(DMPlexGetDepthLabel(dm, &label));
154:     PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_(comm)));
155:   }
156:   /* Create a cell vector */
157:   {
158:     Vec          v;
159:     PetscSection s;
160:     PetscInt     numComp[] = {1};
161:     PetscInt     dof[]     = {0, 0, 0, 1};
162:     PetscInt     N;

164:     PetscCall(DMSetNumFields(dm, 1));
165:     PetscCall(DMPlexCreateSection(dm, NULL, numComp, dof, 0, NULL, NULL, NULL, NULL, &s));
166:     PetscCall(DMSetLocalSection(dm, s));
167:     PetscCall(PetscSectionDestroy(&s));
168:     PetscCall(DMCreateGlobalVector(dm, &v));
169:     PetscCall(VecGetSize(v, &N));
170:     if (N != 2) {
171:       PetscCall(DMView(dm, PETSC_VIEWER_STDOUT_(comm)));
172:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "FAIL: Vector size %" PetscInt_FMT " != 2", N);
173:     }
174:     PetscCall(VecDestroy(&v));
175:   }
176:   PetscCall(DMDestroy(&dm));
177:   PetscFunctionReturn(PETSC_SUCCESS);
178: }

180: static PetscErrorCode TestDistribution(MPI_Comm comm)
181: {
182:   DM               dm, dmDist;
183:   PetscPartitioner part;
184:   DMLabel          label;
185:   char             filename[PETSC_MAX_PATH_LEN];
186:   const char      *name    = "test label";
187:   PetscInt         overlap = 0, cStart, cEnd, c;
188:   PetscMPIInt      rank;
189:   PetscBool        flg;

191:   PetscFunctionBegin;
192:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
193:   PetscCall(PetscOptionsGetString(NULL, NULL, "-filename", filename, sizeof(filename), &flg));
194:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
195:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-overlap", &overlap, NULL));
196:   PetscCall(DMPlexCreateFromFile(comm, filename, "ex11_plex", PETSC_TRUE, &dm));
197:   PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
198:   PetscCall(DMCreateLabel(dm, name));
199:   PetscCall(DMGetLabel(dm, name, &label));
200:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
201:   for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(label, c, c));
202:   PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD));
203:   PetscCall(DMPlexGetPartitioner(dm, &part));
204:   PetscCall(PetscPartitionerSetFromOptions(part));
205:   PetscCall(DMPlexDistribute(dm, overlap, NULL, &dmDist));
206:   if (dmDist) {
207:     PetscCall(DMDestroy(&dm));
208:     dm = dmDist;
209:   }
210:   PetscCall(PetscObjectSetName((PetscObject)dm, "Mesh"));
211:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
212:   PetscCall(DMGetLabel(dm, name, &label));
213:   PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD));
214:   PetscCall(DMDestroy(&dm));
215:   PetscFunctionReturn(PETSC_SUCCESS);
216: }

218: static PetscErrorCode TestUniversalLabel(MPI_Comm comm)
219: {
220:   DM               dm1, dm2;
221:   DMLabel          bd1, bd2, ulabel;
222:   DMUniversalLabel universal;
223:   PetscInt         pStart, pEnd, p;
224:   PetscBool        run = PETSC_FALSE, notFile;

226:   PetscFunctionBeginUser;
227:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-universal", &run, NULL));
228:   if (!run) PetscFunctionReturn(PETSC_SUCCESS);

230:   char      filename[PETSC_MAX_PATH_LEN];
231:   PetscBool flg;

233:   PetscCall(PetscOptionsGetString(NULL, NULL, "-filename", filename, sizeof(filename), &flg));
234:   if (flg) {
235:     PetscCall(DMPlexCreateFromFile(comm, filename, "ex11_plex", PETSC_TRUE, &dm1));
236:   } else {
237:     PetscCall(DMCreate(comm, &dm1));
238:     PetscCall(DMSetType(dm1, DMPLEX));
239:     PetscCall(DMSetFromOptions(dm1));
240:   }
241:   PetscCall(DMHasLabel(dm1, "marker", &notFile));
242:   if (notFile) {
243:     PetscCall(DMCreateLabel(dm1, "Boundary Faces"));
244:     PetscCall(DMGetLabel(dm1, "Boundary Faces", &bd1));
245:     PetscCall(DMPlexMarkBoundaryFaces(dm1, 13, bd1));
246:     PetscCall(DMCreateLabel(dm1, "Boundary"));
247:     PetscCall(DMGetLabel(dm1, "Boundary", &bd2));
248:     PetscCall(DMPlexMarkBoundaryFaces(dm1, 121, bd2));
249:     PetscCall(DMPlexLabelComplete(dm1, bd2));
250:   }
251:   PetscCall(PetscObjectSetName((PetscObject)dm1, "First Mesh"));
252:   PetscCall(DMViewFromOptions(dm1, NULL, "-dm_view"));

254:   PetscCall(DMUniversalLabelCreate(dm1, &universal));
255:   PetscCall(DMUniversalLabelGetLabel(universal, &ulabel));
256:   PetscCall(PetscObjectViewFromOptions((PetscObject)ulabel, NULL, "-universal_view"));

258:   if (!notFile) {
259:     PetscInt Nl, l;

261:     PetscCall(DMClone(dm1, &dm2));
262:     PetscCall(DMGetNumLabels(dm2, &Nl));
263:     for (l = Nl - 1; l >= 0; --l) {
264:       PetscBool   isdepth, iscelltype;
265:       const char *name;

267:       PetscCall(DMGetLabelName(dm2, l, &name));
268:       PetscCall(PetscStrncmp(name, "depth", 6, &isdepth));
269:       PetscCall(PetscStrncmp(name, "celltype", 9, &iscelltype));
270:       if (!isdepth && !iscelltype) PetscCall(DMRemoveLabel(dm2, name, NULL));
271:     }
272:   } else {
273:     PetscCall(DMCreate(comm, &dm2));
274:     PetscCall(DMSetType(dm2, DMPLEX));
275:     PetscCall(DMSetFromOptions(dm2));
276:   }
277:   PetscCall(PetscObjectSetName((PetscObject)dm2, "Second Mesh"));
278:   PetscCall(DMUniversalLabelCreateLabels(universal, PETSC_TRUE, dm2));
279:   PetscCall(DMPlexGetChart(dm2, &pStart, &pEnd));
280:   for (p = pStart; p < pEnd; ++p) {
281:     PetscInt val;

283:     PetscCall(DMLabelGetValue(ulabel, p, &val));
284:     if (val < 0) continue;
285:     PetscCall(DMUniversalLabelSetLabelValue(universal, dm2, PETSC_TRUE, p, val));
286:   }
287:   PetscCall(DMViewFromOptions(dm2, NULL, "-dm_view"));

289:   PetscCall(DMUniversalLabelDestroy(&universal));
290:   PetscCall(DMDestroy(&dm1));
291:   PetscCall(DMDestroy(&dm2));
292:   PetscFunctionReturn(PETSC_SUCCESS);
293: }

295: int main(int argc, char **argv)
296: {
297:   PetscFunctionBeginUser;
298:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
299:   /*PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));*/
300:   PetscCall(TestInsertion());
301:   PetscCall(TestEmptyStrata(PETSC_COMM_WORLD));
302:   PetscCall(TestDistribution(PETSC_COMM_WORLD));
303:   PetscCall(TestUniversalLabel(PETSC_COMM_WORLD));
304:   PetscCall(PetscFinalize());
305:   return 0;
306: }

308: /*TEST

310:   test:
311:     suffix: 0
312:     requires: triangle
313:   test:
314:     suffix: 1
315:     requires: triangle
316:     nsize: 2
317:     args: -petscpartitioner_type simple

319:   testset:
320:     suffix: gmsh
321:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square.msh -petscpartitioner_type simple
322:     test:
323:       suffix: 1
324:       nsize: 1
325:     test:
326:       suffix: 2
327:       nsize: 2

329:   testset:
330:     suffix: exodusii
331:     requires: exodusii
332:     args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/2Dgrd.exo -petscpartitioner_type simple
333:     test:
334:       suffix: 1
335:       nsize: 1
336:     test:
337:       suffix: 2
338:       nsize: 2

340:   test:
341:     suffix: univ
342:     requires: triangle
343:     args: -universal -dm_view -universal_view

345:   test:
346:     # Note that the labels differ because we have multiply-marked some points during EGADS creation
347:     suffix: univ_egads_sphere
348:     requires: egads
349:     args: -universal -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/unit_sphere.egadslite -dm_view -universal_view

351:   test:
352:     # Note that the labels differ because we have multiply-marked some points during EGADS creation
353:     suffix: univ_egads_ball
354:     requires: egads ctetgen
355:     args: -universal -dm_plex_boundary_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/unit_sphere.egadslite -dm_view -universal_view

357: TEST*/