Actual source code: ex77.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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:   PetscInt      debug;             /* The debugging level */
 48:   RunType       runType;           /* Whether to run tests, or solve the full problem */
 49:   PetscLogEvent createMeshEvent;
 50:   PetscBool     showInitial, showSolution;
 51:   /* Domain and mesh definition */
 52:   PetscInt      dim;               /* The topological mesh dimension */
 53:   PetscBool     interpolate;       /* Generate intermediate mesh elements */
 54:   PetscBool     simplex;           /* Use simplices or tensor product cells */
 55:   PetscReal     refinementLimit;   /* The largest allowable cell volume */
 56:   PetscBool     testPartition;     /* Use a fixed partitioning for testing */
 57:   PetscReal     mu;                /* The shear modulus */
 58: } AppCtx;

 60: /* Kronecker-delta */
 61: static const PetscInt delta2D[2*2] = {1,0,0,1};
 62: static const PetscInt delta3D[3*3] = {1,0,0,0,1,0,0,0,1};

 64: /* Levi-Civita symbol */
 65: static const PetscInt epsilon2D[2*2] = {0,1,-1,0};
 66: static const PetscInt epsilon3D[3*3*3] = {0,0,0,0,0,1,0,-1,0,0,0,-1,0,0,0,1,0,0,0,1,0,-1,0,0,0,0,0};

 68: PETSC_STATIC_INLINE void Det2D(PetscReal *detJ, const PetscReal J[])
 69: {
 70:   *detJ = J[0]*J[3] - J[1]*J[2];
 71: }

 73: PETSC_STATIC_INLINE void Det3D(PetscReal *detJ, const PetscReal J[])
 74: {
 75:   *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
 76:            J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
 77:            J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
 78: }

 80: PETSC_STATIC_INLINE void Cof2D(PetscReal C[], const PetscReal A[])
 81: {
 82:   C[0] =  A[3];
 83:   C[1] = -A[2];
 84:   C[2] = -A[1];
 85:   C[3] =  A[0];
 86: }

 88: PETSC_STATIC_INLINE void Cof3D(PetscReal C[], const PetscReal A[])
 89: {
 90:   C[0*3+0] = A[1*3+1]*A[2*3+2] - A[1*3+2]*A[2*3+1];
 91:   C[0*3+1] = A[1*3+2]*A[2*3+0] - A[1*3+0]*A[2*3+2];
 92:   C[0*3+2] = A[1*3+0]*A[2*3+1] - A[1*3+1]*A[2*3+0];
 93:   C[1*3+0] = A[0*3+2]*A[2*3+1] - A[0*3+1]*A[2*3+2];
 94:   C[1*3+1] = A[0*3+0]*A[2*3+2] - A[0*3+2]*A[2*3+0];
 95:   C[1*3+2] = A[0*3+1]*A[2*3+0] - A[0*3+0]*A[2*3+1];
 96:   C[2*3+0] = A[0*3+1]*A[1*3+2] - A[0*3+2]*A[1*3+1];
 97:   C[2*3+1] = A[0*3+2]*A[1*3+0] - A[0*3+0]*A[1*3+2];
 98:   C[2*3+2] = A[0*3+0]*A[1*3+1] - A[0*3+1]*A[1*3+0];
 99: }

101: PetscErrorCode zero_scalar(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
102: {
103:   u[0] = 0.0;
104:   return 0;
105: }

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

111:   PetscInt       comp;
112:   for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0.0;
113:   return 0;
114: }

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

120:   PetscInt       comp;
121:   for (comp = 0; comp < Ncomp; ++comp) u[comp] = x[comp];
122:   return 0;
123: }

125: PetscErrorCode elasticityMaterial(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
126: {
127:   AppCtx *user = (AppCtx *) ctx;
128:   u[0] = user->mu;
129:   return 0;
130: }

132: void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
133:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
134:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
135:           PetscReal t, const PetscReal x[], PetscScalar f0[])
136: {
137:   const PetscInt Ncomp = dim;

139:   PetscInt comp;
140:   for (comp = 0; comp < Ncomp; ++comp) f0[comp] = 0.0;
141: }

