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 {
 45:   RUN_FULL,
 46:   RUN_TEST
 47: } RunType;

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

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

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

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

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

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

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

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

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

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

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

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

128: void f1_u_3d(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 f1[])
129: {
130:   const PetscInt  Ncomp = dim;
131:   const PetscReal mu = PetscRealPart(a[0]), kappa = 3.0;
132:   PetscReal       cofu_x[9 /* Ncomp*dim */], detu_x, p = PetscRealPart(u[Ncomp]);
133:   PetscInt        comp, d;

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

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

145: void g3_uu_3d(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 g3[])
146: {
147:   const PetscInt  Ncomp = dim;
148:   const PetscReal mu = PetscRealPart(a[0]), kappa = 3.0;
149:   PetscReal       cofu_x[9 /* Ncomp*dim */], detu_x, pp, pm, p = PetscRealPart(u[Ncomp]);
150:   PetscInt        compI, compJ, d1, d2;

152:   Cof3D(cofu_x, u_x);
153:   Det3D(&detu_x, u_x);
154:   p += kappa * (detu_x - 1.0);
155:   pp = p / detu_x + kappa;
156:   pm = p / detu_x;

158:   /* 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} */
159:   for (compI = 0; compI < Ncomp; ++compI) {
160:     for (compJ = 0; compJ < Ncomp; ++compJ) {
161:       const PetscReal G = (compI == compJ) ? 1.0 : 0.0;

163:       for (d1 = 0; d1 < dim; ++d1) {
164:         for (d2 = 0; d2 < dim; ++d2) {
165:           const PetscReal g = (d1 == d2) ? 1.0 : 0.0;

167:           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];
168:         }
169:       }
170:     }
171:   }
172: }

174: void f0_bd_u_3d(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[])
175: {
176:   const PetscInt    Ncomp = dim;
177:   const PetscScalar p     = a[aOff[1]];
178:   PetscReal         cofu_x[9 /* Ncomp*dim */];
179:   PetscInt          comp, d;

181:   Cof3D(cofu_x, u_x);
182:   for (comp = 0; comp < Ncomp; ++comp) {
183:     for (d = 0, f0[comp] = 0.0; d < dim; ++d) f0[comp] += cofu_x[comp * dim + d] * n[d];
184:     f0[comp] *= p;
185:   }
186: }

