Actual source code: ex49.c

  1: static char help[] = "Tests dof numberings for external integrators such as LibCEED.\n\n";

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

  6: typedef struct {
  7:   PetscBool useFE;
  8:   PetscInt  check_face;
  9:   PetscBool closure_tensor;
 10: } AppCtx;

 12: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 13: {
 14:   PetscFunctionBeginUser;
 15:   options->useFE          = PETSC_TRUE;
 16:   options->check_face     = 1;
 17:   options->closure_tensor = PETSC_FALSE;
 18:   PetscOptionsBegin(comm, "", "Dof Ordering Options", "DMPLEX");
 19:   PetscCall(PetscOptionsBool("-use_fe", "Use FE or FV discretization", "ex49.c", options->useFE, &options->useFE, NULL));
 20:   PetscCall(PetscOptionsInt("-check_face", "Face set to report on", "ex49.c", options->check_face, &options->check_face, NULL));
 21:   PetscCall(PetscOptionsBool("-closure_tensor", "Use DMPlexSetClosurePermutationTensor()", "ex49.c", options->closure_tensor, &options->closure_tensor, NULL));
 22:   PetscOptionsEnd();
 23:   PetscFunctionReturn(PETSC_SUCCESS);
 24: }

 26: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 27: {
 28:   PetscFunctionBeginUser;
 29:   PetscCall(DMCreate(comm, dm));
 30:   PetscCall(DMSetType(*dm, DMPLEX));
 31:   PetscCall(DMSetFromOptions(*dm));
 32:   PetscCall(DMSetApplicationContext(*dm, user));
 33:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 34:   PetscFunctionReturn(PETSC_SUCCESS);
 35: }

 37: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
 38: {
 39:   DM       cdm = dm;
 40:   PetscInt dim;

 42:   PetscFunctionBeginUser;
 43:   PetscCall(DMGetDimension(dm, &dim));
 44:   if (user->useFE) {
 45:     PetscFE        fe;
 46:     DMPolytopeType ct;
 47:     PetscInt       cStart;

 49:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
 50:     PetscCall(DMPlexGetCellType(dm, cStart, &ct));
 51:     PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, -1, &fe));
 52:     PetscCall(PetscObjectSetName((PetscObject)fe, "scalar"));
 53:     PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
 54:     PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe));
 55:     PetscCall(PetscFEDestroy(&fe));
 56:   } else {
 57:     PetscFV fv;

 59:     PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv));
 60:     PetscCall(PetscFVSetType(fv, PETSCFVLEASTSQUARES));
 61:     PetscCall(PetscFVSetNumComponents(fv, dim));
 62:     PetscCall(PetscFVSetSpatialDimension(fv, dim));
 63:     PetscCall(PetscFVSetFromOptions(fv));
 64:     PetscCall(PetscFVSetUp(fv));
 65:     PetscCall(PetscObjectSetName((PetscObject)fv, "vector"));
 66:     PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fv));
 67:     PetscCall(PetscFVDestroy(&fv));
 68:   }
 69:   PetscCall(DMCreateDS(dm));
 70:   while (cdm) {
 71:     PetscCall(DMCopyDisc(dm, cdm));
 72:     PetscCall(DMGetCoarseDM(cdm, &cdm));
 73:   }
 74:   PetscFunctionReturn(PETSC_SUCCESS);
 75: }

 77: static PetscErrorCode CheckOffsets(DM dm, AppCtx *user, const char *domain_name, PetscInt label_value, PetscInt height)
 78: {
 79:   const char            *height_name[] = {"cells", "faces"};
 80:   DMLabel                domain_label  = NULL;
 81:   DM                     cdm;
 82:   PetscInt               Nf, f;
 83:   ISLocalToGlobalMapping ltog;

 85:   PetscFunctionBeginUser;
 86:   if (domain_name) PetscCall(DMGetLabel(dm, domain_name, &domain_label));
 87:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "## %s: '%s' {%" PetscInt_FMT "}%s\n", height_name[height], domain_name ? domain_name : "default", label_value, domain_name && !domain_label ? " (null label)" : ""));
 88:   if (domain_name && !domain_label) PetscFunctionReturn(PETSC_SUCCESS);
 89:   if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
 90:   // Offsets for cell closures
 91:   PetscCall(DMGetNumFields(dm, &Nf));
 92:   for (f = 0; f < Nf; ++f) {
 93:     PetscObject  obj;
 94:     PetscClassId id;
 95:     char         name[PETSC_MAX_PATH_LEN];

 97:     PetscCall(DMGetField(dm, f, NULL, &obj));
 98:     PetscCall(PetscObjectGetClassId(obj, &id));
 99:     if (id == PETSCFE_CLASSID) {
100:       IS        offIS;
101:       PetscInt *offsets, Ncell, Ncl, Nc, n;

103:       PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, f, &Ncell, &Ncl, &Nc, &n, &offsets));
104:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Ncell * Ncl, offsets, PETSC_OWN_POINTER, &offIS));
105:       PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f));
106:       PetscCall(PetscObjectSetName((PetscObject)offIS, name));
107:       PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view"));
108:       PetscCall(ISDestroy(&offIS));
109:     } else if (id == PETSCFV_CLASSID) {
110:       IS        offIS;
111:       PetscInt *offsets, *offsetsNeg, *offsetsPos, Nface, Nc, n, i = 0;

113:       PetscCall(DMPlexGetLocalOffsetsSupport(dm, domain_label, label_value, &Nface, &Nc, &n, &offsetsNeg, &offsetsPos));
114:       PetscCall(PetscMalloc1(Nface * Nc * 2, &offsets));
115:       for (PetscInt f = 0; f < Nface; ++f) {
116:         for (PetscInt c = 0; c < Nc; ++c) offsets[i++] = offsetsNeg[f] + c;
117:         for (PetscInt c = 0; c < Nc; ++c) offsets[i++] = offsetsPos[f] + c;
118:       }
119:       PetscCheck(i == Nface * Nc * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total offsets %" PetscInt_FMT " != %" PetscInt_FMT, i, Nface * Nc * 2);
120:       PetscCall(PetscFree(offsetsNeg));
121:       PetscCall(PetscFree(offsetsPos));
122:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nface * Nc * 2, offsets, PETSC_OWN_POINTER, &offIS));
123:       PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f));
124:       PetscCall(PetscObjectSetName((PetscObject)offIS, name));
125:       PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view"));
126:       PetscCall(ISDestroy(&offIS));
127:     } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unrecognized type for DM field %" PetscInt_FMT, f);
128:   }
129:   PetscCall(DMGetLocalToGlobalMapping(dm, &ltog));
130:   PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-ltog_view"));