143: void f1_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
144:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
145:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
146:           PetscReal t, const PetscReal x[], PetscScalar f1[])
147: {
148:   const PetscInt  Ncomp = dim;
149:   const PetscReal mu = a[0], p = u[Ncomp], kappa = 3.0;
150:   PetscReal       cofu_x[Ncomp*dim], detu_x;

152:   Cof3D(cofu_x, u_x);
153:   Det3D(&detu_x, u_x);

155:   /* f1 is the first Piola-Kirchhoff tensor */
156:   PetscInt        comp, d;
157:   for (comp = 0; comp < Ncomp; ++comp) {
158:     for (d = 0; d < dim; ++d) {
159:       f1[comp*dim+d] = mu * u_x[comp*dim+d] + cofu_x[comp*dim+d] * (p + kappa * (detu_x - 1.0));
160:     }
161:   }
162: }

164: void g3_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
165:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
166:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
167:           PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
168: {
169:   const PetscInt  Ncomp = dim;
170:   const PetscReal mu = a[0], p = u[Ncomp], kappa = 3.0;
171:   PetscReal       cofu_x[Ncomp*dim], detu_x;

173:   Cof3D(cofu_x, u_x);
174:   Det3D(&detu_x, u_x);

176:   /* 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} */
177:   PetscInt compI, compJ, compK, d1, d2, d3;
178:   for (compI = 0; compI < Ncomp; ++compI) {
179:     for (compJ = 0; compJ < Ncomp; ++compJ) {
180:       for (d1 = 0; d1 < dim; ++d1) {
181:         for (d2 = 0; d2 < dim; ++d2) {
182:           for (compK = 0; compK < Ncomp; ++compK) {
183:             for (d3 = 0; d3 < dim; ++d3) {
184:               g3[((compI*Ncomp+compJ)*dim+d1)*dim+d2] += (p + kappa * (detu_x - 1.0)) * epsilon3D[(compI*Ncomp+compJ)*Ncomp+compK] * epsilon3D[(d1*dim+d2)*dim+d3] * u_x[compK*dim+d3];
185:             }
186:           }
187:           g3[((compI*Ncomp+compJ)*dim+d1)*dim+d2] += kappa * cofu_x[compI*dim+d1] * cofu_x[compJ*dim+d2] + delta3D[compI*Ncomp+compJ] * delta3D[d1*dim+d2] * mu;
188:         }
189:       }
190:     }
191:   }
192: }

194: void f0_bd_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
195:     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
196:     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
197:     PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
198: {
199:   const PetscInt  Ncomp = dim;
200:   const PetscReal p = 0.4; /* Should this be specified as an auxiliary field? */
201:   PetscReal       cofu_x[Ncomp*dim];

203:   Cof3D(cofu_x, u_x);

205:   PetscInt comp, d;
206:   for (comp = 0; comp < Ncomp; ++comp) {
207:     for (d = 0, f0[comp] = 0.0; d < dim; ++d) f0[comp] += cofu_x[comp*dim+d] * n[d];
208:     f0[comp] *= p;
209:   }
210: }

212: void f1_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
213:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
214:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
215:                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[])
216: {
217:   const PetscInt Ncomp = dim;

219:   PetscInt       comp, d;
220:   for (comp = 0; comp < Ncomp; ++comp) {
221:     for (d = 0; d < dim; ++d) {
222:       f1[comp*dim+d] = 0.0;
223:     }
224:   }
225: }

