Actual source code: ex3.c
1: static const char help[] = "Tests for determining whether a new finite element works";
3: /*
4: Use -interpolation_view and -l2_projection_view to look at the interpolants.
5: */
7: #include <petscdmplex.h>
8: #include <petscfe.h>
9: #include <petscds.h>
10: #include <petscsnes.h>
12: static void constant(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
13: {
14: const PetscInt Nc = uOff[1] - uOff[0];
15: for (PetscInt c = 0; c < Nc; ++c) f0[c] += 5.;
16: }
18: static void linear(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
19: {
20: const PetscInt Nc = uOff[1] - uOff[0];
21: for (PetscInt c = 0; c < Nc; ++c) f0[c] += 5. * x[c];
22: }
24: static void quadratic(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
25: {
26: const PetscInt Nc = uOff[1] - uOff[0];
27: for (PetscInt c = 0; c < Nc; ++c) f0[c] += 5. * x[c] * x[c];
28: }
30: static void trig(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
31: {
32: const PetscInt Nc = uOff[1] - uOff[0];
33: for (PetscInt c = 0; c < Nc; ++c) f0[c] += PetscCosReal(2. * PETSC_PI * x[c]);
34: }
36: /*
37: The prime basis for the Wheeler-Yotov-Xue prism.
38: */
39: static void prime(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
40: {
41: PetscReal x = X[0], y = X[1], z = X[2], b = 1 + x + y + z;
42: f0[0] += b + 2.0 * x * z + 2.0 * y * z + x * y + x * x;
43: f0[1] += b + 2.0 * x * z + 2.0 * y * z + x * y + y * y;
44: f0[2] += b - 3.0 * x * z - 3.0 * y * z - 2.0 * z * z;
45: }
47: static const char *names[] = {"constant", "linear", "quadratic", "trig", "prime"};
48: static PetscPointFunc functions[] = {constant, linear, quadratic, trig, prime};
50: typedef struct {
51: PetscPointFunc exactSol;
52: PetscReal shear, flatten;
53: } AppCtx;
55: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
56: {
57: char name[PETSC_MAX_PATH_LEN] = "constant";
58: PetscInt Nfunc = PETSC_STATIC_ARRAY_LENGTH(names), i;
61: options->exactSol = NULL;
62: options->shear = 0.;
63: options->flatten = 1.;
65: PetscOptionsBegin(comm, "", "FE Test Options", "PETSCFE");
66: PetscOptionsString("-func", "Function to project into space", "", name, name, PETSC_MAX_PATH_LEN, NULL);
67: PetscOptionsReal("-shear", "Factor by which to shear along the x-direction", "", options->shear, &(options->shear), NULL);
68: PetscOptionsReal("-flatten", "Factor by which to flatten", "", options->flatten, &(options->flatten), NULL);
69: PetscOptionsEnd();
71: for (i = 0; i < Nfunc; ++i) {
72: PetscBool flg;
74: PetscStrcmp(name, names[i], &flg);
75: if (flg) {
76: options->exactSol = functions[i];
77: break;
78: }
79: }
81: return 0;
82: }
84: /* The exact solution is the negative of the f0 contribution */
85: static PetscErrorCode exactSolution(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
86: {
87: AppCtx *user = (AppCtx *)ctx;
88: PetscInt uOff[2] = {0, Nc};
90: user->exactSol(dim, 1, 0, uOff, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, time, x, 0, NULL, u);
91: for (PetscInt c = 0; c < Nc; ++c) u[c] *= -1.;
92: return 0;
93: }
95: static void f0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
96: {
97: const PetscInt Nc = uOff[1] - uOff[0];
98: for (PetscInt c = 0; c < Nc; ++c) f0[c] += u[c];
99: }
101: static void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
102: {
103: const PetscInt Nc = uOff[1] - uOff[0];
104: for (PetscInt c = 0; c < Nc; ++c) g0[c * Nc + c] = 1.0;
105: }
107: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
108: {
109: PetscDS ds;
110: PetscWeakForm wf;
113: DMGetDS(dm, &ds);
114: PetscDSGetWeakForm(ds, &wf);
115: PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, f0, 0, NULL);
116: PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 1, user->exactSol, 0, NULL);
117: PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, g0, 0, NULL, 0, NULL, 0, NULL);
118: PetscDSSetExactSolution(ds, 0, exactSolution, user);
119: return 0;
120: }
122: static PetscErrorCode SetupDiscretization(DM dm, const char name[], AppCtx *user)
123: {
124: DM cdm = dm;
125: PetscFE fe;
126: char prefix[PETSC_MAX_PATH_LEN];
129: if (name) PetscSNPrintf(prefix, sizeof(prefix), "%s_", name);
130: DMCreateFEDefault(dm, 1, name ? prefix : NULL, -1, &fe);
131: PetscObjectSetName((PetscObject)fe, name ? name : "Solution");
132: /* Set discretization and boundary conditions for each mesh */
133: DMSetField(dm, 0, NULL, (PetscObject)fe);
134: DMCreateDS(dm);
135: SetupProblem(dm, user);
136: while (cdm) {
137: DMCopyDisc(dm, cdm);
138: DMGetCoarseDM(cdm, &cdm);
139: }
140: PetscFEDestroy(&fe);
141: return 0;
142: }
144: /* This test tells us whether the given function is contained in the approximation space */
145: static PetscErrorCode CheckInterpolation(DM dm, AppCtx *user)
146: {
147: PetscSimplePointFunc exactSol[1];
148: void *exactCtx[1];
149: PetscDS ds;
150: Vec u;
151: PetscReal error, tol = PETSC_SMALL;
152: MPI_Comm comm;
155: PetscObjectGetComm((PetscObject)dm, &comm);
156: DMGetDS(dm, &ds);
157: DMGetGlobalVector(dm, &u);
158: PetscDSGetExactSolution(ds, 0, &exactSol[0], &exactCtx[0]);
159: DMProjectFunction(dm, 0.0, exactSol, exactCtx, INSERT_ALL_VALUES, u);
160: VecViewFromOptions(u, NULL, "-interpolation_view");
161: DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error);
162: DMRestoreGlobalVector(dm, &u);
163: if (error > tol) PetscPrintf(comm, "Interpolation tests FAIL at tolerance %g error %g\n", (double)tol, (double)error);
164: else PetscPrintf(comm, "Interpolation tests pass at tolerance %g\n", (double)tol);
165: return 0;
166: }
168: /* This test tells us whether the element is unisolvent (the mass matrix has full rank), and what rate of convergence we achieve */
169: static PetscErrorCode CheckL2Projection(DM dm, AppCtx *user)
170: {
171: PetscSimplePointFunc exactSol[1];
172: void *exactCtx[1];
173: SNES snes;
174: PetscDS ds;
175: Vec u;
176: PetscReal error, tol = PETSC_SMALL;
177: MPI_Comm comm;
180: PetscObjectGetComm((PetscObject)dm, &comm);
181: DMGetDS(dm, &ds);
182: DMGetGlobalVector(dm, &u);
183: PetscDSGetExactSolution(ds, 0, &exactSol[0], &exactCtx[0]);
184: SNESCreate(comm, &snes);
185: SNESSetDM(snes, dm);
186: VecSet(u, 0.0);
187: PetscObjectSetName((PetscObject)u, "solution");
188: DMPlexSetSNESLocalFEM(dm, user, user, user);
189: SNESSetFromOptions(snes);
190: DMSNESCheckFromOptions(snes, u);
191: SNESSolve(snes, NULL, u);
192: SNESDestroy(&snes);
193: VecViewFromOptions(u, NULL, "-l2_projection_view");
194: DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error);
195: DMRestoreGlobalVector(dm, &u);
196: if (error > tol) PetscPrintf(comm, "L2 projection tests FAIL at tolerance %g error %g\n", (double)tol, (double)error);
197: else PetscPrintf(comm, "L2 projection tests pass at tolerance %g\n", (double)tol);
198: return 0;
199: }
201: /* Distorts the mesh by shearing in the x-direction and flattening, factors provided in the options. */
202: static PetscErrorCode DistortMesh(DM dm, AppCtx *user)
203: {
204: Vec coordinates;
205: PetscScalar *ca;
206: PetscInt dE, n, i;
209: DMGetCoordinateDim(dm, &dE);
210: DMGetCoordinates(dm, &coordinates);
211: VecGetLocalSize(coordinates, &n);
212: VecGetArray(coordinates, &ca);
213: for (i = 0; i < (n / dE); ++i) {
214: ca[i * dE + 0] += user->shear * ca[i * dE + 0];
215: ca[i * dE + 1] *= user->flatten;
216: }
217: VecRestoreArray(coordinates, &ca);
218: return 0;
219: }
221: int main(int argc, char **argv)
222: {
223: DM dm;
224: AppCtx user;
225: PetscMPIInt size;
228: PetscInitialize(&argc, &argv, NULL, help);
229: MPI_Comm_size(PETSC_COMM_WORLD, &size);
231: ProcessOptions(PETSC_COMM_WORLD, &user);
232: DMCreate(PETSC_COMM_WORLD, &dm);
233: DMSetType(dm, DMPLEX);
234: DMSetFromOptions(dm);
235: DistortMesh(dm, &user);
236: DMViewFromOptions(dm, NULL, "-dm_view");
237: SetupDiscretization(dm, NULL, &user);
239: CheckInterpolation(dm, &user);
240: CheckL2Projection(dm, &user);
242: DMDestroy(&dm);
243: PetscFinalize();
244: return 0;
245: }
247: /*TEST
249: testset:
250: args: -dm_plex_reference_cell_domain -dm_plex_cell triangle -petscspace_degree 1\
251: -snes_error_if_not_converged -ksp_error_if_not_converged -pc_type lu
253: test:
254: suffix: p1_0
255: args: -func {{constant linear}}
257: # Using -dm_refine 2 -convest_num_refine 4 gives convergence rate 2.0
258: test:
259: suffix: p1_1
260: args: -func {{quadratic trig}} \
261: -snes_convergence_estimate -convest_num_refine 2
263: testset:
264: requires: !complex double
265: args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \
266: -petscspace_type sum \
267: -petscspace_variables 3 \
268: -petscspace_components 3 \
269: -petscspace_sum_spaces 2 \
270: -petscspace_sum_concatenate false \
271: -sumcomp_0_petscspace_variables 3 \
272: -sumcomp_0_petscspace_components 3 \
273: -sumcomp_0_petscspace_degree 1 \
274: -sumcomp_1_petscspace_variables 3 \
275: -sumcomp_1_petscspace_components 3 \
276: -sumcomp_1_petscspace_type wxy \
277: -petscdualspace_form_degree 0 \
278: -petscdualspace_order 1 \
279: -petscdualspace_components 3 \
280: -snes_error_if_not_converged -ksp_error_if_not_converged -pc_type lu
282: test:
283: suffix: wxy_0
284: args: -func constant
286: test:
287: suffix: wxy_1
288: args: -func linear
290: test:
291: suffix: wxy_2
292: args: -func prime
294: test:
295: suffix: wxy_3
296: args: -func linear -shear 1 -flatten 1e-5
298: TEST*/