188: void g1_bd_uu_3d(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
189: {
190:   const PetscInt Ncomp = dim;
191:   PetscScalar    p     = a[aOff[1]];
192:   PetscReal      cofu_x[9 /* Ncomp*dim */], m[3 /* Ncomp */], detu_x;
193:   PetscInt       comp, compI, compJ, d;

195:   Cof3D(cofu_x, u_x);
196:   Det3D(&detu_x, u_x);
197:   p /= detu_x;

199:   for (comp = 0; comp < Ncomp; ++comp)
200:     for (d = 0, m[comp] = 0.0; d < dim; ++d) m[comp] += cofu_x[comp * dim + d] * n[d];
201:   for (compI = 0; compI < Ncomp; ++compI) {
202:     for (compJ = 0; compJ < Ncomp; ++compJ) {
203:       for (d = 0; d < dim; ++d) g1[(compI * Ncomp + compJ) * dim + d] = p * (m[compI] * cofu_x[compJ * dim + d] - cofu_x[compI * dim + d] * m[compJ]);
204:     }
205:   }
206: }

208: void f0_p_3d(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[])
209: {
210:   PetscReal detu_x;
211:   Det3D(&detu_x, u_x);
212:   f0[0] = detu_x - 1.0;
213: }

215: void g1_pu_3d(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[])
216: {
217:   PetscReal cofu_x[9 /* Ncomp*dim */];
218:   PetscInt  compI, d;

220:   Cof3D(cofu_x, u_x);
221:   for (compI = 0; compI < dim; ++compI)
222:     for (d = 0; d < dim; ++d) g1[compI * dim + d] = cofu_x[compI * dim + d];
223: }

225: void g2_up_3d(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 g2[])
226: {
227:   PetscReal cofu_x[9 /* Ncomp*dim */];
228:   PetscInt  compI, d;

230:   Cof3D(cofu_x, u_x);
231:   for (compI = 0; compI < dim; ++compI)
232:     for (d = 0; d < dim; ++d) g2[compI * dim + d] = cofu_x[compI * dim + d];
233: }

235: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
236: {
237:   const char *runTypes[2] = {"full", "test"};
238:   PetscInt    run;

241:   options->runType = RUN_FULL;
242:   options->mu      = 1.0;
243:   options->p_wall  = 0.4;
244:   PetscOptionsBegin(comm, "", "Nonlinear elasticity problem options", "DMPLEX");
245:   run = options->runType;
246:   PetscOptionsEList("-run_type", "The run type", "ex77.c", runTypes, 2, runTypes[options->runType], &run, NULL);
247:   options->runType = (RunType)run;
248:   PetscOptionsReal("-shear_modulus", "The shear modulus", "ex77.c", options->mu, &options->mu, NULL);
249:   PetscOptionsReal("-wall_pressure", "The wall pressure", "ex77.c", options->p_wall, &options->p_wall, NULL);
250:   PetscOptionsEnd();
251:   return 0;
252: }

254: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
255: {
257:   /* TODO The P1 coordinate space gives wrong results when compared to the affine version. Track this down */
258:   if (0) {
259:     DMPlexCreateBoxMesh(comm, 3, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);
260:   } else {
261:     DMCreate(comm, dm);
262:     DMSetType(*dm, DMPLEX);
263:   }
264:   DMSetFromOptions(*dm);
265:   /* Label the faces (bit of a hack here, until it is properly implemented for simplices) */
266:   DMViewFromOptions(*dm, NULL, "-orig_dm_view");
267:   {
268:     DM              cdm;
269:     DMLabel         label;
270:     IS              is;
271:     PetscInt        d, dim, b, f, Nf;
272:     const PetscInt *faces;
273:     PetscInt        csize;
274:     PetscScalar    *coords = NULL;
275:     PetscSection    cs;
276:     Vec             coordinates;

278:     DMGetDimension(*dm, &dim);
279:     DMCreateLabel(*dm, "boundary");
280:     DMGetLabel(*dm, "boundary", &label);
281:     DMPlexMarkBoundaryFaces(*dm, 1, label);
282:     DMGetStratumIS(*dm, "boundary", 1, &is);
283:     if (is) {
284:       PetscReal faceCoord;
285:       PetscInt  v;

287:       ISGetLocalSize(is, &Nf);
288:       ISGetIndices(is, &faces);

290:       DMGetCoordinatesLocal(*dm, &coordinates);
291:       DMGetCoordinateDM(*dm, &cdm);
292:       DMGetLocalSection(cdm, &cs);

294:       /* Check for each boundary face if any component of its centroid is either 0.0 or 1.0 */
295:       for (f = 0; f < Nf; ++f) {
296:         DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
297:         /* Calculate mean coordinate vector */
298:         for (d = 0; d < dim; ++d) {
299:           const PetscInt Nv = csize / dim;
300:           faceCoord         = 0.0;
301:           for (v = 0; v < Nv; ++v) faceCoord += PetscRealPart(coords[v * dim + d]);
302:           faceCoord /= Nv;
303:           for (b = 0; b < 2; ++b) {
304:             if (PetscAbs(faceCoord - b * 1.0) < PETSC_SMALL) DMSetLabelValue(*dm, "Faces", faces[f], d * 2 + b + 1);
305:           }
306:         }
307:         DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
308:       }
309:       ISRestoreIndices(is, &faces);
310:     }
311:     ISDestroy(&is);
312:   }
313:   DMViewFromOptions(*dm, NULL, "-dm_view");
314:   return 0;
315: }

317: PetscErrorCode SetupProblem(DM dm, PetscInt dim, AppCtx *user)
318: {
319:   PetscDS       ds;
320:   PetscWeakForm wf;
321:   DMLabel       label;
322:   PetscInt      bd;

325:   DMGetDS(dm, &ds);
326:   PetscDSSetResidual(ds, 0, NULL, f1_u_3d);
327:   PetscDSSetResidual(ds, 1, f0_p_3d, NULL);
328:   PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu_3d);
329:   PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up_3d, NULL);
330:   PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu_3d, NULL, NULL);

332:   DMGetLabel(dm, "Faces", &label);
333:   DMAddBoundary(dm, DM_BC_NATURAL, "pressure", label, 0, NULL, 0, 0, NULL, NULL, NULL, user, &bd);
334:   PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
335:   PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_u_3d, 0, NULL);
336:   PetscWeakFormSetIndexBdJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, g1_bd_uu_3d, 0, NULL, 0, NULL);

