Actual source code: ex77.c

  1: static char help[] = "Nonlinear elasticity problem in 3d with simplicial finite elements.\n\
  2: We solve a nonlinear elasticity problem, modelled as an incompressible Neo-Hookean solid, \n\
  3:  with pressure loading in a rectangular domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";

  5: /*
  6: Nonlinear elasticity problem, which we discretize using the finite
  7: element method on an unstructured mesh. This uses both Dirichlet boundary conditions (fixed faces)
  8: and nonlinear Neumann boundary conditions (pressure loading).
  9: The Lagrangian density (modulo boundary conditions) for this problem is given by
 10: \begin{equation}
 11:   \frac{\mu}{2} (\mathrm{Tr}{C}-3) + J p + \frac{\kappa}{2} (J-1).
 12: \end{equation}

 14: Discretization:

 16: We use PetscFE to generate a tabulation of the finite element basis functions
 17: at quadrature points. We can currently generate an arbitrary order Lagrange
 18: element.

 20: Field Data:

 22:   DMPLEX data is organized by point, and the closure operation just stacks up the
 23: data from each sieve point in the closure. Thus, for a P_2-P_1 Stokes element, we
 24: have

 26:   cl{e} = {f e_0 e_1 e_2 v_0 v_1 v_2}
 27:   x     = [u_{e_0} v_{e_0} u_{e_1} v_{e_1} u_{e_2} v_{e_2} u_{v_0} v_{v_0} p_{v_0} u_{v_1} v_{v_1} p_{v_1} u_{v_2} v_{v_2} p_{v_2}]

 29: The problem here is that we would like to loop over each field separately for
 30: integration. Therefore, the closure visitor in DMPlexVecGetClosure() reorders
 31: the data so that each field is contiguous

 33:   x'    = [u_{e_0} v_{e_0} u_{e_1} v_{e_1} u_{e_2} v_{e_2} u_{v_0} v_{v_0} u_{v_1} v_{v_1} u_{v_2} v_{v_2} p_{v_0} p_{v_1} p_{v_2}]

 35: Likewise, DMPlexVecSetClosure() takes data partitioned by field, and correctly
 36: puts it into the Sieve ordering.

 38: */

 40: #include <petscdmplex.h>
 41: #include <petscsnes.h>
 42: #include <petscds.h>

 44: typedef enum {RUN_FULL, RUN_TEST} RunType;

 46: typedef struct {
 47:   RunType   runType; /* Whether to run tests, or solve the full problem */
 48:   PetscReal mu;      /* The shear modulus */
 49:   PetscReal p_wall;  /* The wall pressure */
 50: } AppCtx;

 52: #if 0
 53: static inline void Det2D(PetscReal *detJ, const PetscReal J[])
 54: {
 55:   *detJ = J[0]*J[3] - J[1]*J[2];
 56: }
 57: #endif

 59: static inline void Det3D(PetscReal *detJ, const PetscScalar J[])
 60: {
 61:   *detJ = PetscRealPart(J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
 62:                         J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
 63:                         J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
 64: }

 66: #if 0
 67: static inline void Cof2D(PetscReal C[], const PetscReal A[])
 68: {
 69:   C[0] =  A[3];
 70:   C[1] = -A[2];
 71:   C[2] = -A[1];
 72:   C[3] =  A[0];
 73: }
 74: #endif

 76: static inline void Cof3D(PetscReal C[], const PetscScalar A[])
 77: {
 78:   C[0*3+0] = PetscRealPart(A[1*3+1]*A[2*3+2] - A[1*3+2]*A[2*3+1]);
 79:   C[0*3+1] = PetscRealPart(A[1*3+2]*A[2*3+0] - A[1*3+0]*A[2*3+2]);
 80:   C[0*3+2] = PetscRealPart(A[1*3+0]*A[2*3+1] - A[1*3+1]*A[2*3+0]);
 81:   C[1*3+0] = PetscRealPart(A[0*3+2]*A[2*3+1] - A[0*3+1]*A[2*3+2]);
 82:   C[1*3+1] = PetscRealPart(A[0*3+0]*A[2*3+2] - A[0*3+2]*A[2*3+0]);
 83:   C[1*3+2] = PetscRealPart(A[0*3+1]*A[2*3+0] - A[0*3+0]*A[2*3+1]);
 84:   C[2*3+0] = PetscRealPart(A[0*3+1]*A[1*3+2] - A[0*3+2]*A[1*3+1]);
 85:   C[2*3+1] = PetscRealPart(A[0*3+2]*A[1*3+0] - A[0*3+0]*A[1*3+2]);
 86:   C[2*3+2] = PetscRealPart(A[0*3+0]*A[1*3+1] - A[0*3+1]*A[1*3+0]);
 87: }

 89: PetscErrorCode zero_scalar(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 90: {
 91:   u[0] = 0.0;
 92:   return 0;
 93: }

 95: PetscErrorCode zero_vector(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 96: {
 97:   const PetscInt Ncomp = dim;

 99:   PetscInt       comp;
100:   for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0.0;
101:   return 0;
102: }

104: PetscErrorCode coordinates(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
105: {
106:   const PetscInt Ncomp = dim;

108:   PetscInt       comp;
109:   for (comp = 0; comp < Ncomp; ++comp) u[comp] = x[comp];
110:   return 0;
111: }

113: PetscErrorCode elasticityMaterial(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
114: {
115:   AppCtx *user = (AppCtx *) ctx;
116:   u[0] = user->mu;
117:   return 0;
118: }

120: PetscErrorCode wallPressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
121: {
122:   AppCtx *user = (AppCtx *) ctx;
123:   u[0] = user->p_wall;
124:   return 0;
125: }

127: void f1_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
128:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
129:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
130:           PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
131: {
132:   const PetscInt  Ncomp = dim;
133:   const PetscReal mu = PetscRealPart(a[0]), kappa = 3.0;
134:   PetscReal       cofu_x[9/*Ncomp*dim*/], detu_x, p = PetscRealPart(u[Ncomp]);
135:   PetscInt        comp, d;

137:   Cof3D(cofu_x, u_x);
138:   Det3D(&detu_x, u_x);
139:   p += kappa * (detu_x - 1.0);

141:   /* f1 is the first Piola-Kirchhoff tensor */
142:   for (comp = 0; comp < Ncomp; ++comp) {
143:     for (d = 0; d < dim; ++d) {
144:       f1[comp*dim+d] = mu * u_x[comp*dim+d] + p * cofu_x[comp*dim+d];
145:     }
146:   }
147: }

149: void g3_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
150:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
151:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
152:           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
153: {
154:   const PetscInt  Ncomp = dim;
155:   const PetscReal mu = PetscRealPart(a[0]), kappa = 3.0;
156:   PetscReal       cofu_x[9/*Ncomp*dim*/], detu_x, pp, pm, p = PetscRealPart(u[Ncomp]);
157:   PetscInt        compI, compJ, d1, d2;

159:   Cof3D(cofu_x, u_x);
160:   Det3D(&detu_x, u_x);
161:   p += kappa * (detu_x - 1.0);
162:   pp = p/detu_x + kappa;
163:   pm = p/detu_x;

165:   /* g3 is the first elasticity tensor, i.e. A_i^I_j^J = S^{IJ} g_{ij} + C^{KILJ} F^k_K F^l_L g_{kj} g_{li} */
166:   for (compI = 0; compI < Ncomp; ++compI) {
167:     for (compJ = 0; compJ < Ncomp; ++compJ) {
168:       const PetscReal G = (compI == compJ) ? 1.0 : 0.0;

170:       for (d1 = 0; d1 < dim; ++d1) {
171:         for (d2 = 0; d2 < dim; ++d2) {
172:           const PetscReal g = (d1 == d2) ? 1.0 : 0.0;

174:           g3[((compI*Ncomp+compJ)*dim+d1)*dim+d2] = g * G * mu + pp * cofu_x[compI*dim+d1] * cofu_x[compJ*dim+d2] - pm * cofu_x[compI*dim+d2] * cofu_x[compJ*dim+d1];
175:         }
176:       }
177:     }
178:   }
179: }

181: void f0_bd_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
182:     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
183:     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
184:     PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
185: {
186:   const PetscInt    Ncomp = dim;
187:   const PetscScalar p = a[aOff[1]];
188:   PetscReal         cofu_x[9/*Ncomp*dim*/];
189:   PetscInt          comp, d;

191:   Cof3D(cofu_x, u_x);
192:   for (comp = 0; comp < Ncomp; ++comp) {
193:     for (d = 0, f0[comp] = 0.0; d < dim; ++d) f0[comp] += cofu_x[comp*dim+d] * n[d];
194:     f0[comp] *= p;
195:   }
196: }

198: void g1_bd_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
199:     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
200:     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
201:     PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
202: {
203:   const PetscInt Ncomp = dim;
204:   PetscScalar    p = a[aOff[1]];
205:   PetscReal      cofu_x[9/*Ncomp*dim*/], m[3/*Ncomp*/], detu_x;
206:   PetscInt       comp, compI, compJ, d;

208:   Cof3D(cofu_x, u_x);
209:   Det3D(&detu_x, u_x);
210:   p /= detu_x;

212:   for (comp = 0; comp < Ncomp; ++comp) for (d = 0, m[comp] = 0.0; d < dim; ++d) m[comp] += cofu_x[comp*dim+d] * n[d];
213:   for (compI = 0; compI < Ncomp; ++compI) {
214:     for (compJ = 0; compJ < Ncomp; ++compJ) {
215:       for (d = 0; d < dim; ++d) {
216:         g1[(compI*Ncomp+compJ)*dim+d] = p * (m[compI] * cofu_x[compJ*dim+d] - cofu_x[compI*dim+d] * m[compJ]);
217:       }
218:     }
219:   }
220: }

222: void f0_p_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
223:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
224:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
225:           PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
226: {
227:   PetscReal detu_x;
228:   Det3D(&detu_x, u_x);
229:   f0[0] = detu_x - 1.0;
230: }

232: void g1_pu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
233:            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
234:            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
235:            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
236: {
237:   PetscReal cofu_x[9/*Ncomp*dim*/];
238:   PetscInt  compI, d;

240:   Cof3D(cofu_x, u_x);
241:   for (compI = 0; compI < dim; ++compI)
242:     for (d = 0; d < dim; ++d) g1[compI*dim+d] = cofu_x[compI*dim+d];
243: }

245: void g2_up_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
246:            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
247:            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
248:            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
249: {
250:   PetscReal cofu_x[9/*Ncomp*dim*/];
251:   PetscInt  compI, d;

253:   Cof3D(cofu_x, u_x);
254:   for (compI = 0; compI < dim; ++compI)
255:     for (d = 0; d < dim; ++d) g2[compI*dim+d] = cofu_x[compI*dim+d];
256: }

258: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
259: {
260:   const char    *runTypes[2] = {"full", "test"};
261:   PetscInt       run;

265:   options->runType      = RUN_FULL;
266:   options->mu           = 1.0;
267:   options->p_wall       = 0.4;

269:   PetscOptionsBegin(comm, "", "Nonlinear elasticity problem options", "DMPLEX");
270:   run  = options->runType;
271:   PetscOptionsEList("-run_type", "The run type", "ex77.c", runTypes, 2, runTypes[options->runType], &run, NULL);
272:   options->runType = (RunType) run;
273:   PetscOptionsReal("-shear_modulus", "The shear modulus", "ex77.c", options->mu, &options->mu, NULL);
274:   PetscOptionsReal("-wall_pressure", "The wall pressure", "ex77.c", options->p_wall, &options->p_wall, NULL);
275:   PetscOptionsEnd();
276:   return 0;
277: }

279: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
280: {
282:   /* TODO The P1 coordinate space gives wrong results when compared to the affine version. Track this down */
283:   if (0) {
284:     DMPlexCreateBoxMesh(comm, 3, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);
285:   } else {
286:     DMCreate(comm, dm);
287:     DMSetType(*dm, DMPLEX);
288:   }
289:   DMSetFromOptions(*dm);
290:   /* Label the faces (bit of a hack here, until it is properly implemented for simplices) */
291:   DMViewFromOptions(*dm, NULL, "-orig_dm_view");
292:   {
293:     DM              cdm;
294:     DMLabel         label;
295:     IS              is;
296:     PetscInt        d, dim, b, f, Nf;
297:     const PetscInt *faces;
298:     PetscInt        csize;
299:     PetscScalar    *coords = NULL;
300:     PetscSection    cs;
301:     Vec             coordinates ;

303:     DMGetDimension(*dm, &dim);
304:     DMCreateLabel(*dm, "boundary");
305:     DMGetLabel(*dm, "boundary", &label);
306:     DMPlexMarkBoundaryFaces(*dm, 1, label);
307:     DMGetStratumIS(*dm, "boundary", 1,  &is);
308:     if (is) {
309:       PetscReal faceCoord;
310:       PetscInt  v;

312:       ISGetLocalSize(is, &Nf);
313:       ISGetIndices(is, &faces);

315:       DMGetCoordinatesLocal(*dm, &coordinates);
316:       DMGetCoordinateDM(*dm, &cdm);
317:       DMGetLocalSection(cdm, &cs);

319:       /* Check for each boundary face if any component of its centroid is either 0.0 or 1.0 */
320:       for (f = 0; f < Nf; ++f) {
321:         DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
322:         /* Calculate mean coordinate vector */
323:         for (d = 0; d < dim; ++d) {
324:           const PetscInt Nv = csize/dim;
325:           faceCoord = 0.0;
326:           for (v = 0; v < Nv; ++v) faceCoord += PetscRealPart(coords[v*dim+d]);
327:           faceCoord /= Nv;
328:           for (b = 0; b < 2; ++b) {
329:             if (PetscAbs(faceCoord - b*1.0) < PETSC_SMALL) {
330:               DMSetLabelValue(*dm, "Faces", faces[f], d*2+b+1);
331:             }
332:           }
333:         }
334:         DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
335:       }
336:       ISRestoreIndices(is, &faces);
337:     }
338:     ISDestroy(&is);
339:   }
340:   DMViewFromOptions(*dm, NULL, "-dm_view");
341:   return 0;
342: }

344: PetscErrorCode SetupProblem(DM dm, PetscInt dim, AppCtx *user)
345: {
346:   PetscDS        ds;
347:   PetscWeakForm  wf;
348:   DMLabel        label;
349:   PetscInt       bd;

352:   DMGetDS(dm, &ds);
353:   PetscDSSetResidual(ds, 0, NULL, f1_u_3d);
354:   PetscDSSetResidual(ds, 1, f0_p_3d, NULL);
355:   PetscDSSetJacobian(ds, 0, 0, NULL, NULL,  NULL,  g3_uu_3d);
356:   PetscDSSetJacobian(ds, 0, 1, NULL, NULL,  g2_up_3d, NULL);
357:   PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu_3d, NULL,  NULL);

359:   DMGetLabel(dm, "Faces", &label);
360:   DMAddBoundary(dm, DM_BC_NATURAL, "pressure", label, 0, NULL, 0, 0, NULL, NULL, NULL, user, &bd);
361:   PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
362:   PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_u_3d, 0, NULL);
363:   PetscWeakFormSetIndexBdJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, g1_bd_uu_3d, 0, NULL, 0, NULL);

365:   DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed", label, 0, NULL, 0, 0, NULL, (void (*)(void)) coordinates, NULL, user, NULL);
366:   return 0;
367: }

369: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
370: {
371:   PetscErrorCode (*matFuncs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {elasticityMaterial, wallPressure};
372:   Vec            A;
373:   void          *ctxs[2];

375:   ctxs[0] = user; ctxs[1] = user;
376:   DMCreateLocalVector(dmAux, &A);
377:   DMProjectFunctionLocal(dmAux, 0.0, matFuncs, ctxs, INSERT_ALL_VALUES, A);
378:   DMSetAuxiliaryVec(dm, NULL, 0, 0, A);
379:   VecDestroy(&A);
380:   return 0;
381: }

383: PetscErrorCode SetupNearNullSpace(DM dm, AppCtx *user)
384: {
385:   /* Set up the near null space (a.k.a. rigid body modes) that will be used by the multigrid preconditioner */
386:   DM             subdm;
387:   MatNullSpace   nearNullSpace;
388:   PetscInt       fields = 0;
389:   PetscObject    deformation;

392:   DMCreateSubDM(dm, 1, &fields, NULL, &subdm);
393:   DMPlexCreateRigidBody(subdm, 0, &nearNullSpace);
394:   DMGetField(dm, 0, NULL, &deformation);
395:   PetscObjectCompose(deformation, "nearnullspace", (PetscObject) nearNullSpace);
396:   DMDestroy(&subdm);
397:   MatNullSpaceDestroy(&nearNullSpace);
398:   return 0;
399: }

401: static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
402: {
403:   DM             dmAux, coordDM;
404:   PetscInt       f;

406:   /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
407:   DMGetCoordinateDM(dm, &coordDM);
408:   if (!feAux) return 0;
409:   DMClone(dm, &dmAux);
410:   DMSetCoordinateDM(dmAux, coordDM);
411:   for (f = 0; f < NfAux; ++f) DMSetField(dmAux, f, NULL, (PetscObject) feAux[f]);
412:   DMCreateDS(dmAux);
413:   SetupMaterial(dm, dmAux, user);
414:   DMDestroy(&dmAux);
415:   return 0;
416: }

418: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
419: {
420:   DM              cdm = dm;
421:   PetscFE         fe[2], feAux[2];
422:   PetscBool       simplex;
423:   PetscInt        dim;
424:   MPI_Comm        comm;

427:   PetscObjectGetComm((PetscObject) dm, &comm);
428:   DMGetDimension(dm, &dim);
429:   DMPlexIsSimplex(dm, &simplex);
430:   /* Create finite element */
431:   PetscFECreateDefault(comm, dim, dim, simplex, "def_", PETSC_DEFAULT, &fe[0]);
432:   PetscObjectSetName((PetscObject) fe[0], "deformation");
433:   PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);
434:   PetscFECopyQuadrature(fe[0], fe[1]);

436:   PetscObjectSetName((PetscObject) fe[1], "pressure");

438:   PetscFECreateDefault(comm, dim, 1, simplex, "elastMat_", PETSC_DEFAULT, &feAux[0]);
439:   PetscObjectSetName((PetscObject) feAux[0], "elasticityMaterial");
440:   PetscFECopyQuadrature(fe[0], feAux[0]);
441:   /* It is not yet possible to define a field on a submesh (e.g. a boundary), so we will use a normal finite element */
442:   PetscFECreateDefault(comm, dim, 1, simplex, "wall_pres_", PETSC_DEFAULT, &feAux[1]);
443:   PetscObjectSetName((PetscObject) feAux[1], "wall_pressure");
444:   PetscFECopyQuadrature(fe[0], feAux[1]);

446:   /* Set discretization and boundary conditions for each mesh */
447:   DMSetField(dm, 0, NULL, (PetscObject) fe[0]);
448:   DMSetField(dm, 1, NULL, (PetscObject) fe[1]);
449:   DMCreateDS(dm);
450:   SetupProblem(dm, dim, user);
451:   while (cdm) {
452:     SetupAuxDM(cdm, 2, feAux, user);
453:     DMCopyDisc(dm, cdm);
454:     DMGetCoarseDM(cdm, &cdm);
455:   }
456:   PetscFEDestroy(&fe[0]);
457:   PetscFEDestroy(&fe[1]);
458:   PetscFEDestroy(&feAux[0]);
459:   PetscFEDestroy(&feAux[1]);
460:   return 0;
461: }

463: int main(int argc, char **argv)
464: {
465:   SNES           snes;                 /* nonlinear solver */
466:   DM             dm;                   /* problem definition */
467:   Vec            u,r;                  /* solution, residual vectors */
468:   Mat            A,J;                  /* Jacobian matrix */
469:   AppCtx         user;                 /* user-defined work context */
470:   PetscInt       its;                  /* iterations for convergence */

472:   PetscInitialize(&argc, &argv, NULL,help);
473:   ProcessOptions(PETSC_COMM_WORLD, &user);
474:   SNESCreate(PETSC_COMM_WORLD, &snes);
475:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
476:   SNESSetDM(snes, dm);
477:   DMSetApplicationContext(dm, &user);

479:   SetupDiscretization(dm, &user);
480:   DMPlexCreateClosureIndex(dm, NULL);
481:   SetupNearNullSpace(dm, &user);

483:   DMCreateGlobalVector(dm, &u);
484:   VecDuplicate(u, &r);

486:   DMSetMatType(dm,MATAIJ);
487:   DMCreateMatrix(dm, &J);
488:   A = J;

490:   DMPlexSetSNESLocalFEM(dm,&user,&user,&user);
491:   SNESSetJacobian(snes, A, J, NULL, NULL);

493:   SNESSetFromOptions(snes);

495:   {
496:     PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void* ctx);
497:     initialGuess[0] = coordinates; initialGuess[1] = zero_scalar;
498:     DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
499:   }

501:   if (user.runType == RUN_FULL) {
502:     SNESSolve(snes, NULL, u);
503:     SNESGetIterationNumber(snes, &its);
504:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
505:   } else {
506:     PetscReal res = 0.0;

508:     /* Check initial guess */
509:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
510:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
511:     /* Check residual */
512:     SNESComputeFunction(snes, u, r);
513:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
514:     VecChop(r, 1.0e-10);
515:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
516:     VecNorm(r, NORM_2, &res);
517:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", (double)res);
518:     /* Check Jacobian */
519:     {
520:       Vec b;

522:       SNESComputeJacobian(snes, u, A, A);
523:       VecDuplicate(u, &b);
524:       VecSet(r, 0.0);
525:       SNESComputeFunction(snes, r, b);
526:       MatMult(A, u, r);
527:       VecAXPY(r, 1.0, b);
528:       VecDestroy(&b);
529:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
530:       VecChop(r, 1.0e-10);
531:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
532:       VecNorm(r, NORM_2, &res);
533:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", (double)res);
534:     }
535:   }
536:   VecViewFromOptions(u, NULL, "-sol_vec_view");

538:   if (A != J) MatDestroy(&A);
539:   MatDestroy(&J);
540:   VecDestroy(&u);
541:   VecDestroy(&r);
542:   SNESDestroy(&snes);
543:   DMDestroy(&dm);
544:   PetscFinalize();
545:   return 0;
546: }

548: /*TEST

550:   build:
551:     requires: !complex

553:   testset:
554:     requires: ctetgen
555:     args: -run_type full -dm_plex_dim 3 \
556:           -def_petscspace_degree 2 -pres_petscspace_degree 1 -elastMat_petscspace_degree 0 -wall_pres_petscspace_degree 0 \
557:           -snes_rtol 1e-05 -snes_monitor_short -snes_converged_reason \
558:           -ksp_type fgmres -ksp_rtol 1e-10 -ksp_monitor_short -ksp_converged_reason \
559:           -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper \
560:             -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu \
561:             -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi

563:     test:
564:       suffix: 0
565:       requires: !single
566:       args: -dm_refine 2 \
567:             -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4

569:     test:
570:       suffix: 1
571:       requires: superlu_dist
572:       nsize: 2
573:       args: -dm_refine 0 -petscpartitioner_type simple \
574:             -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
575:       timeoutfactor: 2

577:     test:
578:       suffix: 4
579:       requires: superlu_dist
580:       nsize: 2
581:       args: -dm_refine 0 -petscpartitioner_type simple \
582:             -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
583:       output_file: output/ex77_1.out

585:     test:
586:       suffix: 2
587:       requires: !single
588:       args: -dm_refine 2 \
589:             -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0

591:     test:
592:       suffix: 2_par
593:       requires: superlu_dist !single
594:       nsize: 4
595:       args: -dm_refine 2 -petscpartitioner_type simple \
596:             -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0
597:       output_file: output/ex77_2.out

599: TEST*/