1 static char help[] = "Tests 1D nested mesh refinement.\n\n";
2
3 #include <petscdmplex.h>
4 #include <petscds.h>
5
6 typedef struct {
7 PetscInt Nr; /* Number of refinements */
8 PetscSimplePointFn *funcs[2]; /* Functions to test */
9 } AppCtx;
10
constant(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)11 static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
12 {
13 u[0] = 1.;
14 return PETSC_SUCCESS;
15 }
16
linear(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)17 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
18 {
19 u[0] = 2. * x[0] + 1.;
20 return PETSC_SUCCESS;
21 }
22
quadratic(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)23 static PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
24 {
25 u[0] = 3. * x[0] * x[0] + 2. * x[0] + 1.;
26 return PETSC_SUCCESS;
27 }
28
cubic(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)29 static PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx 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 }
34
ProcessOptions(MPI_Comm comm,AppCtx * options)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 }
44
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)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 }
55
SetupDiscretization(DM dm,AppCtx * user)56 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
57 {
58 DM cdm = dm;
59 PetscFE fe;
60 PetscSpace sp;
61 PetscInt dim, deg;
62
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 }
96
CheckError(DM dm,Vec u,PetscSimplePointFn * funcs[])97 static PetscErrorCode CheckError(DM dm, Vec u, PetscSimplePointFn *funcs[])
98 {
99 PetscReal error, tol = PETSC_SMALL;
100 MPI_Comm comm;
101
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 }
109
main(int argc,char ** argv)110 int main(int argc, char **argv)
111 {
112 DM dm;
113 Vec u;
114 AppCtx user;
115 PetscInt cStart, cEnd, c, r;
116
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;
130
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"));
141
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 }
157
158 /*TEST
159
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
164
165 TEST*/
166