338:   DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed", label, 0, NULL, 0, 0, NULL, (void (*)(void))coordinates, NULL, user, NULL);
339:   return 0;
340: }

342: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
343: {
344:   PetscErrorCode (*matFuncs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {elasticityMaterial, wallPressure};
345:   Vec   A;
346:   void *ctxs[2];

348:   ctxs[0] = user;
349:   ctxs[1] = user;
350:   DMCreateLocalVector(dmAux, &A);
351:   DMProjectFunctionLocal(dmAux, 0.0, matFuncs, ctxs, INSERT_ALL_VALUES, A);
352:   DMSetAuxiliaryVec(dm, NULL, 0, 0, A);
353:   VecDestroy(&A);
354:   return 0;
355: }

357: PetscErrorCode SetupNearNullSpace(DM dm, AppCtx *user)
358: {
359:   /* Set up the near null space (a.k.a. rigid body modes) that will be used by the multigrid preconditioner */
360:   DM           subdm;
361:   MatNullSpace nearNullSpace;
362:   PetscInt     fields = 0;
363:   PetscObject  deformation;

366:   DMCreateSubDM(dm, 1, &fields, NULL, &subdm);
367:   DMPlexCreateRigidBody(subdm, 0, &nearNullSpace);
368:   DMGetField(dm, 0, NULL, &deformation);
369:   PetscObjectCompose(deformation, "nearnullspace", (PetscObject)nearNullSpace);
370:   DMDestroy(&subdm);
371:   MatNullSpaceDestroy(&nearNullSpace);
372:   return 0;
373: }

375: static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
376: {
377:   DM       dmAux, coordDM;
378:   PetscInt f;

380:   /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
381:   DMGetCoordinateDM(dm, &coordDM);
382:   if (!feAux) return 0;
383:   DMClone(dm, &dmAux);
384:   DMSetCoordinateDM(dmAux, coordDM);
385:   for (f = 0; f < NfAux; ++f) DMSetField(dmAux, f, NULL, (PetscObject)feAux[f]);
386:   DMCreateDS(dmAux);
387:   SetupMaterial(dm, dmAux, user);
388:   DMDestroy(&dmAux);
389:   return 0;
390: }

392: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
393: {
394:   DM        cdm = dm;
395:   PetscFE   fe[2], feAux[2];
396:   PetscBool simplex;
397:   PetscInt  dim;
398:   MPI_Comm  comm;

401:   PetscObjectGetComm((PetscObject)dm, &comm);
402:   DMGetDimension(dm, &dim);
403:   DMPlexIsSimplex(dm, &simplex);
404:   /* Create finite element */
405:   PetscFECreateDefault(comm, dim, dim, simplex, "def_", PETSC_DEFAULT, &fe[0]);
406:   PetscObjectSetName((PetscObject)fe[0], "deformation");
407:   PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);
408:   PetscFECopyQuadrature(fe[0], fe[1]);

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

412:   PetscFECreateDefault(comm, dim, 1, simplex, "elastMat_", PETSC_DEFAULT, &feAux[0]);
413:   PetscObjectSetName((PetscObject)feAux[0], "elasticityMaterial");
414:   PetscFECopyQuadrature(fe[0], feAux[0]);
415:   /* It is not yet possible to define a field on a submesh (e.g. a boundary), so we will use a normal finite element */
416:   PetscFECreateDefault(comm, dim, 1, simplex, "wall_pres_", PETSC_DEFAULT, &feAux[1]);
417:   PetscObjectSetName((PetscObject)feAux[1], "wall_pressure");
418:   PetscFECopyQuadrature(fe[0], feAux[1]);

420:   /* Set discretization and boundary conditions for each mesh */
421:   DMSetField(dm, 0, NULL, (PetscObject)fe[0]);
422:   DMSetField(dm, 1, NULL, (PetscObject)fe[1]);
423:   DMCreateDS(dm);
424:   SetupProblem(dm, dim, user);
425:   while (cdm) {
426:     SetupAuxDM(cdm, 2, feAux, user);
427:     DMCopyDisc(dm, cdm);
428:     DMGetCoarseDM(cdm, &cdm);
429:   }
430:   PetscFEDestroy(&fe[0]);
431:   PetscFEDestroy(&fe[1]);
432:   PetscFEDestroy(&feAux[0]);
433:   PetscFEDestroy(&feAux[1]);
434:   return 0;
435: }