227: void g1_bd_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
228:     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
229:     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
230:     PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscScalar g1[])
231: {
232:   const PetscInt  Ncomp = dim;
233:   const PetscReal p = 0.4; /* Should this be specified as an auxiliary field? */

235:   PetscInt compI, compJ, compK, d1, d2, d3;
236:   for (compI = 0; compI < Ncomp; ++compI) {
237:     for (compJ = 0; compJ < Ncomp; ++compJ) {
238:       for (d2 = 0; d2 < dim; ++d2) {
239:         for (compK = 0; compK < Ncomp; ++compK) {
240:           for (d3 = 0; d3 < dim; ++d3) {
241:             for (d1 = 0; d1 < dim; ++d1) {
242:               g1[(compI*Ncomp+compJ)*dim+d2] += epsilon3D[(compI*Ncomp+compJ)*Ncomp+compK] * epsilon3D[(d2*dim+d3)*dim+d1] * u_x[compK*dim+d3] * n[d1];
243:             }
244:           }
245:         }
246:         g1[(compI*Ncomp+compJ)*dim+d2] *= p;
247:       }
248:     }
249:   }
250: }

252: void f0_p_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
253:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
254:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
255:           PetscReal t, const PetscReal x[], PetscScalar f0[])
256: {
257:   PetscReal detu_x;
258:   Det3D(&detu_x, u_x);
259:   f0[0] = detu_x - 1.0;
260: }

262: void f1_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
263:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
264:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
265:           PetscReal t, const PetscReal x[], PetscScalar f1[])
266: {
267:   PetscInt d;
268:   for (d = 0; d < dim; ++d) f1[d] = 0.0;
269: }

271: void f0_bd_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
272:     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
273:     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
274:     PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
275: {
276:   f0[0] = 0.0;
277: }
278: void f1_bd_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
279:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
280:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
281:                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[])
282: {
283:   PetscInt d;
284:   for (d = 0; d < dim; ++d) f1[d] = 0.0;
285: }

287: void g1_pu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
288:            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
289:            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
290:            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g1[])
291: {
292:   Cof3D(g1, u_x);
293: }

295: void g2_up_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
296:            const PetscInt uOff[], const PetscInt uOff_x[], const PetscReal u[], const PetscReal u_t[], const PetscReal u_x[],
297:            const PetscInt aOff[], const PetscInt aOff_x[], const PetscReal a[], const PetscReal a_t[], const PetscReal a_x[],
298:            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscReal g2[])
299: {
300:   Cof3D(g2, u_x);
301: }

305: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
306: {
307:   const char    *runTypes[2] = {"full", "test"};
308:   PetscInt       run;

312:   options->debug           = 0;
313:   options->runType         = RUN_FULL;
314:   options->dim             = 3;
315:   options->interpolate     = PETSC_FALSE;
316:   options->simplex         = PETSC_TRUE;
317:   options->refinementLimit = 0.0;
318:   options->mu              = 1.0;
319:   options->testPartition   = PETSC_FALSE;
320:   options->showInitial     = PETSC_FALSE;
321:   options->showSolution    = PETSC_TRUE;

323:   PetscOptionsBegin(comm, "", "Nonlinear elasticity problem options", "DMPLEX");
324:   PetscOptionsInt("-debug", "The debugging level", "ex77.c", options->debug, &options->debug, NULL);
325:   run  = options->runType;
326:   PetscOptionsEList("-run_type", "The run type", "ex77.c", runTypes, 2, runTypes[options->runType], &run, NULL);

328:   options->runType = (RunType) run;

330:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex77.c", options->dim, &options->dim, NULL);
331:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex77.c", options->interpolate, &options->interpolate, NULL);
332:   PetscOptionsBool("-simplex", "Use simplices or tensor product cells", "ex77.c", options->simplex, &options->simplex, NULL);
333:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex77.c", options->refinementLimit, &options->refinementLimit, NULL);
334:   PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex77.c", options->testPartition, &options->testPartition, NULL);
335:   PetscOptionsReal("-shear_modulus", "The shear modulus", "ex77.c", options->mu, &options->mu, NULL);

337:   PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex77.c", options->showInitial, &options->showInitial, NULL);
338:   PetscOptionsBool("-show_solution", "Output the solution for verification", "ex77.c", options->showSolution, &options->showSolution, NULL);
339:   PetscOptionsEnd();

341:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
342:   return(0);
343: }