132:   // Offsets for coordinates
133:   {
134:     Vec                X;
135:     PetscSection       s;
136:     const PetscScalar *x;
137:     const char        *cname;
138:     PetscInt           cdim, *offsets, Ncell, Ncl, Nc, n;
139:     PetscBool          isDG = PETSC_FALSE;

141:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
142:     if (!cdm) {
143:       PetscCall(DMGetCoordinateDM(dm, &cdm));
144:       cname = "Coordinates";
145:       PetscCall(DMGetCoordinatesLocal(dm, &X));
146:     } else {
147:       isDG  = PETSC_TRUE;
148:       cname = "DG Coordinates";
149:       PetscCall(DMGetCellCoordinatesLocal(dm, &X));
150:     }
151:     if (isDG && height) PetscFunctionReturn(PETSC_SUCCESS);
152:     if (domain_name) PetscCall(DMGetLabel(cdm, domain_name, &domain_label));
153:     if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
154:     PetscCall(DMPlexGetLocalOffsets(cdm, domain_label, label_value, height, 0, &Ncell, &Ncl, &Nc, &n, &offsets));
155:     PetscCall(DMGetCoordinateDim(dm, &cdim));
156:     PetscCheck(Nc == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometric dimension %" PetscInt_FMT " should be %" PetscInt_FMT, Nc, cdim);
157:     PetscCall(DMGetLocalSection(cdm, &s));
158:     PetscCall(VecGetArrayRead(X, &x));
159:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s by element in %s order\n", cname, user->closure_tensor ? "tensor" : "bfs"));
160:     for (PetscInt c = 0; c < Ncell; ++c) {
161:       for (PetscInt v = 0; v < Ncl; ++v) {
162:         PetscInt           off = offsets[c * Ncl + v], dgdof;
163:         const PetscScalar *vx  = &x[off];

165:         if (isDG) {
166:           PetscCall(PetscSectionGetDof(s, c, &dgdof));
167:           PetscCheck(Ncl * Nc == dgdof, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Offset size %" PetscInt_FMT " should be %" PetscInt_FMT, Ncl * Nc, dgdof);
168:         }
169:         switch (cdim) {
170:         case 1:
171:           PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0])));
172:           break;
173:         case 2:
174:           PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f, % 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]), (double)PetscRealPart(vx[1])));
175:           break;
176:         case 3:
177:           PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f, % 4.2f, % 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]), (double)PetscRealPart(vx[1]), (double)PetscRealPart(vx[2])));
178:         }
179:       }
180:     }
181:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
182:     PetscCall(VecRestoreArrayRead(X, &x));
183:     PetscCall(PetscFree(offsets));
184:     PetscCall(DMGetLocalToGlobalMapping(cdm, &ltog));
185:     PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-coord_ltog_view"));
186:     {
187:       DM                 clonedm;
188:       Vec                cloneX, X;
189:       PetscInt           clone_num_x, num_x;
190:       const PetscScalar *clonex, *x;

192:       PetscCall(DMClone(dm, &clonedm));
193:       { // Force recreation of local coordinate vector
194:         Vec X_global;

196:         PetscCall(DMGetCoordinates(dm, &X_global));
197:         PetscCall(DMSetCoordinates(clonedm, X_global));
198:       }
199:       PetscCall(DMGetCoordinatesLocal(dm, &X));
200:       PetscCall(DMGetCoordinatesLocal(clonedm, &cloneX));
201:       PetscCall(VecGetLocalSize(X, &num_x));
202:       PetscCall(VecGetLocalSize(cloneX, &clone_num_x));
203:       PetscCheck(num_x == clone_num_x, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "Cloned DM coordinate size (%" PetscInt_FMT ") different from original DM coordinate size (%" PetscInt_FMT ")", clone_num_x, num_x);

205:       PetscCall(VecGetArrayRead(X, &x));
206:       PetscCall(VecGetArrayRead(cloneX, &clonex));

208:       for (PetscInt i = 0; i < num_x; i++) {
209:         PetscCheck(PetscIsCloseAtTolScalar(x[i], clonex[i], 1e-13, 1e-13), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Original coordinate (%4.2f) and cloned coordinate (%4.2f) are different", (double)PetscRealPart(x[i]), (double)PetscRealPart(clonex[i]));
210:       }

212:       PetscCall(VecRestoreArrayRead(X, &x));
213:       PetscCall(VecRestoreArrayRead(cloneX, &clonex));
214:       PetscCall(DMDestroy(&clonedm));
215:     }
216:   }
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: int main(int argc, char **argv)
221: {
222:   DM       dm;
223:   AppCtx   user;
224:   PetscInt depth;

226:   PetscFunctionBeginUser;
227:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
228:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
229:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
230:   PetscCall(SetupDiscretization(dm, &user));
231:   PetscCall(CheckOffsets(dm, &user, NULL, 0, 0));
232:   PetscCall(DMPlexGetDepth(dm, &depth));
233:   if (depth > 1) PetscCall(CheckOffsets(dm, &user, "Face Sets", user.check_face, 1));
234:   PetscCall(DMDestroy(&dm));
235:   PetscCall(PetscFinalize());
236:   return 0;
237: }