437: int main(int argc, char **argv)
438: {
439:   SNES     snes; /* nonlinear solver */
440:   DM       dm;   /* problem definition */
441:   Vec      u, r; /* solution, residual vectors */
442:   Mat      A, J; /* Jacobian matrix */
443:   AppCtx   user; /* user-defined work context */
444:   PetscInt its;  /* iterations for convergence */

447:   PetscInitialize(&argc, &argv, NULL, help);
448:   ProcessOptions(PETSC_COMM_WORLD, &user);
449:   SNESCreate(PETSC_COMM_WORLD, &snes);
450:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
451:   SNESSetDM(snes, dm);
452:   DMSetApplicationContext(dm, &user);

454:   SetupDiscretization(dm, &user);
455:   DMPlexCreateClosureIndex(dm, NULL);
456:   SetupNearNullSpace(dm, &user);

458:   DMCreateGlobalVector(dm, &u);
459:   VecDuplicate(u, &r);

461:   DMSetMatType(dm, MATAIJ);
462:   DMCreateMatrix(dm, &J);
463:   A = J;

465:   DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
466:   SNESSetJacobian(snes, A, J, NULL, NULL);

468:   SNESSetFromOptions(snes);

470:   {
471:     PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
472:     initialGuess[0] = coordinates;
473:     initialGuess[1] = zero_scalar;
474:     DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
475:   }

477:   if (user.runType == RUN_FULL) {
478:     SNESSolve(snes, NULL, u);
479:     SNESGetIterationNumber(snes, &its);
480:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its);
481:   } else {
482:     PetscReal res = 0.0;

484:     /* Check initial guess */
485:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
486:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
487:     /* Check residual */
488:     SNESComputeFunction(snes, u, r);
489:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
490:     VecChop(r, 1.0e-10);
491:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
492:     VecNorm(r, NORM_2, &res);
493:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", (double)res);
494:     /* Check Jacobian */
495:     {
496:       Vec b;

498:       SNESComputeJacobian(snes, u, A, A);
499:       VecDuplicate(u, &b);
500:       VecSet(r, 0.0);
501:       SNESComputeFunction(snes, r, b);
502:       MatMult(A, u, r);
503:       VecAXPY(r, 1.0, b);
504:       VecDestroy(&b);
505:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
506:       VecChop(r, 1.0e-10);
507:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
508:       VecNorm(r, NORM_2, &res);
509:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", (double)res);
510:     }
511:   }
512:   VecViewFromOptions(u, NULL, "-sol_vec_view");

514:   if (A != J) MatDestroy(&A);
515:   MatDestroy(&J);
516:   VecDestroy(&u);
517:   VecDestroy(&r);
518:   SNESDestroy(&snes);
519:   DMDestroy(&dm);
520:   PetscFinalize();
521:   return 0;
522: }

524: /*TEST

526:   build:
527:     requires: !complex

529:   testset:
530:     requires: ctetgen
531:     args: -run_type full -dm_plex_dim 3 \
532:           -def_petscspace_degree 2 -pres_petscspace_degree 1 -elastMat_petscspace_degree 0 -wall_pres_petscspace_degree 0 \
533:           -snes_rtol 1e-05 -snes_monitor_short -snes_converged_reason \
534:           -ksp_type fgmres -ksp_rtol 1e-10 -ksp_monitor_short -ksp_converged_reason \
535:           -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper \
536:             -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu \
537:             -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi

539:     test:
540:       suffix: 0
541:       requires: !single
542:       args: -dm_refine 2 \
543:             -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4

545:     test:
546:       suffix: 1
547:       requires: superlu_dist
548:       nsize: 2
549:       args: -dm_refine 0 -petscpartitioner_type simple \
550:             -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
551:       timeoutfactor: 2

553:     test:
554:       suffix: 4
555:       requires: superlu_dist
556:       nsize: 2
557:       args: -dm_refine 0 -petscpartitioner_type simple \
558:             -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
559:       output_file: output/ex77_1.out

561:     test:
562:       suffix: 2
563:       requires: !single
564:       args: -dm_refine 2 \
565:             -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0

567:     test:
568:       suffix: 2_par
569:       requires: superlu_dist !single
570:       nsize: 4
571:       args: -dm_refine 2 -petscpartitioner_type simple \
572:             -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0
573:       output_file: output/ex77_2.out

575: TEST*/