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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
33: }
35: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
36: {
37: PetscFunctionBeginUser;
38: options->Nr = 1;
39: PetscOptionsBegin(comm, "", "1D Refinement Options", "DMPLEX");
40: PetscCall(PetscOptionsInt("-num_refine", "Refine cycles", "ex46.c", options->Nr, &options->Nr, NULL));
41: PetscOptionsEnd();
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
46: {
47: PetscFunctionBeginUser;
48: PetscCall(DMCreate(comm, dm));
49: PetscCall(DMSetType(*dm, DMPLEX));
50: PetscCall(DMSetFromOptions(*dm));
51: PetscCall(DMSetApplicationContext(*dm, user));
52: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
57: {
58: DM cdm = dm;
59: PetscFE fe;
60: PetscSpace sp;
61: PetscInt dim, deg;
63: PetscFunctionBeginUser;
64: PetscCall(DMGetDimension(dm, &dim));
65: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, -1, &fe));
66: PetscCall(PetscObjectSetName((PetscObject)fe, "scalar"));
67: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
68: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe));
69: PetscCall(DMCreateDS(dm));
70: while (cdm) {
71: PetscCall(DMCopyDisc(dm, cdm));
72: PetscCall(DMGetCoarseDM(cdm, &cdm));
73: }
74: PetscCall(PetscFEGetBasisSpace(fe, &sp));
75: PetscCall(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: PetscCall(PetscFEDestroy(&fe));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode CheckError(DM dm, Vec u, PetscSimplePointFunc funcs[])
98: {
99: PetscReal error, tol = PETSC_SMALL;
100: MPI_Comm comm;
102: PetscFunctionBeginUser;
103: PetscCall(DMComputeL2Diff(dm, 0.0, funcs, NULL, u, &error));
104: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
105: if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL at tolerance %g error %g\n", (double)tol, (double)error));
106: else PetscCall(PetscPrintf(comm, "Function tests pass at tolerance %g\n", (double)tol));
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: int main(int argc, char **argv)
111: {
112: DM dm;
113: Vec u;
114: AppCtx user;
115: PetscInt cStart, cEnd, c, r;
117: PetscFunctionBeginUser;
118: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
119: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
120: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
121: PetscCall(SetupDiscretization(dm, &user));
122: PetscCall(DMGetGlobalVector(dm, &u));
123: PetscCall(DMProjectFunction(dm, 0.0, user.funcs, NULL, INSERT_ALL_VALUES, u));
124: PetscCall(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: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adapt));
132: PetscCall(DMLabelSetDefaultValue(adapt, DM_ADAPT_COARSEN));
133: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
134: for (c = cStart; c < cEnd; ++c) {
135: if (c % 2) PetscCall(DMLabelSetValue(adapt, c, DM_ADAPT_REFINE));
136: }
137: PetscCall(DMAdaptLabel(dm, adapt, &adm));
138: PetscCall(DMLabelDestroy(&adapt));
139: PetscCall(PetscObjectSetName((PetscObject)adm, "Adapted Mesh"));
140: PetscCall(DMViewFromOptions(adm, NULL, "-dm_view"));
142: PetscCall(DMCreateInterpolation(dm, adm, &Interp, NULL));
143: PetscCall(DMGetGlobalVector(adm, &au));
144: PetscCall(MatInterpolate(Interp, u, au));
145: PetscCall(CheckError(adm, au, user.funcs));
146: PetscCall(MatDestroy(&Interp));
147: PetscCall(DMRestoreGlobalVector(dm, &u));
148: PetscCall(DMDestroy(&dm));
149: dm = adm;
150: u = au;
151: }
152: PetscCall(DMRestoreGlobalVector(dm, &u));
153: PetscCall(DMDestroy(&dm));
154: PetscCall(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*/