347: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
348: {
349:   Vec            lv;
350:   PetscInt       p;
351:   PetscMPIInt    rank, numProcs;

355:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
356:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
357:   DMGetLocalVector(dm, &lv);
358:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
359:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
360:   PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
361:   for (p = 0; p < numProcs; ++p) {
362:     if (p == rank) {VecView(lv, PETSC_VIEWER_STDOUT_SELF);}
363:     PetscBarrier((PetscObject) dm);
364:   }
365:   DMRestoreLocalVector(dm, &lv);
366:   return(0);
367: }

371: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
372: {
373:   PetscInt       dim             = user->dim;
374:   PetscBool      interpolate     = user->interpolate;
375:   PetscReal      refinementLimit = user->refinementLimit;
376:   const PetscInt cells[3]        = {3, 3, 3};

380:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
381:   if (user->simplex) {DMPlexCreateBoxMesh(comm, dim, interpolate, dm);}
382:   else               {DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);}
383:   /* Label the faces (bit of a hack here, until it is properly implemented for simplices) */
384:   {
385:     DM              cdm;
386:     DMLabel         label;
387:     IS              is;
388:     PetscInt        d, dim = user->dim, b, f, Nf;
389:     const PetscInt *faces;
390:     PetscInt        csize;
391:     PetscReal      *coords = NULL;
392:     PetscSection    cs;
393:     Vec             coordinates ;

395:     DMCreateLabel(*dm, "boundary");
396:     DMGetLabel(*dm, "boundary", &label);
397:     DMPlexMarkBoundaryFaces(*dm, label);
398:     DMGetStratumIS(*dm, "boundary", 1,  &is);
399:     DMCreateLabel(*dm, "Faces");
400:     if (is) {
401:       ISGetLocalSize(is, &Nf);
402:       ISGetIndices(is, &faces);

404:       DMGetCoordinatesLocal(*dm, &coordinates);
405:       DMGetCoordinateDM(*dm, &cdm);
406:       DMGetDefaultSection(cdm, &cs);

408:       /* Check for each boundary face if any component of its centroid is either 0.0 or 1.0 */
409:       PetscReal faceCoord;
410:       PetscInt  v;
411:       for (f = 0; f < Nf; ++f) {
412:         DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
413:         /* Calculate mean coordinate vector */
414:         const PetscInt Nv = csize/dim;
415:         for (d = 0; d < dim; ++d) {
416:           faceCoord = 0.0;
417:           for (v = 0; v < Nv; ++v) faceCoord += coords[v*dim+d];
418:           faceCoord /= Nv;
419:           for (b = 0; b < 2; ++b) {
420:             if (PetscAbs(faceCoord - b*1.0) < PETSC_SMALL) {
421:               DMSetLabelValue(*dm, "Faces", faces[f], d*2+b+1);
422:             }
423:           }
424:         }
425:         DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
426:       }
427:       ISRestoreIndices(is, &faces);
428:     }
429:     ISDestroy(&is);
430:     DMGetLabel(*dm, "Faces", &label);
431:     DMPlexLabelComplete(*dm, label);
432:   }
433:   {
434:     DM refinedMesh     = NULL;
435:     DM distributedMesh = NULL;

437:     /* Refine mesh using a volume constraint */
438:     DMPlexSetRefinementLimit(*dm, refinementLimit);
439:     if (user->simplex) {DMRefine(*dm, comm, &refinedMesh);}
440:     if (refinedMesh) {
441:       DMDestroy(dm);
442:       *dm  = refinedMesh;
443:     }
444:     /* Setup test partitioning */
445:     if (user->testPartition) {
446:       PetscInt         triSizes_n2[2]       = {4, 4};
447:       PetscInt         triPoints_n2[8]      = {3, 5, 6, 7, 0, 1, 2, 4};
448:       PetscInt         triSizes_n3[3]       = {2, 3, 3};
449:       PetscInt         triPoints_n3[8]      = {3, 5, 1, 6, 7, 0, 2, 4};
450:       PetscInt         triSizes_n5[5]       = {1, 2, 2, 1, 2};
451:       PetscInt         triPoints_n5[8]      = {3, 5, 6, 4, 7, 0, 1, 2};
452:       PetscInt         triSizes_ref_n2[2]   = {8, 8};
453:       PetscInt         triPoints_ref_n2[16] = {1, 5, 6, 7, 10, 11, 14, 15, 0, 2, 3, 4, 8, 9, 12, 13};
454:       PetscInt         triSizes_ref_n3[3]   = {5, 6, 5};
455:       PetscInt         triPoints_ref_n3[16] = {1, 7, 10, 14, 15, 2, 6, 8, 11, 12, 13, 0, 3, 4, 5, 9};
456:       PetscInt         triSizes_ref_n5[5]   = {3, 4, 3, 3, 3};
457:       PetscInt         triPoints_ref_n5[16] = {1, 7, 10, 2, 11, 13, 14, 5, 6, 15, 0, 8, 9, 3, 4, 12};
458:       const PetscInt  *sizes = NULL;
459:       const PetscInt  *points = NULL;
460:       PetscPartitioner part;
461:       PetscInt         cEnd;
462:       PetscMPIInt      rank, numProcs;

464:       MPI_Comm_rank(comm, &rank);
465:       MPI_Comm_size(comm, &numProcs);
466:       DMPlexGetHeightStratum(*dm, 0, NULL, &cEnd);
467:       if (!rank) {
468:         if (dim == 2 && user->simplex && numProcs == 2 && cEnd == 8) {
469:            sizes = triSizes_n2; points = triPoints_n2;
470:         } else if (dim == 2 && user->simplex && numProcs == 3 && cEnd == 8) {
471:           sizes = triSizes_n3; points = triPoints_n3;
472:         } else if (dim == 2 && user->simplex && numProcs == 5 && cEnd == 8) {
473:           sizes = triSizes_n5; points = triPoints_n5;
474:         } else if (dim == 2 && user->simplex && numProcs == 2 && cEnd == 16) {
475:            sizes = triSizes_ref_n2; points = triPoints_ref_n2;
476:         } else if (dim == 2 && user->simplex && numProcs == 3 && cEnd == 16) {
477:           sizes = triSizes_ref_n3; points = triPoints_ref_n3;
478:         } else if (dim == 2 && user->simplex && numProcs == 5 && cEnd == 16) {
479:           sizes = triSizes_ref_n5; points = triPoints_ref_n5;
480:         } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "No stored partition matching run parameters");
481:       }
482:       DMPlexGetPartitioner(*dm, &part);
483:       PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);
484:       PetscPartitionerShellSetPartition(part, numProcs, sizes, points);
485:     }
486:     /* Distribute mesh over processes */
487:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
488:     if (distributedMesh) {
489:       DMDestroy(dm);
490:       *dm  = distributedMesh;
491:     }
492:   }
493:   DMSetFromOptions(*dm);
494:   DMViewFromOptions(*dm, NULL, "-dm_view");
495: 
496:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
497:   return(0);
498: }

