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 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, void *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, void *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, void *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, 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 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: 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: PetscFunctionBegin;
206: /* Default to 2D, unperturbed triangle mesh and Linear solution.*/
207: user->mesh_transform = NONE;
208: user->sol_form = LINEAR;
210: PetscOptionsBegin(comm, "", "H-div Test Options", "DMPLEX");
211: 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));
212: 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));
213: PetscOptionsEnd();
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: /* Perturb the position of each mesh vertex by a small amount.*/
218: static PetscErrorCode PerturbMesh(DM *mesh, PetscScalar *coordVals, PetscInt npoints, PetscInt dim)
219: {
220: PetscInt i, j, k;
221: PetscReal minCoords[3], maxCoords[3], maxPert[3], randVal, amp;
222: PetscRandom ran;
224: PetscFunctionBegin;
225: PetscCall(DMGetCoordinateDim(*mesh, &dim));
226: PetscCall(DMGetLocalBoundingBox(*mesh, minCoords, maxCoords));
227: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &ran));
229: /* Compute something approximately equal to half an edge length. This is the
230: * most we can perturb points and guarantee that there won't be any topology
231: * issues. */
232: for (k = 0; k < dim; ++k) maxPert[k] = 0.025 * (maxCoords[k] - minCoords[k]) / (PetscPowReal(npoints, 1. / dim) - 1);
233: /* For each mesh vertex */
234: for (i = 0; i < npoints; ++i) {
235: /* For each coordinate of the vertex */
236: for (j = 0; j < dim; ++j) {
237: /* Generate a random amplitude in [-0.5*maxPert, 0.5*maxPert] */
238: PetscCall(PetscRandomGetValueReal(ran, &randVal));
239: amp = maxPert[j] * (randVal - 0.5);
240: /* Add the perturbation to the vertex*/
241: coordVals[dim * i + j] += amp;
242: }
243: }
245: PetscCall(PetscRandomDestroy(&ran));
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: /* Apply a global skew transformation to the mesh. */
250: static PetscErrorCode SkewMesh(DM *mesh, PetscScalar *coordVals, PetscInt npoints, PetscInt dim)
251: {
252: PetscInt i, j, k, l;
253: PetscScalar *transMat;
254: PetscScalar tmpcoord[3];
255: PetscRandom ran;
256: PetscReal randVal;
258: PetscFunctionBegin;
259: PetscCall(PetscCalloc1(dim * dim, &transMat));
260: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &ran));
262: /* Make a matrix representing a skew transformation */
263: for (i = 0; i < dim; ++i) {
264: for (j = 0; j < dim; ++j) {
265: PetscCall(PetscRandomGetValueReal(ran, &randVal));
266: if (i == j) transMat[i * dim + j] = 1.;
267: else if (j < i) transMat[i * dim + j] = 2 * (j + i) * randVal;
268: else transMat[i * dim + j] = 0;
269: }
270: }
272: /* Multiply each coordinate vector by our transformation.*/
273: for (i = 0; i < npoints; ++i) {
274: for (j = 0; j < dim; ++j) {
275: tmpcoord[j] = 0;
276: for (k = 0; k < dim; ++k) tmpcoord[j] += coordVals[dim * i + k] * transMat[dim * k + j];
277: }
278: for (l = 0; l < dim; ++l) coordVals[dim * i + l] = tmpcoord[l];
279: }
280: PetscCall(PetscFree(transMat));
281: PetscCall(PetscRandomDestroy(&ran));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /* Accesses the mesh coordinate array and performs the transformation operations
286: * specified by the user options */
287: static PetscErrorCode TransformMesh(UserCtx *user, DM *mesh)
288: {
289: PetscInt dim, npoints;
290: PetscScalar *coordVals;
291: Vec coords;
293: PetscFunctionBegin;
294: PetscCall(DMGetCoordinates(*mesh, &coords));
295: PetscCall(VecGetArray(coords, &coordVals));
296: PetscCall(VecGetLocalSize(coords, &npoints));
297: PetscCall(DMGetCoordinateDim(*mesh, &dim));
298: npoints = npoints / dim;
300: switch (user->mesh_transform) {
301: case PERTURB:
302: PetscCall(PerturbMesh(mesh, coordVals, npoints, dim));
303: break;
304: case SKEW:
305: PetscCall(SkewMesh(mesh, coordVals, npoints, dim));
306: break;
307: case SKEW_PERTURB:
308: PetscCall(SkewMesh(mesh, coordVals, npoints, dim));
309: PetscCall(PerturbMesh(mesh, coordVals, npoints, dim));
310: break;
311: default:
312: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "invalid mesh transformation");
313: }
314: PetscCall(VecRestoreArray(coords, &coordVals));
315: PetscCall(DMSetCoordinates(*mesh, coords));
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: static PetscErrorCode CreateMesh(MPI_Comm comm, UserCtx *user, DM *mesh)
320: {
321: PetscFunctionBegin;
322: PetscCall(DMCreate(comm, mesh));
323: PetscCall(DMSetType(*mesh, DMPLEX));
324: PetscCall(DMSetFromOptions(*mesh));
326: /* Perform any mesh transformations if specified by user */
327: if (user->mesh_transform != NONE) PetscCall(TransformMesh(user, mesh));
329: /* Get any other mesh options from the command line */
330: PetscCall(DMSetApplicationContext(*mesh, user));
331: PetscCall(DMViewFromOptions(*mesh, NULL, "-dm_view"));
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: /* Setup the system of equations that we wish to solve */
336: static PetscErrorCode SetupProblem(DM dm, UserCtx *user)
337: {
338: PetscDS prob;
339: DMLabel label;
340: const PetscInt id = 1;
342: PetscFunctionBegin;
343: PetscCall(DMGetDS(dm, &prob));
344: /* All of these are independent of the user's choice of solution */
345: PetscCall(PetscDSSetResidual(prob, 1, f0_q, NULL));
346: PetscCall(PetscDSSetResidual(prob, 2, f0_w, NULL));
347: PetscCall(PetscDSSetJacobian(prob, 0, 0, g0_vu, NULL, NULL, NULL));
348: PetscCall(PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL));
349: PetscCall(PetscDSSetJacobian(prob, 1, 1, g0_qp, NULL, NULL, NULL));
350: PetscCall(PetscDSSetJacobian(prob, 2, 0, NULL, g1_wu, NULL, NULL));
351: PetscCall(PetscDSSetJacobian(prob, 2, 1, g0_wp, NULL, NULL, NULL));
352: PetscCall(PetscDSSetJacobian(prob, 2, 2, g0_wd, NULL, NULL, NULL));
354: /* Field 2 is the error between \div{u} and pressure in a higher dimensional
355: * space. If all is right this should be machine zero. */
356: PetscCall(PetscDSSetExactSolution(prob, 2, zero_func, NULL));
358: switch (user->sol_form) {
359: case LINEAR:
360: PetscCall(PetscDSSetResidual(prob, 0, f0_v_linear, NULL));
361: PetscCall(PetscDSSetBdResidual(prob, 0, f0_bd_u_linear, NULL));
362: PetscCall(PetscDSSetExactSolution(prob, 0, linear_u, NULL));
363: PetscCall(PetscDSSetExactSolution(prob, 1, linear_p, NULL));
364: break;
365: case SINUSOIDAL:
366: PetscCall(PetscDSSetResidual(prob, 0, f0_v_sinusoid, NULL));
367: PetscCall(PetscDSSetBdResidual(prob, 0, f0_bd_u_sinusoid, NULL));
368: PetscCall(PetscDSSetExactSolution(prob, 0, sinusoid_u, NULL));
369: PetscCall(PetscDSSetExactSolution(prob, 1, sinusoid_p, NULL));
370: break;
371: default:
372: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "invalid solution form");
373: }
375: PetscCall(DMGetLabel(dm, "marker", &label));
376: PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL, "Boundary Integral", label, 1, &id, 0, 0, NULL, (void (*)(void))NULL, NULL, user, NULL));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /* Create the finite element spaces we will use for this system */
381: static PetscErrorCode SetupDiscretization(DM mesh, PetscErrorCode (*setup)(DM, UserCtx *), UserCtx *user)
382: {
383: DM cdm = mesh;
384: PetscFE fevel, fepres, fedivErr;
385: PetscInt dim;
386: PetscBool simplex;
388: PetscFunctionBegin;
389: PetscCall(DMGetDimension(mesh, &dim));
390: PetscCall(DMPlexIsSimplex(mesh, &simplex));
391: /* Create FE objects and give them names so that options can be set from
392: * command line */
393: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, dim, simplex, "velocity_", -1, &fevel));
394: PetscCall(PetscObjectSetName((PetscObject)fevel, "velocity"));
396: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, 1, simplex, "pressure_", -1, &fepres));
397: PetscCall(PetscObjectSetName((PetscObject)fepres, "pressure"));
399: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh), dim, 1, simplex, "divErr_", -1, &fedivErr));
400: PetscCall(PetscObjectSetName((PetscObject)fedivErr, "divErr"));
402: PetscCall(PetscFECopyQuadrature(fevel, fepres));
403: PetscCall(PetscFECopyQuadrature(fevel, fedivErr));
405: /* Associate the FE objects with the mesh and setup the system */
406: PetscCall(DMSetField(mesh, 0, NULL, (PetscObject)fevel));
407: PetscCall(DMSetField(mesh, 1, NULL, (PetscObject)fepres));
408: PetscCall(DMSetField(mesh, 2, NULL, (PetscObject)fedivErr));
409: PetscCall(DMCreateDS(mesh));
410: PetscCall((*setup)(mesh, user));
412: while (cdm) {
413: PetscCall(DMCopyDisc(mesh, cdm));
414: PetscCall(DMGetCoarseDM(cdm, &cdm));
415: }
417: /* The Mesh now owns the fields, so we can destroy the FEs created in this
418: * function */
419: PetscCall(PetscFEDestroy(&fevel));
420: PetscCall(PetscFEDestroy(&fepres));
421: PetscCall(PetscFEDestroy(&fedivErr));
422: PetscCall(DMDestroy(&cdm));
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: int main(int argc, char **argv)
427: {
428: PetscInt i;
429: UserCtx user;
430: DM mesh;
431: SNES snes;
432: Vec computed, divErr;
433: PetscReal divErrNorm;
434: IS *fieldIS;
435: PetscBool exampleSuccess = PETSC_FALSE;
436: const PetscReal errTol = 10. * PETSC_SMALL;
438: char stdFormat[] = "L2 Norm of the Divergence Error is: %g\n H(div) elements working correctly: %s\n";
440: /* Initialize PETSc */
441: PetscFunctionBeginUser;
442: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
443: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
445: /* Set up the system, we need to create a solver and a mesh and then assign
446: * the correct spaces into the mesh*/
447: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
448: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &mesh));
449: PetscCall(SNESSetDM(snes, mesh));
450: PetscCall(SetupDiscretization(mesh, SetupProblem, &user));
451: PetscCall(DMPlexSetSNESLocalFEM(mesh, &user, &user, &user));
452: PetscCall(SNESSetFromOptions(snes));
454: /* Grab field IS so that we can view the solution by field */
455: PetscCall(DMCreateFieldIS(mesh, NULL, NULL, &fieldIS));
457: /* Create a vector to store the SNES solution, solve the system and grab the
458: * solution from SNES */
459: PetscCall(DMCreateGlobalVector(mesh, &computed));
460: PetscCall(PetscObjectSetName((PetscObject)computed, "computedSolution"));
461: PetscCall(VecSet(computed, 0.0));
462: PetscCall(SNESSolve(snes, NULL, computed));
463: PetscCall(SNESGetSolution(snes, &computed));
464: PetscCall(VecViewFromOptions(computed, NULL, "-computedSolution_view"));
466: /* Now we pull out the portion of the vector corresponding to the 3rd field
467: * which is the error between \div{u} computed in a higher dimensional space
468: * and p = \div{u} computed in a low dimension space. We report the L2 norm of
469: * this vector which should be zero if the H(div) spaces are implemented
470: * correctly. */
471: PetscCall(VecGetSubVector(computed, fieldIS[2], &divErr));
472: PetscCall(VecNorm(divErr, NORM_2, &divErrNorm));
473: PetscCall(VecRestoreSubVector(computed, fieldIS[2], &divErr));
474: exampleSuccess = (PetscBool)(divErrNorm <= errTol);
476: PetscCall(PetscPrintf(PETSC_COMM_WORLD, stdFormat, divErrNorm, exampleSuccess ? "true" : "false"));
478: /* Tear down */
479: PetscCall(VecDestroy(&divErr));
480: PetscCall(VecDestroy(&computed));
481: for (i = 0; i < 3; ++i) PetscCall(ISDestroy(&fieldIS[i]));
482: PetscCall(PetscFree(fieldIS));
483: PetscCall(SNESDestroy(&snes));
484: PetscCall(DMDestroy(&mesh));
485: PetscCall(PetscFinalize());
486: return 0;
487: }
489: /*TEST
490: testset:
491: suffix: 2d_bdm
492: requires: triangle
493: args: -velocity_petscfe_default_quadrature_order 1 \
494: -velocity_petscspace_degree 1 \
495: -velocity_petscdualspace_type bdm \
496: -divErr_petscspace_degree 1 \
497: -divErr_petscdualspace_lagrange_continuity false \
498: -snes_error_if_not_converged \
499: -ksp_rtol 1e-10 \
500: -ksp_error_if_not_converged \
501: -pc_type fieldsplit\
502: -pc_fieldsplit_detect_saddle_point\
503: -pc_fieldsplit_type schur\
504: -pc_fieldsplit_schur_precondition full
505: test:
506: suffix: linear
507: args: -sol_form linear -mesh_transform none
508: test:
509: suffix: sinusoidal
510: args: -sol_form sinusoidal -mesh_transform none
511: test:
512: suffix: sinusoidal_skew
513: args: -sol_form sinusoidal -mesh_transform skew
514: test:
515: suffix: sinusoidal_perturb
516: args: -sol_form sinusoidal -mesh_transform perturb
517: test:
518: suffix: sinusoidal_skew_perturb
519: args: -sol_form sinusoidal -mesh_transform skew_perturb
521: testset:
522: TODO: broken
523: suffix: 2d_bdmq
524: args: -dm_plex_simplex false \
525: -velocity_petscspace_degree 1 \
526: -velocity_petscdualspace_type bdm \
527: -velocity_petscdualspace_lagrange_tensor 1 \
528: -divErr_petscspace_degree 1 \
529: -divErr_petscdualspace_lagrange_continuity false \
530: -snes_error_if_not_converged \
531: -ksp_rtol 1e-10 \
532: -ksp_error_if_not_converged \
533: -pc_type fieldsplit\
534: -pc_fieldsplit_detect_saddle_point\
535: -pc_fieldsplit_type schur\
536: -pc_fieldsplit_schur_precondition full
537: test:
538: suffix: linear
539: args: -sol_form linear -mesh_transform none
540: test:
541: suffix: sinusoidal
542: args: -sol_form sinusoidal -mesh_transform none
543: test:
544: suffix: sinusoidal_skew
545: args: -sol_form sinusoidal -mesh_transform skew
546: test:
547: suffix: sinusoidal_perturb
548: args: -sol_form sinusoidal -mesh_transform perturb
549: test:
550: suffix: sinusoidal_skew_perturb
551: args: -sol_form sinusoidal -mesh_transform skew_perturb
553: testset:
554: suffix: 3d_bdm
555: requires: ctetgen
556: args: -dm_plex_dim 3 \
557: -velocity_petscspace_degree 1 \
558: -velocity_petscdualspace_type bdm \
559: -divErr_petscspace_degree 1 \
560: -divErr_petscdualspace_lagrange_continuity false \
561: -snes_error_if_not_converged \
562: -ksp_rtol 1e-10 \
563: -ksp_error_if_not_converged \
564: -pc_type fieldsplit \
565: -pc_fieldsplit_detect_saddle_point \
566: -pc_fieldsplit_type schur \
567: -pc_fieldsplit_schur_precondition full
568: test:
569: suffix: linear
570: args: -sol_form linear -mesh_transform none
571: test:
572: suffix: sinusoidal
573: args: -sol_form sinusoidal -mesh_transform none
574: test:
575: suffix: sinusoidal_skew
576: args: -sol_form sinusoidal -mesh_transform skew
577: test:
578: suffix: sinusoidal_perturb
579: args: -sol_form sinusoidal -mesh_transform perturb
580: test:
581: suffix: sinusoidal_skew_perturb
582: args: -sol_form sinusoidal -mesh_transform skew_perturb
584: testset:
585: TODO: broken
586: suffix: 3d_bdmq
587: requires: ctetgen
588: args: -dm_plex_dim 3 \
589: -dm_plex_simplex false \
590: -velocity_petscspace_degree 1 \
591: -velocity_petscdualspace_type bdm \
592: -velocity_petscdualspace_lagrange_tensor 1 \
593: -divErr_petscspace_degree 1 \
594: -divErr_petscdualspace_lagrange_continuity false \
595: -snes_error_if_not_converged \
596: -ksp_rtol 1e-10 \
597: -ksp_error_if_not_converged \
598: -pc_type fieldsplit \
599: -pc_fieldsplit_detect_saddle_point \
600: -pc_fieldsplit_type schur \
601: -pc_fieldsplit_schur_precondition full
602: test:
603: suffix: linear
604: args: -sol_form linear -mesh_transform none
605: test:
606: suffix: sinusoidal
607: args: -sol_form sinusoidal -mesh_transform none
608: test:
609: suffix: sinusoidal_skew
610: args: -sol_form sinusoidal -mesh_transform skew
611: test:
612: suffix: sinusoidal_perturb
613: args: -sol_form sinusoidal -mesh_transform perturb
614: test:
615: suffix: sinusoidal_skew_perturb
616: args: -sol_form sinusoidal -mesh_transform skew_perturb
618: test:
619: suffix: quad_rt_0
620: args: -dm_plex_simplex false -mesh_transform skew \
621: -divErr_petscspace_degree 1 \
622: -divErr_petscdualspace_lagrange_continuity false \
623: -snes_error_if_not_converged \
624: -ksp_rtol 1e-10 \
625: -ksp_error_if_not_converged \
626: -pc_type fieldsplit\
627: -pc_fieldsplit_detect_saddle_point\
628: -pc_fieldsplit_type schur\
629: -pc_fieldsplit_schur_precondition full \
630: -velocity_petscfe_default_quadrature_order 1 \
631: -velocity_petscspace_type sum \
632: -velocity_petscspace_variables 2 \
633: -velocity_petscspace_components 2 \
634: -velocity_petscspace_sum_spaces 2 \
635: -velocity_petscspace_sum_concatenate true \
636: -velocity_sumcomp_0_petscspace_variables 2 \
637: -velocity_sumcomp_0_petscspace_type tensor \
638: -velocity_sumcomp_0_petscspace_tensor_spaces 2 \
639: -velocity_sumcomp_0_petscspace_tensor_uniform false \
640: -velocity_sumcomp_0_tensorcomp_0_petscspace_degree 1 \
641: -velocity_sumcomp_0_tensorcomp_1_petscspace_degree 0 \
642: -velocity_sumcomp_1_petscspace_variables 2 \
643: -velocity_sumcomp_1_petscspace_type tensor \
644: -velocity_sumcomp_1_petscspace_tensor_spaces 2 \
645: -velocity_sumcomp_1_petscspace_tensor_uniform false \
646: -velocity_sumcomp_1_tensorcomp_0_petscspace_degree 0 \
647: -velocity_sumcomp_1_tensorcomp_1_petscspace_degree 1 \
648: -velocity_petscdualspace_form_degree -1 \
649: -velocity_petscdualspace_order 1 \
650: -velocity_petscdualspace_lagrange_trimmed true
651: TEST*/