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

 19:     DMLabelGetValue(label, i, &val);
 21:   }
 22:   /* Test stratum */
 23:   for (v = 0; v < 5; ++v) {
 24:     IS              stratum;
 25:     const PetscInt *points;
 26:     PetscInt        n;

 28:     DMLabelGetStratumIS(label, values[v], &stratum);
 30:     ISGetIndices(stratum, &points);
 31:     ISGetLocalSize(stratum, &n);
 33:     ISRestoreIndices(stratum, &points);
 34:     ISDestroy(&stratum);
 35:   }
 36:   /* Test get in array mode */
 37:   for (i = 0; i < N; ++i) {
 38:     PetscInt val;

 40:     DMLabelGetValue(label, i, &val);
 42:   }
 43:   /* Test Duplicate */
 44:   DMLabelDuplicate(label, &label2);
 45:   for (i = 0; i < N; ++i) {
 46:     PetscInt val;

 48:     DMLabelGetValue(label2, i, &val);
 50:   }
 51:   DMLabelDestroy(&label2);
 52:   DMLabelDestroy(&label);
 53:   return 0;
 54: }

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

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

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

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

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

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

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

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

224:   PetscOptionsGetBool(NULL, NULL, "-universal", &run, NULL);
225:   if (!run) return 0;

227:   char      filename[PETSC_MAX_PATH_LEN];
228:   PetscBool flg;

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

251:   DMUniversalLabelCreate(dm1, &universal);
252:   DMUniversalLabelGetLabel(universal, &ulabel);
253:   PetscObjectViewFromOptions((PetscObject)ulabel, NULL, "-universal_view");

255:   if (!notFile) {
256:     PetscInt Nl, l;

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

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

280:     DMLabelGetValue(ulabel, p, &val);
281:     if (val < 0) continue;
282:     DMUniversalLabelSetLabelValue(universal, dm2, PETSC_TRUE, p, val);
283:   }
284:   DMViewFromOptions(dm2, NULL, "-dm_view");

286:   DMUniversalLabelDestroy(&universal);
287:   DMDestroy(&dm1);
288:   DMDestroy(&dm2);
289:   return 0;
290: }

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

305: /*TEST

307:   test:
308:     suffix: 0
309:     requires: triangle
310:   test:
311:     suffix: 1
312:     requires: triangle
313:     nsize: 2
314:     args: -petscpartitioner_type simple

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

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

337:   test:
338:     suffix: univ
339:     requires: triangle
340:     args: -universal -dm_view -universal_view

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

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

354: TEST*/