239: /*TEST

241:   test:
242:     suffix: 0
243:     requires: triangle
244:     args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view

246:   test:
247:     suffix: 1
248:     args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \
249:           -dm_view -offsets_view

251:   test:
252:     suffix: cg_2d
253:     args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \
254:           -dm_view -offsets_view

256:   test:
257:     suffix: 1d_sfc
258:     args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_shape zbox -dm_plex_box_faces 3 1 -dm_view -coord_ltog_view

260:   test:
261:     suffix: 2d_sfc
262:     nsize: 2
263:     args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 4,3 -dm_distribute 0 -petscspace_degree 1 -dm_view

265:   test:
266:     suffix: 2d_sfc_periodic
267:     nsize: 2
268:     args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 4,3 -dm_distribute 0 -petscspace_degree 1 -dm_plex_box_bd periodic,none -dm_view ::ascii_info_detail

270:   testset:
271:     args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 3,2 -petscspace_degree 1 -dm_view ::ascii_info_detail -closure_tensor
272:     nsize: 2
273:     test:
274:       suffix: 2d_sfc_periodic_stranded
275:       args: -dm_distribute 0 -dm_plex_box_bd none,periodic
276:     test:
277:       suffix: 2d_sfc_periodic_stranded_dist
278:       args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_bd none,periodic
279:     test:
280:       suffix: 2d_sfc_biperiodic_stranded
281:       args: -dm_distribute 0 -dm_plex_box_bd periodic,periodic
282:     test:
283:       suffix: 2d_sfc_biperiodic_stranded_dist
284:       args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_bd periodic,periodic

286:   test:
287:     suffix: fv_0
288:     requires: triangle
289:     args: -dm_refine 1 -use_fe 0 -dm_view -offsets_view

291: TEST*/