xref: /petsc/src/dm/impls/plex/tests/ex49.c (revision 19620d22b41d3b8f978e0856cb1f132f4c411490) 
1 static char help[] = "Tests dof numberings for external integrators such as LibCEED.\n\n";
2 
3 #include <petscdmplex.h>
4 #include <petscds.h>
5 
6 typedef struct {
7   PetscBool useFE;
8   PetscInt  check_face;
9   PetscBool closure_tensor;
10 } AppCtx;
11 
ProcessOptions(MPI_Comm comm,AppCtx * options)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 }
25 
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)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 }
36 
SetupDiscretization(DM dm,AppCtx * user)37 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
38 {
39   DM       cdm = dm;
40   PetscInt dim;
41 
42   PetscFunctionBeginUser;
43   PetscCall(DMGetDimension(dm, &dim));
44   if (user->useFE) {
45     PetscFE        fe;
46     DMPolytopeType ct;
47     PetscInt       cStart;
48 
49     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
50     PetscCall(DMPlexGetCellType(dm, cStart, &ct));
51     PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, PETSC_DETERMINE, &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;
58 
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 }
76 
CheckOffsets(DM dm,AppCtx * user,const char * domain_name,PetscInt label_value,PetscInt height)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;
84 
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];
96 
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;
102 
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;
112 
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"));
131 
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;
140 
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];
164 
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;
191 
192       PetscCall(DMClone(dm, &clonedm));
193       { // Force recreation of local coordinate vector
194         Vec X_global;
195 
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);
204 
205       PetscCall(VecGetArrayRead(X, &x));
206       PetscCall(VecGetArrayRead(cloneX, &clonex));
207 
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       }
211 
212       PetscCall(VecRestoreArrayRead(X, &x));
213       PetscCall(VecRestoreArrayRead(cloneX, &clonex));
214       PetscCall(DMDestroy(&clonedm));
215     }
216   }
217   PetscFunctionReturn(PETSC_SUCCESS);
218 }
219 
main(int argc,char ** argv)220 int main(int argc, char **argv)
221 {
222   DM        dm;
223   AppCtx    user;
224   PetscInt  depth;
225   PetscBool flg;
226 
227   PetscFunctionBeginUser;
228   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
229   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
230   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
231   PetscCall(SetupDiscretization(dm, &user));
232   PetscCall(CheckOffsets(dm, &user, NULL, 0, 0));
233   PetscCall(DMPlexGetDepth(dm, &depth));
234   if (depth > 1) PetscCall(CheckOffsets(dm, &user, "Face Sets", user.check_face, 1));
235 
236   PetscCall(PetscOptionsHasName(NULL, NULL, "-view_mat", &flg));
237   if (flg) {
238     Mat A;
239     PetscCall(DMCreateMatrix(dm, &A));
240     PetscCall(MatViewFromOptions(A, NULL, "-view_mat"));
241     PetscCall(MatDestroy(&A));
242   }
243   PetscCall(DMDestroy(&dm));
244   PetscCall(PetscFinalize());
245   return 0;
246 }
247 
248 /*TEST
249 
250   test:
251     suffix: 0
252     requires: triangle
253     args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view
254 
255   test:
256     suffix: 1
257     args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \
258           -dm_view -offsets_view
259 
260   test:
261     suffix: cg_2d
262     args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \
263           -dm_view -offsets_view
264 
265   test:
266     suffix: 1d_sfc
267     args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_shape zbox -dm_plex_box_faces 3 1 -dm_view -coord_ltog_view
268 
269   test:
270     suffix: 1d_sfc_periodic
271     args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_shape zbox -dm_plex_box_faces 3 1 -dm_view -coord_ltog_view -petscspace_degree 1 -view_mat -dm_plex_box_bd periodic
272 
273   test:
274     suffix: 2d_sfc
275     nsize: 2
276     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
277 
278   test:
279     suffix: 2d_sfc_periodic
280     nsize: 2
281     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 -view_mat
282 
283   test:
284     suffix: 2d_sfc_periodic_mat
285     args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 1,1 -dm_distribute 0 -petscspace_degree 2 -dm_plex_box_bd periodic,periodic -dm_view ::ascii_info_detail -view_mat
286 
287   testset:
288     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
289     nsize: 2
290     test:
291       suffix: 2d_sfc_periodic_stranded
292       args: -dm_distribute 0 -dm_plex_box_bd none,periodic
293     test:
294       suffix: 2d_sfc_periodic_stranded_dist
295       args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_bd none,periodic
296     test:
297       suffix: 2d_sfc_biperiodic_stranded
298       args: -dm_distribute 0 -dm_plex_box_bd periodic,periodic
299     test:
300       suffix: 2d_sfc_biperiodic_stranded_dist
301       args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_bd periodic,periodic
302     test:
303       suffix: 2d_sfc_biperiodic_stranded_dist_box_label
304       args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_label_bd periodic,periodic -dm_plex_box_label
305 
306   test:
307     suffix: fv_0
308     requires: triangle
309     args: -dm_refine 1 -use_fe 0 -dm_view -offsets_view
310 
311 TEST*/
312