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: PetscInt dummy;
8: } AppCtx;
10: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
11: {
13: options->dummy = 1;
14: PetscOptionsBegin(comm, "", "Dof Ordering Options", "DMPLEX");
15: //PetscOptionsInt("-num_refine", "Refine cycles", "ex46.c", options->Nr, &options->Nr, NULL);
16: PetscOptionsEnd();
17: return 0;
18: }
20: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
21: {
23: DMCreate(comm, dm);
24: DMSetType(*dm, DMPLEX);
25: DMSetFromOptions(*dm);
26: DMSetApplicationContext(*dm, user);
27: DMViewFromOptions(*dm, NULL, "-dm_view");
28: return 0;
29: }
31: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
32: {
33: DM cdm = dm;
34: PetscFE fe;
35: PetscInt dim;
38: DMGetDimension(dm, &dim);
39: PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, -1, &fe);
40: PetscObjectSetName((PetscObject)fe, "scalar");
41: DMSetField(dm, 0, NULL, (PetscObject)fe);
42: DMSetField(dm, 1, NULL, (PetscObject)fe);
43: PetscFEDestroy(&fe);
44: DMCreateDS(dm);
45: while (cdm) {
46: DMCopyDisc(dm, cdm);
47: DMGetCoarseDM(cdm, &cdm);
48: }
49: return 0;
50: }
52: static PetscErrorCode CheckOffsets(DM dm, const char *domain_name, PetscInt label_value, PetscInt height)
53: {
54: const char *height_name[] = {"cells", "faces"};
55: DMLabel domain_label = NULL;
56: DM cdm;
57: IS offIS;
58: PetscInt *offsets, Ncell, Ncl, Nc, n;
59: PetscInt Nf, f;
62: if (domain_name) DMGetLabel(dm, domain_name, &domain_label);
63: PetscPrintf(PETSC_COMM_WORLD, "## %s: '%s' {%" PetscInt_FMT "}\n", height_name[height], domain_name ? domain_name : "default", label_value);
64: // Offsets for cell closures
65: DMGetNumFields(dm, &Nf);
66: for (f = 0; f < Nf; ++f) {
67: char name[PETSC_MAX_PATH_LEN];
69: DMPlexGetLocalOffsets(dm, domain_label, label_value, height, f, &Ncell, &Ncl, &Nc, &n, &offsets);
70: ISCreateGeneral(PETSC_COMM_SELF, Ncell * Ncl, offsets, PETSC_OWN_POINTER, &offIS);
71: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f);
72: PetscObjectSetName((PetscObject)offIS, name);
73: ISViewFromOptions(offIS, NULL, "-offsets_view");
74: ISDestroy(&offIS);
75: }
76: // Offsets for coordinates
77: {
78: Vec X;
79: PetscSection s;
80: const PetscScalar *x;
81: const char *cname;
82: PetscInt cdim;
83: PetscBool isDG = PETSC_FALSE;
85: DMGetCellCoordinateDM(dm, &cdm);
86: if (!cdm) {
87: DMGetCoordinateDM(dm, &cdm);
88: cname = "Coordinates";
89: DMGetCoordinatesLocal(dm, &X);
90: } else {
91: isDG = PETSC_TRUE;
92: cname = "DG Coordinates";
93: DMGetCellCoordinatesLocal(dm, &X);
94: }
95: if (isDG && height) return 0;
96: if (domain_name) DMGetLabel(cdm, domain_name, &domain_label);
97: DMPlexGetLocalOffsets(cdm, domain_label, label_value, height, 0, &Ncell, &Ncl, &Nc, &n, &offsets);
98: DMGetCoordinateDim(dm, &cdim);
100: DMGetLocalSection(cdm, &s);
101: VecGetArrayRead(X, &x);
102: PetscPrintf(PETSC_COMM_WORLD, "%s by element\n", cname);
103: for (PetscInt c = 0; c < Ncell; ++c) {
104: for (PetscInt v = 0; v < Ncl; ++v) {
105: PetscInt off = offsets[c * Ncl + v], dgdof;
106: const PetscScalar *vx = &x[off];
108: if (isDG) {
109: PetscSectionGetDof(s, c, &dgdof);
111: }
112: switch (cdim) {
113: case 2:
114: PetscPrintf(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]));
115: break;
116: case 3:
117: PetscPrintf(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]));
118: break;
119: }
120: }
121: }
122: VecRestoreArrayRead(X, &x);
123: PetscFree(offsets);
124: }
125: return 0;
126: }
128: int main(int argc, char **argv)
129: {
130: DM dm;
131: AppCtx user;
134: PetscInitialize(&argc, &argv, NULL, help);
135: ProcessOptions(PETSC_COMM_WORLD, &user);
136: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
137: SetupDiscretization(dm, &user);
138: CheckOffsets(dm, NULL, 0, 0);
139: CheckOffsets(dm, "Face Sets", 1, 1);
140: DMDestroy(&dm);
141: PetscFinalize();
142: return 0;
143: }
145: /*TEST
147: test:
148: suffix: 0
149: requires: triangle
150: args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view
152: test:
153: suffix: 1
154: args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \
155: -dm_view -offsets_view
157: test:
158: suffix: cg_2d
159: args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \
160: -dm_view -offsets_view
161: TEST*/