Actual source code: ex39.c
1: const char help[] = "A test of H-div conforming discretizations on different cell types.\n";
3: #include <petscdmplex.h>
4: #include <petscds.h>
5: #include <petscsnes.h>
6: #include <petscconvest.h>
7: #include <petscfe.h>
8: #include <petsc/private/petscfeimpl.h>
10: /*
11: We are using the system
13: \vec{u} = \vec{\hat{u}}
14: p = \div{\vec{u}} in low degree approximation space
15: d = \div{\vec{u}} - p == 0 in higher degree approximation space
17: That is, we are using the field d to compute the error between \div{\vec{u}}
18: computed in a space 1 degree higher than p and the value of p which is
19: \div{u} computed in the low degree space. If H-div
20: elements are implemented correctly then this should be identically zero since
21: the divergence of a function in H(div) should be exactly representable in L_2
22: by definition.
23: */
24: static PetscErrorCode zero_func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
25: {
26: PetscInt c;
27: for (c = 0; c < Nc; ++c) u[c] = 0;
28: return 0;
29: }
30: /* Linear Exact Functions
31: \vec{u} = \vec{x};
32: p = dim;
33: */
34: static PetscErrorCode linear_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
35: {
36: PetscInt c;
37: for (c = 0; c < Nc; ++c) u[c] = x[c];
38: return 0;
39: }
40: static PetscErrorCode linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
41: {
42: u[0] = dim;
43: return 0;
44: }
46: /* Sinusoidal Exact Functions
47: * u_i = \sin{2*\pi*x_i}
48: * p = \Sum_{i=1}^{dim} 2*\pi*cos{2*\pi*x_i}
49: * */
51: static PetscErrorCode sinusoid_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
52: {
53: PetscInt c;
54: for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(2 * PETSC_PI * x[c]);
55: return 0;
56: }
57: static PetscErrorCode sinusoid_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
58: {
59: PetscInt d;
60: u[0] = 0;
61: for (d = 0; d < dim; ++d) u[0] += 2 * PETSC_PI * PetscCosReal(2 * PETSC_PI * x[d]);
62: return 0;
63: }
65: /* Pointwise residual for u = u*. Need one of these for each possible u* */
66: static void f0_v_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[])
67: {
68: PetscInt i;
69: PetscScalar *u_rhs;
71: PetscCalloc1(dim, &u_rhs);
72: (void)linear_u(dim, t, x, dim, u_rhs, NULL);
73: for (i = 0; i < dim; ++i) f0[i] = u[uOff[0] + i] - u_rhs[i];
74: PetscFree(u_rhs);
75: }
77: static void f0_v_sinusoid(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[])
78: {
79: PetscInt i;
80: PetscScalar *u_rhs;
82: PetscCalloc1(dim, &u_rhs);
83: (void)sinusoid_u(dim, t, x, dim, u_rhs, NULL);
84: for (i = 0; i < dim; ++i) f0[i] = u[uOff[0] + i] - u_rhs[i];
85: PetscFree(u_rhs);
86: }
88: /* Residual function for enforcing p = \div{u}. */
89: static void f0_q(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[])
90: {
91: PetscInt i;
92: PetscScalar divu;
94: divu = 0.;
95: for (i = 0; i < dim; ++i) divu += u_x[uOff_x[0] + i * dim + i];
96: f0[0] = u[uOff[1]] - divu;
97: }
99: /* Residual function for p_err = \div{u} - p. */
100: static void f0_w(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[])
101: {
102: PetscInt i;
103: PetscScalar divu;
105: divu = 0.;
106: for (i = 0; i < dim; ++i) divu += u_x[uOff_x[0] + i * dim + i];
107: f0[0] = u[uOff[2]] - u[uOff[1]] + divu;
108: }
110: /* Boundary residual for the embedding system. Need one for each form of
111: * solution. These enforce u = \hat{u} at the boundary. */
112: static void f0_bd_u_sinusoid(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
113: {
114: PetscInt d;
115: PetscScalar *u_rhs;
116: PetscCalloc1(dim, &u_rhs);
117: (void)sinusoid_u(dim, t, x, dim, u_rhs, NULL);
119: for (d = 0; d < dim; ++d) f0[d] = u_rhs[d];
120: PetscFree(u_rhs);
121: }
123: static void f0_bd_u_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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
124: {
125: PetscInt d;
126: PetscScalar *u_rhs;
127: PetscCalloc1(dim, &u_rhs);
128: (void)linear_u(dim, t, x, dim, u_rhs, NULL);
130: for (d = 0; d < dim; ++d) f0[d] = u_rhs[d];
131: PetscFree(u_rhs);
132: }
133: /* Jacobian functions. For the following, v is the test function associated with
134: * u, q the test function associated with p, and w the test function associated
135: * with d. */
136: /* <v, u> */
137: static void g0_vu(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[])
138: {
139: PetscInt c;
140: for (c = 0; c < dim; ++c) g0[c * dim + c] = 1.0;
141: }
143: /* <q, p> */
144: static void g0_qp(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[])
145: {
146: PetscInt d;
147: for (d = 0; d < dim; ++d) g0[d * dim + d] = 1.0;
148: }
150: /* -<q, \div{u}> For the embedded system. This is different from the method of
151: * manufactured solution because instead of computing <q,\div{u}> - <q,f> we
152: * need <q,p> - <q,\div{u}.*/
153: static void g1_qu(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 g1[])
154: {
155: PetscInt d;
156: for (d = 0; d < dim; ++d) g1[d * dim + d] = -1.0;
157: }
159: /* <w, p> This is only used by the embedded system. Where we need to compute
160: * <w,d> - <w,p> + <w, \div{u}>*/
161: static void g0_wp(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[])
162: {
163: PetscInt d;
164: for (d = 0; d < dim; ++d) g0[d * dim + d] = -1.0;
165: }
167: /* <w, d> */
168: static void g0_wd(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[])
169: {
170: PetscInt c;
171: for (c = 0; c < dim; ++c) g0[c * dim + c] = 1.0;
172: }
174: /* <w, \div{u}> for the embedded system. */
175: static void g1_wu(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 g1[])
176: {
177: PetscInt d;
178: for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
179: }
181: /* Enum and string array for selecting mesh perturbation options */
182: typedef enum {
183: NONE = 0,
184: PERTURB = 1,
185: SKEW = 2,
186: SKEW_PERTURB = 3
187: } Transform;
188: const char *const TransformTypes[] = {"none", "perturb", "skew", "skew_perturb", "Perturbation", "", NULL};
190: /* Enum and string array for selecting the form of the exact solution*/
191: typedef enum {
192: LINEAR = 0,
193: SINUSOIDAL = 1
194: } Solution;
195: const char *const SolutionTypes[] = {"linear", "sinusoidal", "Solution", "", NULL};
197: typedef struct {
198: Transform mesh_transform;
199: Solution sol_form;
200: } UserCtx;
202: /* Process command line options and initialize the UserCtx struct */
203: static PetscErrorCode ProcessOptions(MPI_Comm comm, UserCtx *user)
204: {
205: /* Default to 2D, unperturbed triangle mesh and Linear solution.*/
206: user->mesh_transform = NONE;
207: user->sol_form = LINEAR;
209: PetscOptionsBegin(comm, "", "H-div Test Options", "DMPLEX");
210: PetscOptionsEnum("-mesh_transform", "Method used to perturb the mesh vertices. Options are skew, perturb, skew_perturb,or none", "ex39.c", TransformTypes, (PetscEnum)user->mesh_transform, (PetscEnum *)&user->mesh_transform, NULL);
211: PetscOptionsEnum("-sol_form", "Form of the exact solution. Options are Linear or Sinusoidal", "ex39.c", SolutionTypes, (PetscEnum)user->sol_form, (PetscEnum *)&user->sol_form, NULL);
212: PetscOptionsEnd();
213: return 0;
214: }
216: /* Perturb the position of each mesh vertex by a small amount.*/
217: static PetscErrorCode PerturbMesh(DM *mesh, PetscScalar *coordVals, PetscInt npoints, PetscInt dim)
218: {
219: PetscInt i, j, k;
220: PetscReal minCoords[3], maxCoords[3], maxPert[3], randVal, amp;
221: PetscRandom ran;
223: DMGetCoordinateDim(*mesh, &dim);
224: DMGetLocalBoundingBox(*mesh, minCoords, maxCoords);
225: PetscRandomCreate(PETSC_COMM_WORLD, &ran);
227: /* Compute something approximately equal to half an edge length. This is the
228: * most we can perturb points and guarantee that there won't be any topology
229: * issues. */
230: for (k = 0; k < dim; ++k) maxPert[k] = 0.025 * (maxCoords[k] - minCoords[k]) / (PetscPowReal(npoints, 1. / dim) - 1);
231: /* For each mesh vertex */
232: for (i = 0; i < npoints; ++i) {
233: /* For each coordinate of the vertex */
234: for (j = 0; j < dim; ++j) {
235: /* Generate a random amplitude in [-0.5*maxPert, 0.5*maxPert] */
236: PetscRandomGetValueReal(ran, &randVal);
237: amp = maxPert[j] * (randVal - 0.5);
238: /* Add the perturbation to the vertex*/
239: coordVals[dim * i + j] += amp;
240: }
241: }
243: PetscRandomDestroy(&ran);
244: return 0;
245: }
247: /* Apply a global skew transformation to the mesh. */
248: static PetscErrorCode SkewMesh(DM *mesh, PetscScalar *coordVals, PetscInt npoints, PetscInt dim)
249: {
250: PetscInt i, j, k, l;
251: PetscScalar *transMat;
252: PetscScalar tmpcoord[3];
253: PetscRandom ran;
254: PetscReal randVal;
256: PetscCalloc1(dim * dim, &transMat);
257: PetscRandomCreate(PETSC_COMM_WORLD, &ran);
259: /* Make a matrix representing a skew transformation */
260: for (i = 0; i < dim; ++i) {
261: for (j = 0; j < dim; ++j) {
262: PetscRandomGetValueReal(ran, &randVal);
263: if (i == j) transMat[i * dim + j] = 1.;
264: else if (j < i) transMat[i * dim + j] = 2 * (j + i) * randVal;
265: else transMat[i * dim + j] = 0;
266: }
267: }
269: /* Multiply each coordinate vector by our transformation.*/
270: for (i = 0; i < npoints; ++i) {
271: for (j = 0; j < dim; ++j) {
272: tmpcoord[j] = 0;
273: for (k = 0; k < dim; ++k) tmpcoord[j] += coordVals[dim * i + k] * transMat[dim * k + j];
274: }
275: for (l = 0; l < dim; ++l) coordVals[dim * i + l] = tmpcoord[l];
276: }
277: PetscFree(transMat);
278: PetscRandomDestroy(&ran);
279: return 0;
280: }
282: /* Accesses the mesh coordinate array and performs the transformation operations
283: * specified by the user options */
284: static PetscErrorCode TransformMesh(UserCtx *user, DM *mesh)
285: {
286: PetscInt dim, npoints;
287: PetscScalar *coordVals;
288: Vec coords;
290: DMGetCoordinates(*mesh, &coords);
291: VecGetArray(coords, &coordVals);
292: VecGetLocalSize(coords, &npoints);
293: DMGetCoordinateDim(*mesh, &dim);
294: npoints = npoints / dim;
296: switch (user->mesh_transform) {
297: case PERTURB:
298: PerturbMesh(mesh, coordVals, npoints, dim);
299: break;
300: case SKEW:
301: SkewMesh(mesh, coordVals, npoints, dim);
302: break;
303: case SKEW_PERTURB:
304: SkewMesh(mesh, coordVals, npoints, dim);
305: PerturbMesh(mesh, coordVals, npoints, dim);
306: break;
307: default:
308: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "invalid mesh transformation");
309: }
310: VecRestoreArray(coords, &coordVals);
311: DMSetCoordinates(*mesh, coords);
312: return 0;
313: }
315: static PetscErrorCode CreateMesh(MPI_Comm comm, UserCtx *user, DM *mesh)
316: {
317: DMCreate(comm, mesh);
318: DMSetType(*mesh, DMPLEX);
319: DMSetFromOptions(*mesh);
321: /* Perform any mesh transformations if specified by user */
322: if (user->mesh_transform != NONE) TransformMesh(user, mesh);
324: /* Get any other mesh options from the command line */
325: DMSetApplicationContext(*mesh, user);
326: DMViewFromOptions(*mesh, NULL, "-dm_view");
327: return 0;
328: }
330: /* Setup the system of equations that we wish to solve */
331: static PetscErrorCode SetupProblem(DM dm, UserCtx *user)
332: {
333: PetscDS prob;
334: DMLabel label;
335: const PetscInt id = 1;
337: DMGetDS(dm, &prob);
338: /* All of these are independent of the user's choice of solution */
339: PetscDSSetResidual(prob, 1, f0_q, NULL);
340: PetscDSSetResidual(prob, 2, f0_w, NULL);
341: PetscDSSetJacobian(prob, 0, 0, g0_vu, NULL, NULL, NULL);
342: PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL);
343: PetscDSSetJacobian(prob, 1, 1, g0_qp, NULL, NULL, NULL);
344: PetscDSSetJacobian(prob, 2, 0, NULL, g1_wu, NULL, NULL);
345: PetscDSSetJacobian(prob, 2, 1, g0_wp, NULL, NULL, NULL);
346: PetscDSSetJacobian(prob, 2, 2, g0_wd, NULL, NULL, NULL);
348: /* Field 2 is the error between \div{u} and pressure in a higher dimensional
349: * space. If all is right this should be machine zero. */
350: PetscDSSetExactSolution(prob, 2, zero_func, NULL);
352: switch (user->sol_form) {
353: case LINEAR:
354: PetscDSSetResidual(prob, 0, f0_v_linear, NULL);
355: PetscDSSetBdResidual(prob, 0, f0_bd_u_linear, NULL);
356: PetscDSSetExactSolution(prob, 0, linear_u, NULL);
357: PetscDSSetExactSolution(prob, 1, linear_p, NULL);
358: break;
359: case SINUSOIDAL:
360: PetscDSSetResidual(prob, 0, f0_v_sinusoid, NULL);
361: PetscDSSetBdResidual(prob, 0, f0_bd_u_sinusoid, NULL);
362: PetscDSSetExactSolution(prob, 0, sinusoid_u, NULL);
363: PetscDSSetExactSolution(prob, 1, sinusoid_p, NULL);
364: break;
365: default:
366: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "invalid solution form");
367: }
369: DMGetLabel(dm, "marker", &label);
370: PetscDSAddBoundary(prob, DM_BC_NATURAL, "Boundary Integral", label, 1, &id, 0, 0, NULL, (void (*)(void))NULL, NULL, user, NULL);
371: return 0;
372: }
374: /* Create the finite element spaces we will use for this system */
375: static PetscErrorCode SetupDiscretization(DM mesh, PetscErrorCode (*setup)(DM, UserCtx *), UserCtx *user)
376: {
377: DM cdm = mesh;
378: PetscFE fevel, fepres, fedivErr;
379: PetscInt dim;
380: PetscBool simplex;
382: DMGetDimension(mesh, &dim);
383: DMPlexIsSimplex(mesh, &simplex);
384: /* Create FE objects and give them names so that options can be set from
385: * command line */
386: PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, dim, simplex, "velocity_", -1, &fevel);
387: PetscObjectSetName((PetscObject)fevel, "velocity");
389: PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, 1, simplex, "pressure_", -1, &fepres);
390: PetscObjectSetName((PetscObject)fepres, "pressure");
392: PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, 1, simplex, "divErr_", -1, &fedivErr);
393: PetscObjectSetName((PetscObject)fedivErr, "divErr");
395: PetscFECopyQuadrature(fevel, fepres);
396: PetscFECopyQuadrature(fevel, fedivErr);
398: /* Associate the FE objects with the mesh and setup the system */
399: DMSetField(mesh, 0, NULL, (PetscObject)fevel);
400: DMSetField(mesh, 1, NULL, (PetscObject)fepres);
401: DMSetField(mesh, 2, NULL, (PetscObject)fedivErr);
402: DMCreateDS(mesh);
403: (*setup)(mesh, user);
405: while (cdm) {
406: DMCopyDisc(mesh, cdm);
407: DMGetCoarseDM(cdm, &cdm);
408: }
410: /* The Mesh now owns the fields, so we can destroy the FEs created in this
411: * function */
412: PetscFEDestroy(&fevel);
413: PetscFEDestroy(&fepres);
414: PetscFEDestroy(&fedivErr);
415: DMDestroy(&cdm);
416: return 0;
417: }
419: int main(int argc, char **argv)
420: {
421: PetscInt i;
422: UserCtx user;
423: DM mesh;
424: SNES snes;
425: Vec computed, divErr;
426: PetscReal divErrNorm;
427: IS *fieldIS;
428: PetscBool exampleSuccess = PETSC_FALSE;
429: const PetscReal errTol = 10. * PETSC_SMALL;
431: char stdFormat[] = "L2 Norm of the Divergence Error is: %g\n H(div) elements working correctly: %s\n";
433: /* Initialize PETSc */
435: PetscInitialize(&argc, &argv, NULL, help);
436: ProcessOptions(PETSC_COMM_WORLD, &user);
438: /* Set up the system, we need to create a solver and a mesh and then assign
439: * the correct spaces into the mesh*/
440: SNESCreate(PETSC_COMM_WORLD, &snes);
441: CreateMesh(PETSC_COMM_WORLD, &user, &mesh);
442: SNESSetDM(snes, mesh);
443: SetupDiscretization(mesh, SetupProblem, &user);
444: DMPlexSetSNESLocalFEM(mesh, &user, &user, &user);
445: SNESSetFromOptions(snes);
447: /* Grab field IS so that we can view the solution by field */
448: DMCreateFieldIS(mesh, NULL, NULL, &fieldIS);
450: /* Create a vector to store the SNES solution, solve the system and grab the
451: * solution from SNES */
452: DMCreateGlobalVector(mesh, &computed);
453: PetscObjectSetName((PetscObject)computed, "computedSolution");
454: VecSet(computed, 0.0);
455: SNESSolve(snes, NULL, computed);
456: SNESGetSolution(snes, &computed);
457: VecViewFromOptions(computed, NULL, "-computedSolution_view");
459: /* Now we pull out the portion of the vector corresponding to the 3rd field
460: * which is the error between \div{u} computed in a higher dimensional space
461: * and p = \div{u} computed in a low dimension space. We report the L2 norm of
462: * this vector which should be zero if the H(div) spaces are implemented
463: * correctly. */
464: VecGetSubVector(computed, fieldIS[2], &divErr);
465: VecNorm(divErr, NORM_2, &divErrNorm);
466: VecRestoreSubVector(computed, fieldIS[2], &divErr);
467: exampleSuccess = (PetscBool)(divErrNorm <= errTol);
469: PetscPrintf(PETSC_COMM_WORLD, stdFormat, divErrNorm, exampleSuccess ? "true" : "false");
471: /* Tear down */
472: VecDestroy(&divErr);
473: VecDestroy(&computed);
474: for (i = 0; i < 3; ++i) ISDestroy(&fieldIS[i]);
475: PetscFree(fieldIS);
476: SNESDestroy(&snes);
477: DMDestroy(&mesh);
478: PetscFinalize();
479: return 0;
480: }
482: /*TEST
483: testset:
484: suffix: 2d_bdm
485: requires: triangle
486: args: -velocity_petscfe_default_quadrature_order 1 \
487: -velocity_petscspace_degree 1 \
488: -velocity_petscdualspace_type bdm \
489: -divErr_petscspace_degree 1 \
490: -divErr_petscdualspace_lagrange_continuity false \
491: -snes_error_if_not_converged \
492: -ksp_rtol 1e-10 \
493: -ksp_error_if_not_converged \
494: -pc_type fieldsplit\
495: -pc_fieldsplit_detect_saddle_point\
496: -pc_fieldsplit_type schur\
497: -pc_fieldsplit_schur_precondition full
498: test:
499: suffix: linear
500: args: -sol_form linear -mesh_transform none
501: test:
502: suffix: sinusoidal
503: args: -sol_form sinusoidal -mesh_transform none
504: test:
505: suffix: sinusoidal_skew
506: args: -sol_form sinusoidal -mesh_transform skew
507: test:
508: suffix: sinusoidal_perturb
509: args: -sol_form sinusoidal -mesh_transform perturb
510: test:
511: suffix: sinusoidal_skew_perturb
512: args: -sol_form sinusoidal -mesh_transform skew_perturb
514: testset:
515: TODO: broken
516: suffix: 2d_bdmq
517: args: -dm_plex_simplex false \
518: -velocity_petscspace_degree 1 \
519: -velocity_petscdualspace_type bdm \
520: -velocity_petscdualspace_lagrange_tensor 1 \
521: -divErr_petscspace_degree 1 \
522: -divErr_petscdualspace_lagrange_continuity false \
523: -snes_error_if_not_converged \
524: -ksp_rtol 1e-10 \
525: -ksp_error_if_not_converged \
526: -pc_type fieldsplit\
527: -pc_fieldsplit_detect_saddle_point\
528: -pc_fieldsplit_type schur\
529: -pc_fieldsplit_schur_precondition full
530: test:
531: suffix: linear
532: args: -sol_form linear -mesh_transform none
533: test:
534: suffix: sinusoidal
535: args: -sol_form sinusoidal -mesh_transform none
536: test:
537: suffix: sinusoidal_skew
538: args: -sol_form sinusoidal -mesh_transform skew
539: test:
540: suffix: sinusoidal_perturb
541: args: -sol_form sinusoidal -mesh_transform perturb
542: test:
543: suffix: sinusoidal_skew_perturb
544: args: -sol_form sinusoidal -mesh_transform skew_perturb
546: testset:
547: suffix: 3d_bdm
548: requires: ctetgen
549: args: -dm_plex_dim 3 \
550: -velocity_petscspace_degree 1 \
551: -velocity_petscdualspace_type bdm \
552: -divErr_petscspace_degree 1 \
553: -divErr_petscdualspace_lagrange_continuity false \
554: -snes_error_if_not_converged \
555: -ksp_rtol 1e-10 \
556: -ksp_error_if_not_converged \
557: -pc_type fieldsplit \
558: -pc_fieldsplit_detect_saddle_point \
559: -pc_fieldsplit_type schur \
560: -pc_fieldsplit_schur_precondition full
561: test:
562: suffix: linear
563: args: -sol_form linear -mesh_transform none
564: test:
565: suffix: sinusoidal
566: args: -sol_form sinusoidal -mesh_transform none
567: test:
568: suffix: sinusoidal_skew
569: args: -sol_form sinusoidal -mesh_transform skew
570: test:
571: suffix: sinusoidal_perturb
572: args: -sol_form sinusoidal -mesh_transform perturb
573: test:
574: suffix: sinusoidal_skew_perturb
575: args: -sol_form sinusoidal -mesh_transform skew_perturb
577: testset:
578: TODO: broken
579: suffix: 3d_bdmq
580: requires: ctetgen
581: args: -dm_plex_dim 3 \
582: -dm_plex_simplex false \
583: -velocity_petscspace_degree 1 \
584: -velocity_petscdualspace_type bdm \
585: -velocity_petscdualspace_lagrange_tensor 1 \
586: -divErr_petscspace_degree 1 \
587: -divErr_petscdualspace_lagrange_continuity false \
588: -snes_error_if_not_converged \
589: -ksp_rtol 1e-10 \
590: -ksp_error_if_not_converged \
591: -pc_type fieldsplit \
592: -pc_fieldsplit_detect_saddle_point \
593: -pc_fieldsplit_type schur \
594: -pc_fieldsplit_schur_precondition full
595: test:
596: suffix: linear
597: args: -sol_form linear -mesh_transform none
598: test:
599: suffix: sinusoidal
600: args: -sol_form sinusoidal -mesh_transform none
601: test:
602: suffix: sinusoidal_skew
603: args: -sol_form sinusoidal -mesh_transform skew
604: test:
605: suffix: sinusoidal_perturb
606: args: -sol_form sinusoidal -mesh_transform perturb
607: test:
608: suffix: sinusoidal_skew_perturb
609: args: -sol_form sinusoidal -mesh_transform skew_perturb
611: test:
612: suffix: quad_rt_0
613: args: -dm_plex_simplex false -mesh_transform skew \
614: -divErr_petscspace_degree 1 \
615: -divErr_petscdualspace_lagrange_continuity false \
616: -snes_error_if_not_converged \
617: -ksp_rtol 1e-10 \
618: -ksp_error_if_not_converged \
619: -pc_type fieldsplit\
620: -pc_fieldsplit_detect_saddle_point\
621: -pc_fieldsplit_type schur\
622: -pc_fieldsplit_schur_precondition full \
623: -velocity_petscfe_default_quadrature_order 1 \
624: -velocity_petscspace_type sum \
625: -velocity_petscspace_variables 2 \
626: -velocity_petscspace_components 2 \
627: -velocity_petscspace_sum_spaces 2 \
628: -velocity_petscspace_sum_concatenate true \
629: -velocity_sumcomp_0_petscspace_variables 2 \
630: -velocity_sumcomp_0_petscspace_type tensor \
631: -velocity_sumcomp_0_petscspace_tensor_spaces 2 \
632: -velocity_sumcomp_0_petscspace_tensor_uniform false \
633: -velocity_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
634: -velocity_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
635: -velocity_sumcomp_1_petscspace_variables 2 \
636: -velocity_sumcomp_1_petscspace_type tensor \
637: -velocity_sumcomp_1_petscspace_tensor_spaces 2 \
638: -velocity_sumcomp_1_petscspace_tensor_uniform false \
639: -velocity_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
640: -velocity_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
641: -velocity_petscdualspace_form_degree -1 \
642: -velocity_petscdualspace_order 1 \
643: -velocity_petscdualspace_lagrange_trimmed true
644: TEST*/