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, PetscCtx ctx)
25: {
26: PetscInt c;
27: for (c = 0; c < Nc; ++c) u[c] = 0;
28: return PETSC_SUCCESS;
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, PetscCtx ctx)
35: {
36: PetscInt c;
37: for (c = 0; c < Nc; ++c) u[c] = x[c];
38: return PETSC_SUCCESS;
39: }
40: static PetscErrorCode linear_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
41: {
42: u[0] = dim;
43: return PETSC_SUCCESS;
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, PetscCtx ctx)
52: {
53: PetscInt c;
54: for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(2 * PETSC_PI * x[c]);
55: return PETSC_SUCCESS;
56: }
57: static PetscErrorCode sinusoid_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx 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 PETSC_SUCCESS;
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: PetscCallAbort(PETSC_COMM_SELF, PetscCalloc1(dim, &u_rhs));
72: PetscCallAbort(PETSC_COMM_SELF, 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: PetscCallAbort(PETSC_COMM_SELF, 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: PetscCallAbort(PETSC_COMM_SELF, PetscCalloc1(dim, &u_rhs));
83: PetscCallAbort(PETSC_COMM_SELF, 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: PetscCallAbort(PETSC_COMM_SELF, 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;
117: PetscCallAbort(PETSC_COMM_SELF, PetscCalloc1(dim, &u_rhs));
118: PetscCallAbort(PETSC_COMM_SELF, sinusoid_u(dim, t, x, dim, u_rhs, NULL));
119: for (d = 0; d < dim; ++d) f0[d] = u_rhs[d];
120: PetscCallAbort(PETSC_COMM_SELF, 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;
128: PetscCallAbort(PETSC_COMM_SELF, PetscCalloc1(dim, &u_rhs));
129: PetscCallAbort(PETSC_COMM_SELF, linear_u(dim, t, x, dim, u_rhs, NULL));
130: for (d = 0; d < dim; ++d) f0[d] = u_rhs[d];
131: PetscCallAbort(PETSC_COMM_SELF, 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: for (PetscInt c = 0; c < dim; ++c) g0[c * dim + c] = 1.0;
140: }
142: /* <q, p> */
143: 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[])
144: {
145: PetscInt d;
146: for (d = 0; d < dim; ++d) g0[d * dim + d] = 1.0;
147: }
149: /* -<q, \div{u}> For the embedded system. This is different from the method of
150: * manufactured solution because instead of computing <q,\div{u}> - <q,f> we
151: * need <q,p> - <q,\div{u}.*/
152: 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[])
153: {
154: for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] = -1.0;
155: }
157: /* <w, p> This is only used by the embedded system. Where we need to compute
158: * <w,d> - <w,p> + <w, \div{u}>*/
159: 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[])
160: {
161: for (PetscInt d = 0; d < dim; ++d) g0[d * dim + d] = -1.0;
162: }
164: /* <w, d> */
165: 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[])
166: {
167: for (PetscInt c = 0; c < dim; ++c) g0[c * dim + c] = 1.0;
168: }
170: /* <w, \div{u}> for the embedded system. */
171: 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[])
172: {
173: for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
174: }
176: /* Enum and string array for selecting mesh perturbation options */
177: typedef enum {
178: NONE = 0,
179: PERTURB = 1,
180: SKEW = 2,
181: SKEW_PERTURB = 3
182: } Transform;
183: const char *const TransformTypes[] = {"none", "perturb", "skew", "skew_perturb", "Perturbation", "", NULL};
185: /* Enum and string array for selecting the form of the exact solution*/
186: typedef enum {
187: LINEAR = 0,
188: SINUSOIDAL = 1
189: } Solution;
190: const char *const SolutionTypes[] = {"linear", "sinusoidal", "Solution", "", NULL};
192: typedef struct {
193: Transform mesh_transform;
194: Solution sol_form;
195: } UserCtx;
197: /* Process command line options and initialize the UserCtx struct */
198: static PetscErrorCode ProcessOptions(MPI_Comm comm, UserCtx *user)
199: {
200: PetscFunctionBegin;
201: /* Default to 2D, unperturbed triangle mesh and Linear solution.*/
202: user->mesh_transform = NONE;
203: user->sol_form = LINEAR;
205: PetscOptionsBegin(comm, "", "H-div Test Options", "DMPLEX");
206: PetscCall(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));
207: PetscCall(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));
208: PetscOptionsEnd();
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /* Perturb the position of each mesh vertex by a small amount.*/
213: static PetscErrorCode PerturbMesh(DM *mesh, PetscScalar *coordVals, PetscInt npoints, PetscInt dim)
214: {
215: PetscInt i, j, k;
216: PetscReal minCoords[3], maxCoords[3], maxPert[3], randVal, amp;
217: PetscRandom ran;
219: PetscFunctionBegin;
220: PetscCall(DMGetCoordinateDim(*mesh, &dim));
221: PetscCall(DMGetLocalBoundingBox(*mesh, minCoords, maxCoords));
222: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &ran));
224: /* Compute something approximately equal to half an edge length. This is the
225: * most we can perturb points and guarantee that there won't be any topology
226: * issues. */
227: for (k = 0; k < dim; ++k) maxPert[k] = 0.025 * (maxCoords[k] - minCoords[k]) / (PetscPowReal(npoints, 1. / dim) - 1);
228: /* For each mesh vertex */
229: for (i = 0; i < npoints; ++i) {
230: /* For each coordinate of the vertex */
231: for (j = 0; j < dim; ++j) {
232: /* Generate a random amplitude in [-0.5*maxPert, 0.5*maxPert] */
233: PetscCall(PetscRandomGetValueReal(ran, &randVal));
234: amp = maxPert[j] * (randVal - 0.5);
235: /* Add the perturbation to the vertex*/
236: coordVals[dim * i + j] += amp;
237: }
238: }
240: PetscCall(PetscRandomDestroy(&ran));
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: /* Apply a global skew transformation to the mesh. */
245: static PetscErrorCode SkewMesh(DM *mesh, PetscScalar *coordVals, PetscInt npoints, PetscInt dim)
246: {
247: PetscInt i, j, k, l;
248: PetscScalar *transMat;
249: PetscScalar tmpcoord[3];
250: PetscRandom ran;
251: PetscReal randVal;
253: PetscFunctionBegin;
254: PetscCall(PetscCalloc1(dim * dim, &transMat));
255: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &ran));
257: /* Make a matrix representing a skew transformation */
258: for (i = 0; i < dim; ++i) {
259: for (j = 0; j < dim; ++j) {
260: PetscCall(PetscRandomGetValueReal(ran, &randVal));
261: if (i == j) transMat[i * dim + j] = 1.;
262: else if (j < i) transMat[i * dim + j] = 2 * (j + i) * randVal;
263: else transMat[i * dim + j] = 0;
264: }
265: }
267: /* Multiply each coordinate vector by our transformation.*/
268: for (i = 0; i < npoints; ++i) {
269: for (j = 0; j < dim; ++j) {
270: tmpcoord[j] = 0;
271: for (k = 0; k < dim; ++k) tmpcoord[j] += coordVals[dim * i + k] * transMat[dim * k + j];
272: }
273: for (l = 0; l < dim; ++l) coordVals[dim * i + l] = tmpcoord[l];
274: }
275: PetscCall(PetscFree(transMat));
276: PetscCall(PetscRandomDestroy(&ran));
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /* Accesses the mesh coordinate array and performs the transformation operations
281: * specified by the user options */
282: static PetscErrorCode TransformMesh(UserCtx *user, DM *mesh)
283: {
284: PetscInt dim, npoints;
285: PetscScalar *coordVals;
286: Vec coords;
288: PetscFunctionBegin;
289: PetscCall(DMGetCoordinates(*mesh, &coords));
290: PetscCall(VecGetArray(coords, &coordVals));
291: PetscCall(VecGetLocalSize(coords, &npoints));
292: PetscCall(DMGetCoordinateDim(*mesh, &dim));
293: npoints = npoints / dim;
295: switch (user->mesh_transform) {
296: case PERTURB:
297: PetscCall(PerturbMesh(mesh, coordVals, npoints, dim));
298: break;
299: case SKEW:
300: PetscCall(SkewMesh(mesh, coordVals, npoints, dim));
301: break;
302: case SKEW_PERTURB:
303: PetscCall(SkewMesh(mesh, coordVals, npoints, dim));
304: PetscCall(PerturbMesh(mesh, coordVals, npoints, dim));
305: break;
306: default:
307: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "invalid mesh transformation");
308: }
309: PetscCall(VecRestoreArray(coords, &coordVals));
310: PetscCall(DMSetCoordinates(*mesh, coords));
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: static PetscErrorCode CreateMesh(MPI_Comm comm, UserCtx *user, DM *mesh)
315: {
316: PetscFunctionBegin;
317: PetscCall(DMCreate(comm, mesh));
318: PetscCall(DMSetType(*mesh, DMPLEX));
319: PetscCall(DMSetFromOptions(*mesh));
321: /* Perform any mesh transformations if specified by user */
322: if (user->mesh_transform != NONE) PetscCall(TransformMesh(user, mesh));
324: /* Get any other mesh options from the command line */
325: PetscCall(DMSetApplicationContext(*mesh, user));
326: PetscCall(DMViewFromOptions(*mesh, NULL, "-dm_view"));
327: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
338: PetscCall(DMGetDS(dm, &prob));
339: /* All of these are independent of the user's choice of solution */
340: PetscCall(PetscDSSetResidual(prob, 1, f0_q, NULL));
341: PetscCall(PetscDSSetResidual(prob, 2, f0_w, NULL));
342: PetscCall(PetscDSSetJacobian(prob, 0, 0, g0_vu, NULL, NULL, NULL));
343: PetscCall(PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL));
344: PetscCall(PetscDSSetJacobian(prob, 1, 1, g0_qp, NULL, NULL, NULL));
345: PetscCall(PetscDSSetJacobian(prob, 2, 0, NULL, g1_wu, NULL, NULL));
346: PetscCall(PetscDSSetJacobian(prob, 2, 1, g0_wp, NULL, NULL, NULL));
347: PetscCall(PetscDSSetJacobian(prob, 2, 2, g0_wd, NULL, NULL, NULL));
349: /* Field 2 is the error between \div{u} and pressure in a higher dimensional
350: * space. If all is right this should be machine zero. */
351: PetscCall(PetscDSSetExactSolution(prob, 2, zero_func, NULL));
353: switch (user->sol_form) {
354: case LINEAR:
355: PetscCall(PetscDSSetResidual(prob, 0, f0_v_linear, NULL));
356: PetscCall(PetscDSSetBdResidual(prob, 0, f0_bd_u_linear, NULL));
357: PetscCall(PetscDSSetExactSolution(prob, 0, linear_u, NULL));
358: PetscCall(PetscDSSetExactSolution(prob, 1, linear_p, NULL));
359: break;
360: case SINUSOIDAL:
361: PetscCall(PetscDSSetResidual(prob, 0, f0_v_sinusoid, NULL));
362: PetscCall(PetscDSSetBdResidual(prob, 0, f0_bd_u_sinusoid, NULL));
363: PetscCall(PetscDSSetExactSolution(prob, 0, sinusoid_u, NULL));
364: PetscCall(PetscDSSetExactSolution(prob, 1, sinusoid_p, NULL));
365: break;
366: default:
367: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "invalid solution form");
368: }
370: PetscCall(DMGetLabel(dm, "marker", &label));
371: PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL, "Boundary Integral", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)NULL, NULL, user, NULL));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: /* Create the finite element spaces we will use for this system */
376: static PetscErrorCode SetupDiscretization(DM mesh, PetscErrorCode (*setup)(DM, UserCtx *), UserCtx *user)
377: {
378: DM cdm = mesh;
379: PetscFE fevel, fepres, fedivErr;
380: PetscInt dim;
381: PetscBool simplex;
383: PetscFunctionBegin;
384: PetscCall(DMGetDimension(mesh, &dim));
385: PetscCall(DMPlexIsSimplex(mesh, &simplex));
386: /* Create FE objects and give them names so that options can be set from
387: * command line */
388: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, dim, simplex, "velocity_", -1, &fevel));
389: PetscCall(PetscObjectSetName((PetscObject)fevel, "velocity"));
391: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, 1, simplex, "pressure_", -1, &fepres));
392: PetscCall(PetscObjectSetName((PetscObject)fepres, "pressure"));
394: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, 1, simplex, "divErr_", -1, &fedivErr));
395: PetscCall(PetscObjectSetName((PetscObject)fedivErr, "divErr"));
397: PetscCall(PetscFECopyQuadrature(fevel, fepres));
398: PetscCall(PetscFECopyQuadrature(fevel, fedivErr));
400: /* Associate the FE objects with the mesh and setup the system */
401: PetscCall(DMSetField(mesh, 0, NULL, (PetscObject)fevel));
402: PetscCall(DMSetField(mesh, 1, NULL, (PetscObject)fepres));
403: PetscCall(DMSetField(mesh, 2, NULL, (PetscObject)fedivErr));
404: PetscCall(DMCreateDS(mesh));
405: PetscCall((*setup)(mesh, user));
407: while (cdm) {
408: PetscCall(DMCopyDisc(mesh, cdm));
409: PetscCall(DMGetCoarseDM(cdm, &cdm));
410: }
412: /* The Mesh now owns the fields, so we can destroy the FEs created in this
413: * function */
414: PetscCall(PetscFEDestroy(&fevel));
415: PetscCall(PetscFEDestroy(&fepres));
416: PetscCall(PetscFEDestroy(&fedivErr));
417: PetscCall(DMDestroy(&cdm));
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: int main(int argc, char **argv)
422: {
423: UserCtx user;
424: DM mesh;
425: SNES snes;
426: Vec computed, divErr;
427: PetscReal divErrNorm;
428: IS *fieldIS;
429: PetscBool exampleSuccess = PETSC_FALSE;
430: const PetscReal errTol = 10. * PETSC_SMALL;
432: char stdFormat[] = "L2 Norm of the Divergence Error is: %g\n H(div) elements working correctly: %s\n";
434: /* Initialize PETSc */
435: PetscFunctionBeginUser;
436: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
437: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
439: /* Set up the system, we need to create a solver and a mesh and then assign
440: * the correct spaces into the mesh*/
441: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
442: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &mesh));
443: PetscCall(SNESSetDM(snes, mesh));
444: PetscCall(SetupDiscretization(mesh, SetupProblem, &user));
445: PetscCall(DMPlexSetSNESLocalFEM(mesh, PETSC_FALSE, &user));
446: PetscCall(SNESSetFromOptions(snes));
448: /* Grab field IS so that we can view the solution by field */
449: PetscCall(DMCreateFieldIS(mesh, NULL, NULL, &fieldIS));
451: /* Create a vector to store the SNES solution, solve the system and grab the
452: * solution from SNES */
453: PetscCall(DMCreateGlobalVector(mesh, &computed));
454: PetscCall(PetscObjectSetName((PetscObject)computed, "computedSolution"));
455: PetscCall(SNESSolve(snes, NULL, computed));
456: PetscCall(SNESGetSolution(snes, &computed));
457: PetscCall(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: PetscCall(VecGetSubVector(computed, fieldIS[2], &divErr));
465: PetscCall(VecNorm(divErr, NORM_2, &divErrNorm));
466: PetscCall(VecRestoreSubVector(computed, fieldIS[2], &divErr));
467: exampleSuccess = (PetscBool)(divErrNorm <= errTol);
469: PetscCall(PetscPrintf(PETSC_COMM_WORLD, stdFormat, divErrNorm, exampleSuccess ? "true" : "false"));
471: /* Tear down */
472: PetscCall(VecDestroy(&divErr));
473: PetscCall(VecDestroy(&computed));
474: for (PetscInt i = 0; i < 3; ++i) PetscCall(ISDestroy(&fieldIS[i]));
475: PetscCall(PetscFree(fieldIS));
476: PetscCall(SNESDestroy(&snes));
477: PetscCall(DMDestroy(&mesh));
478: PetscCall(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: output_file: output/empty.out
518: args: -dm_plex_simplex false \
519: -velocity_petscspace_degree 1 \
520: -velocity_petscdualspace_type bdm \
521: -velocity_petscdualspace_lagrange_tensor 1 \
522: -divErr_petscspace_degree 1 \
523: -divErr_petscdualspace_lagrange_continuity false \
524: -snes_error_if_not_converged \
525: -ksp_rtol 1e-10 \
526: -ksp_error_if_not_converged \
527: -pc_type fieldsplit\
528: -pc_fieldsplit_detect_saddle_point\
529: -pc_fieldsplit_type schur\
530: -pc_fieldsplit_schur_precondition full
531: test:
532: suffix: linear
533: args: -sol_form linear -mesh_transform none
534: test:
535: suffix: sinusoidal
536: args: -sol_form sinusoidal -mesh_transform none
537: test:
538: suffix: sinusoidal_skew
539: args: -sol_form sinusoidal -mesh_transform skew
540: test:
541: suffix: sinusoidal_perturb
542: args: -sol_form sinusoidal -mesh_transform perturb
543: test:
544: suffix: sinusoidal_skew_perturb
545: args: -sol_form sinusoidal -mesh_transform skew_perturb
547: testset:
548: suffix: 3d_bdm
549: requires: ctetgen
550: args: -dm_plex_dim 3 \
551: -velocity_petscspace_degree 1 \
552: -velocity_petscdualspace_type bdm \
553: -divErr_petscspace_degree 1 \
554: -divErr_petscdualspace_lagrange_continuity false \
555: -snes_error_if_not_converged \
556: -ksp_rtol 1e-10 \
557: -ksp_error_if_not_converged \
558: -pc_type fieldsplit \
559: -pc_fieldsplit_detect_saddle_point \
560: -pc_fieldsplit_type schur \
561: -pc_fieldsplit_schur_precondition full
562: test:
563: suffix: linear
564: args: -sol_form linear -mesh_transform none
565: test:
566: suffix: sinusoidal
567: args: -sol_form sinusoidal -mesh_transform none
568: test:
569: suffix: sinusoidal_skew
570: args: -sol_form sinusoidal -mesh_transform skew
571: test:
572: suffix: sinusoidal_perturb
573: args: -sol_form sinusoidal -mesh_transform perturb
574: test:
575: suffix: sinusoidal_skew_perturb
576: args: -sol_form sinusoidal -mesh_transform skew_perturb
578: testset:
579: TODO: broken
580: suffix: 3d_bdmq
581: output_file: output/empty.out
582: requires: ctetgen
583: args: -dm_plex_dim 3 \
584: -dm_plex_simplex false \
585: -velocity_petscspace_degree 1 \
586: -velocity_petscdualspace_type bdm \
587: -velocity_petscdualspace_lagrange_tensor 1 \
588: -divErr_petscspace_degree 1 \
589: -divErr_petscdualspace_lagrange_continuity false \
590: -snes_error_if_not_converged \
591: -ksp_rtol 1e-10 \
592: -ksp_error_if_not_converged \
593: -pc_type fieldsplit \
594: -pc_fieldsplit_detect_saddle_point \
595: -pc_fieldsplit_type schur \
596: -pc_fieldsplit_schur_precondition full
597: test:
598: suffix: linear
599: args: -sol_form linear -mesh_transform none
600: test:
601: suffix: sinusoidal
602: args: -sol_form sinusoidal -mesh_transform none
603: test:
604: suffix: sinusoidal_skew
605: args: -sol_form sinusoidal -mesh_transform skew
606: test:
607: suffix: sinusoidal_perturb
608: args: -sol_form sinusoidal -mesh_transform perturb
609: test:
610: suffix: sinusoidal_skew_perturb
611: args: -sol_form sinusoidal -mesh_transform skew_perturb
613: test:
614: suffix: quad_rt_0
615: args: -dm_plex_simplex false -mesh_transform skew \
616: -divErr_petscspace_degree 1 \
617: -divErr_petscdualspace_lagrange_continuity false \
618: -snes_error_if_not_converged \
619: -ksp_rtol 1e-10 \
620: -ksp_error_if_not_converged \
621: -pc_type fieldsplit\
622: -pc_fieldsplit_detect_saddle_point\
623: -pc_fieldsplit_type schur\
624: -pc_fieldsplit_schur_precondition full \
625: -velocity_petscfe_default_quadrature_order 1 \
626: -velocity_petscspace_type sum \
627: -velocity_petscspace_variables 2 \
628: -velocity_petscspace_components 2 \
629: -velocity_petscspace_sum_spaces 2 \
630: -velocity_petscspace_sum_concatenate true \
631: -velocity_sumcomp_0_petscspace_variables 2 \
632: -velocity_sumcomp_0_petscspace_type tensor \
633: -velocity_sumcomp_0_petscspace_tensor_spaces 2 \
634: -velocity_sumcomp_0_petscspace_tensor_uniform false \
635: -velocity_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
636: -velocity_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
637: -velocity_sumcomp_1_petscspace_variables 2 \
638: -velocity_sumcomp_1_petscspace_type tensor \
639: -velocity_sumcomp_1_petscspace_tensor_spaces 2 \
640: -velocity_sumcomp_1_petscspace_tensor_uniform false \
641: -velocity_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
642: -velocity_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
643: -velocity_petscdualspace_form_degree -1 \
644: -velocity_petscdualspace_order 1 \
645: -velocity_petscdualspace_lagrange_trimmed true
646: TEST*/