502: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
503: {
504:   PetscDS        prob;

508:   DMGetDS(dm, &prob);
509:   PetscDSSetResidual(prob, 0, f0_u, f1_u_3d);
510:   PetscDSSetResidual(prob, 1, f0_p_3d, f1_p);
511:   PetscDSSetJacobian(prob, 0, 0, NULL, NULL,  NULL,  g3_uu_3d);
512:   PetscDSSetJacobian(prob, 0, 1, NULL, NULL,  g2_up_3d, NULL);
513:   PetscDSSetJacobian(prob, 1, 0, NULL, g1_pu_3d, NULL,  NULL);

515:   PetscDSSetBdResidual(prob, 0, f0_bd_u_3d, f1_bd_u);
516:   PetscDSSetBdResidual(prob, 1, f0_bd_p, f1_bd_p);
517:   PetscDSSetBdJacobian(prob, 0, 0, NULL, g1_bd_uu_3d, NULL, NULL);
518:   return(0);
519: }

523: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
524: {
525:   PetscErrorCode (*matFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {elasticityMaterial};
526:   Vec            nu;
527:   void *ctxs[] = {user};

531:   DMCreateLocalVector(dmAux, &nu);
532:   DMProjectFunctionLocal(dmAux, 0.0, matFuncs, ctxs, INSERT_ALL_VALUES, nu);
533:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
534:   VecDestroy(&nu);
535:   return(0);
536: }

540: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
541: {
542:   DM              cdm   = dm, coordDM, dmAux;
543:   const PetscInt  dim   = user->dim;
544:   const PetscBool simplex = user->simplex;
545:   PetscFE         fe[2], feBd[2], feAux;
546:   PetscQuadrature q;
547:   PetscDS         prob, probAux;
548:   PetscInt        order;
549:   PetscErrorCode  ierr;

552:   /* Create finite element */
553:   PetscFECreateDefault(dm, dim, dim, simplex, "def_", -1, &fe[0]);
554:   PetscObjectSetName((PetscObject) fe[0], "deformation");
555:   PetscFEGetQuadrature(fe[0], &q);
556:   PetscQuadratureGetOrder(q, &order);
557:   PetscFECreateDefault(dm, dim, 1, simplex, "pres_", order, &fe[1]);
558:   PetscObjectSetName((PetscObject) fe[1], "pressure");
559:   PetscFECreateDefault(dm, dim-1, dim, simplex, "bd_def_", order, &feBd[0]);
560:   PetscObjectSetName((PetscObject) feBd[0], "deformation");
561:   PetscFECreateDefault(dm, dim-1, 1, simplex, "bd_pres_", order, &feBd[1]);
562:   PetscObjectSetName((PetscObject) feBd[1], "pressure");

564:   PetscFECreateDefault(dm, dim, 1, simplex, "elastMat_", order, &feAux);
565:   PetscObjectSetName((PetscObject) feAux, "elasticityMaterial");
566: 

568:   /* Set discretization and boundary conditions for each mesh */
569:   while (cdm) {
570:     DMGetDS(cdm, &prob);
571:     PetscDSSetDiscretization(prob, 0, (PetscObject) fe[0]);
572:     PetscDSSetDiscretization(prob, 1, (PetscObject) fe[1]);
573:     PetscDSSetBdDiscretization(prob, 0, (PetscObject) feBd[0]);
574:     PetscDSSetBdDiscretization(prob, 1, (PetscObject) feBd[1]);
575:     SetupProblem(cdm, user);

577:     DMClone(cdm, &dmAux);
578:     DMGetCoordinateDM(cdm, &coordDM);
579:     DMSetCoordinateDM(dmAux, coordDM);
580:     DMGetDS(dmAux, &probAux);
581:     PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
582:     PetscObjectCompose((PetscObject) cdm, "dmAux", (PetscObject) dmAux);
583:     SetupMaterial(cdm, dmAux, user);
584:     DMDestroy(&dmAux);

586:     const PetscInt Ncomp = dim;
587:     const PetscInt components[] = {0,1,2};
588:     const PetscInt Nfid = 1;
589:     const PetscInt fid[] = {1}; /* The fixed faces */
590:     const PetscInt Npid = 1;
591:     const PetscInt pid[] = {2}; /* The faces with pressure loading */
592:     DMAddBoundary(cdm, PETSC_TRUE, "fixed", "Faces", 0, Ncomp, components, (void (*)()) coordinates, Nfid, fid, user);
593:     DMAddBoundary(cdm, PETSC_FALSE, "pressure", "Faces", 0, Ncomp, components, NULL, Npid, pid, user);
594:     DMGetCoarseDM(cdm, &cdm);
595:   }

597:   {
598:     /* Set up the near null space (a.k.a. rigid body modes) that will be used by the multigrid preconditioner */
599:     DM           subdm;
600:     MatNullSpace nearNullSpace;
601:     PetscInt     fields = 0;
602:     PetscObject  deformation;
603:     DMCreateSubDM(dm, 1, &fields, NULL, &subdm);
604:     DMPlexCreateRigidBody(subdm, &nearNullSpace);
605:     DMGetField(dm, 0, &deformation);
606:     PetscObjectCompose(deformation, "nearnullspace", (PetscObject) nearNullSpace);
607:     DMDestroy(&subdm);
608:     MatNullSpaceDestroy(&nearNullSpace);
609:   }

611:   PetscFEDestroy(&fe[0]);
612:   PetscFEDestroy(&fe[1]);
613:   PetscFEDestroy(&feBd[0]);
614:   PetscFEDestroy(&feBd[1]);
615:   PetscFEDestroy(&feAux);
616:   return(0);
617: }


622: int main(int argc, char **argv)
623: {
624:   SNES           snes;                 /* nonlinear solver */
625:   DM             dm;                   /* problem definition */
626:   Vec            u,r;                  /* solution, residual vectors */
627:   Mat            A,J;                  /* Jacobian matrix */
628:   MatNullSpace   nullSpace;            /* May be necessary for pressure */
629:   AppCtx         user;                 /* user-defined work context */
630:   PetscInt       its;                  /* iterations for convergence */

633:   PetscInitialize(&argc, &argv, NULL, help);
634:   ProcessOptions(PETSC_COMM_WORLD, &user);
635:   SNESCreate(PETSC_COMM_WORLD, &snes);
636:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
637:   SNESSetDM(snes, dm);
638:   DMSetApplicationContext(dm, &user);

640:   SetupDiscretization(dm, &user);
641:   DMPlexCreateClosureIndex(dm, NULL);

643:   DMCreateGlobalVector(dm, &u);
644:   VecDuplicate(u, &r);

646:   DMSetMatType(dm,MATAIJ);
647:   DMCreateMatrix(dm, &J);
648:   A = J;

650:   DMPlexSetSNESLocalFEM(dm,&user,&user,&user);
651:   SNESSetJacobian(snes, A, J, NULL, NULL);

653:   SNESSetFromOptions(snes);

655:   PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void* ctx) = {coordinates, zero_scalar};
656:   DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
657:   if (user.showInitial) {DMVecViewLocal(dm, u, PETSC_VIEWER_STDOUT_SELF);}

659:   if (user.runType == RUN_FULL) {
660:     if (user.debug) {
661:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
662:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
663:     }
664:     SNESSolve(snes, NULL, u);
665:     SNESGetIterationNumber(snes, &its);
666:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
667:     if (user.showSolution) {
668:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
669:       VecChop(u, 3.0e-9);
670:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
671:     }
672:   } else {
673:     PetscReal res = 0.0;

675:     /* Check initial guess */
676:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
677:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
678:     /* Check residual */
679:     SNESComputeFunction(snes, u, r);
680:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
681:     VecChop(r, 1.0e-10);
682:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
683:     VecNorm(r, NORM_2, &res);
684:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
685:     /* Check Jacobian */
686:     {
687:       Vec          b;
688:       PetscBool    isNull;

690:       SNESComputeJacobian(snes, u, A, A);
691:       VecDuplicate(u, &b);
692:       VecSet(r, 0.0);
693:       SNESComputeFunction(snes, r, b);
694:       MatMult(A, u, r);
695:       VecAXPY(r, 1.0, b);
696:       VecDestroy(&b);
697:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
698:       VecChop(r, 1.0e-10);
699:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
700:       VecNorm(r, NORM_2, &res);
701:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
702:     }
703:   }
704:   VecViewFromOptions(u, NULL, "-sol_vec_view");

706:   if (A != J) {MatDestroy(&A);}
707:   MatDestroy(&J);
708:   VecDestroy(&u);
709:   VecDestroy(&r);
710:   SNESDestroy(&snes);
711:   DMDestroy(&dm);
712:   PetscFinalize();
713:   return 0;
714: }