Actual source code: ex46.c
1: static char help[] = "Tests 1D nested mesh refinement.\n\n";
3: #include <petscdmplex.h>
4: #include <petscds.h>
6: typedef struct {
7: PetscInt Nr; /* Number of refinements */
8: PetscSimplePointFunc funcs[2]; /* Functions to test */
9: } AppCtx;
11: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
12: {
13: u[0] = 1.;
14: return 0;
15: }
17: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
18: {
19: u[0] = 2. * x[0] + 1.;
20: return 0;
21: }
23: static PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
24: {
25: u[0] = 3. * x[0] * x[0] + 2. * x[0] + 1.;
26: return 0;
27: }
29: static PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
30: {
31: u[0] = 4. * x[0] * x[0] * x[0] + 3. * x[0] * x[0] + 2. * x[0] + 1.;
32: return 0;
33: }
35: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
36: {
38: options->Nr = 1;
39: PetscOptionsBegin(comm, "", "1D Refinement Options", "DMPLEX");
40: PetscOptionsInt("-num_refine", "Refine cycles", "ex46.c", options->Nr, &options->Nr, NULL);
41: PetscOptionsEnd();
42: return 0;
43: }
45: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
46: {
48: DMCreate(comm, dm);
49: DMSetType(*dm, DMPLEX);
50: DMSetFromOptions(*dm);
51: DMSetApplicationContext(*dm, user);
52: DMViewFromOptions(*dm, NULL, "-dm_view");
53: return 0;
54: }
56: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
57: {
58: DM cdm = dm;
59: PetscFE fe;
60: PetscSpace sp;
61: PetscInt dim, deg;
64: DMGetDimension(dm, &dim);
65: PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, -1, &fe);
66: PetscObjectSetName((PetscObject)fe, "scalar");
67: DMSetField(dm, 0, NULL, (PetscObject)fe);
68: DMSetField(dm, 1, NULL, (PetscObject)fe);
69: DMCreateDS(dm);
70: while (cdm) {
71: DMCopyDisc(dm, cdm);
72: DMGetCoarseDM(cdm, &cdm);
73: }
74: PetscFEGetBasisSpace(fe, &sp);
75: PetscSpaceGetDegree(sp, °, NULL);
76: switch (deg) {
77: case 0:
78: user->funcs[0] = constant;
79: break;
80: case 1:
81: user->funcs[0] = linear;
82: break;
83: case 2:
84: user->funcs[0] = quadratic;
85: break;
86: case 3:
87: user->funcs[0] = cubic;
88: break;
89: default:
90: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Could not determine function to test for degree %" PetscInt_FMT, deg);
91: }
92: user->funcs[1] = user->funcs[0];
93: PetscFEDestroy(&fe);
94: return 0;
95: }
97: static PetscErrorCode CheckError(DM dm, Vec u, PetscSimplePointFunc funcs[])
98: {
99: PetscReal error, tol = PETSC_SMALL;
100: MPI_Comm comm;
103: DMComputeL2Diff(dm, 0.0, funcs, NULL, u, &error);
104: PetscObjectGetComm((PetscObject)dm, &comm);
105: if (error > tol) PetscPrintf(comm, "Function tests FAIL at tolerance %g error %g\n", (double)tol, (double)error);
106: else PetscPrintf(comm, "Function tests pass at tolerance %g\n", (double)tol);
107: return 0;
108: }
110: int main(int argc, char **argv)
111: {
112: DM dm;
113: Vec u;
114: AppCtx user;
115: PetscInt cStart, cEnd, c, r;
118: PetscInitialize(&argc, &argv, NULL, help);
119: ProcessOptions(PETSC_COMM_WORLD, &user);
120: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
121: SetupDiscretization(dm, &user);
122: DMGetGlobalVector(dm, &u);
123: DMProjectFunction(dm, 0.0, user.funcs, NULL, INSERT_ALL_VALUES, u);
124: CheckError(dm, u, user.funcs);
125: for (r = 0; r < user.Nr; ++r) {
126: DM adm;
127: DMLabel adapt;
128: Vec au;
129: Mat Interp;
131: DMLabelCreate(PETSC_COMM_SELF, "adapt", &adapt);
132: DMLabelSetDefaultValue(adapt, DM_ADAPT_COARSEN);
133: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
134: for (c = cStart; c < cEnd; ++c) {
135: if (c % 2) DMLabelSetValue(adapt, c, DM_ADAPT_REFINE);
136: }
137: DMAdaptLabel(dm, adapt, &adm);
138: DMLabelDestroy(&adapt);
139: PetscObjectSetName((PetscObject)adm, "Adapted Mesh");
140: DMViewFromOptions(adm, NULL, "-dm_view");
142: DMCreateInterpolation(dm, adm, &Interp, NULL);
143: DMGetGlobalVector(adm, &au);
144: MatInterpolate(Interp, u, au);
145: CheckError(adm, au, user.funcs);
146: MatDestroy(&Interp);
147: DMRestoreGlobalVector(dm, &u);
148: DMDestroy(&dm);
149: dm = adm;
150: u = au;
151: }
152: DMRestoreGlobalVector(dm, &u);
153: DMDestroy(&dm);
154: PetscFinalize();
155: return 0;
156: }
158: /*TEST
160: test:
161: suffix: 0
162: args: -num_refine 4 -petscspace_degree 3 \
163: -dm_plex_dim 1 -dm_plex_box_faces 5 -dm_plex_transform_type refine_1d -dm_plex_hash_location -dm_view
165: TEST*/