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