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*/