Actual source code: ex16.c


  2: static char help[] = "Large-deformation Elasticity Buckling Example";

  4: /*F-----------------------------------------------------------------------

  6:     This example solves the 3D large deformation elasticity problem

  8: \begin{equation}
  9:  \int_{\Omega}F \cdot S : \nabla v d\Omega + \int_{\Omega} (loading)\mathbf{e}_y \cdot v d\Omega = 0
 10: \end{equation}

 12:     F is the deformation gradient, and S is the second Piola-Kirchhoff tensor from the Saint Venant-Kirchhoff model of
 13:     hyperelasticity.  \Omega is a (arc) angle subsection of a cylindrical shell of thickness (height), inner radius
 14:     (rad) and width (width).  The problem is discretized using Q1 finite elements on a logically structured grid.
 15:     Homogenous Dirichlet boundary conditions are applied at the centers of the ends of the sphere.

 17:     This example is tunable with the following options:
 18:     -rad : the radius of the circle
 19:     -arc : set the angle of the arch represented
 20:     -loading : set the bulk loading (the mass)
 21:     -ploading : set the point loading at the top
 22:     -height : set the height of the arch
 23:     -width : set the width of the arch
 24:     -view_line : print initial and final offsets of the centerline of the
 25:                  beam along the x direction

 27:     The material properties may be modified using either:
 28:     -mu      : the first lame' parameter
 29:     -lambda  : the second lame' parameter

 31:     Or:
 32:     -young   : the Young's modulus
 33:     -poisson : the poisson ratio

 35:     This example is meant to show the strain placed upon the nonlinear solvers when trying to "snap through" the arch
 36:     using the loading.  Under certain parameter regimes, the arch will invert under the load, and the number of Newton
 37:     steps will jump considerably.  Composed nonlinear solvers may be used to mitigate this difficulty.

 39:     The initial setup follows the example in pg. 268 of "Nonlinear Finite Element Methods" by Peter Wriggers, but is a
 40:     3D extension.

 42:   ------------------------------------------------------------------------F*/

 44: #include <petscsnes.h>
 45: #include <petscdm.h>
 46: #include <petscdmda.h>

 48: #define QP0 0.2113248654051871
 49: #define QP1 0.7886751345948129
 50: #define NQ  2
 51: #define NB  2
 52: #define NEB 8
 53: #define NEQ 8
 54: #define NPB 24

 56: #define NVALS NEB *NEQ
 57: const PetscReal pts[NQ] = {QP0, QP1};
 58: const PetscReal wts[NQ] = {0.5, 0.5};

 60: PetscScalar vals[NVALS];
 61: PetscScalar grad[3 * NVALS];

 63: typedef PetscScalar Field[3];
 64: typedef PetscScalar CoordField[3];

 66: typedef PetscScalar JacField[9];

 68: PetscErrorCode FormFunctionLocal(DMDALocalInfo *, Field ***, Field ***, void *);
 69: PetscErrorCode FormJacobianLocal(DMDALocalInfo *, Field ***, Mat, Mat, void *);
 70: PetscErrorCode DisplayLine(SNES, Vec);
 71: PetscErrorCode FormElements(void);

 73: typedef struct {
 74:   PetscReal loading;
 75:   PetscReal mu;
 76:   PetscReal lambda;
 77:   PetscReal rad;
 78:   PetscReal height;
 79:   PetscReal width;
 80:   PetscReal arc;
 81:   PetscReal ploading;
 82: } AppCtx;

 84: PetscErrorCode        InitialGuess(DM, AppCtx *, Vec);
 85: PetscErrorCode        FormRHS(DM, AppCtx *, Vec);
 86: PetscErrorCode        FormCoordinates(DM, AppCtx *);
 87: extern PetscErrorCode NonlinearGS(SNES, Vec, Vec, void *);

 89: int main(int argc, char **argv)
 90: {
 91:   AppCtx    user; /* user-defined work context */
 92:   PetscInt  mx, my, its;
 93:   MPI_Comm  comm;
 94:   SNES      snes;
 95:   DM        da;
 96:   Vec       x, X, b;
 97:   PetscBool youngflg, poissonflg, muflg, lambdaflg, view = PETSC_FALSE, viewline = PETSC_FALSE;
 98:   PetscReal poisson = 0.2, young = 4e4;
 99:   char      filename[PETSC_MAX_PATH_LEN]     = "ex16.vts";
100:   char      filename_def[PETSC_MAX_PATH_LEN] = "ex16_def.vts";

103:   PetscInitialize(&argc, &argv, (char *)0, help);
104:   FormElements();
105:   comm = PETSC_COMM_WORLD;
106:   SNESCreate(comm, &snes);
107:   DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 11, 2, 2, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, NULL, NULL, &da);
108:   DMSetFromOptions(da);
109:   DMSetUp(da);
110:   SNESSetDM(snes, (DM)da);

112:   SNESSetNGS(snes, NonlinearGS, &user);

114:   DMDAGetInfo(da, 0, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE);
115:   user.loading  = 0.0;
116:   user.arc      = PETSC_PI / 3.;
117:   user.mu       = 4.0;
118:   user.lambda   = 1.0;
119:   user.rad      = 100.0;
120:   user.height   = 3.;
121:   user.width    = 1.;
122:   user.ploading = -5e3;

124:   PetscOptionsGetReal(NULL, NULL, "-arc", &user.arc, NULL);
125:   PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, &muflg);
126:   PetscOptionsGetReal(NULL, NULL, "-lambda", &user.lambda, &lambdaflg);
127:   PetscOptionsGetReal(NULL, NULL, "-rad", &user.rad, NULL);
128:   PetscOptionsGetReal(NULL, NULL, "-height", &user.height, NULL);
129:   PetscOptionsGetReal(NULL, NULL, "-width", &user.width, NULL);
130:   PetscOptionsGetReal(NULL, NULL, "-loading", &user.loading, NULL);
131:   PetscOptionsGetReal(NULL, NULL, "-ploading", &user.ploading, NULL);
132:   PetscOptionsGetReal(NULL, NULL, "-poisson", &poisson, &poissonflg);
133:   PetscOptionsGetReal(NULL, NULL, "-young", &young, &youngflg);
134:   if ((youngflg || poissonflg) || !(muflg || lambdaflg)) {
135:     /* set the lame' parameters based upon the poisson ratio and young's modulus */
136:     user.lambda = poisson * young / ((1. + poisson) * (1. - 2. * poisson));
137:     user.mu     = young / (2. * (1. + poisson));
138:   }
139:   PetscOptionsGetBool(NULL, NULL, "-view", &view, NULL);
140:   PetscOptionsGetBool(NULL, NULL, "-view_line", &viewline, NULL);

142:   DMDASetFieldName(da, 0, "x_disp");
143:   DMDASetFieldName(da, 1, "y_disp");
144:   DMDASetFieldName(da, 2, "z_disp");

146:   DMSetApplicationContext(da, &user);
147:   DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode(*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user);
148:   DMDASNESSetJacobianLocal(da, (DMDASNESJacobian)FormJacobianLocal, &user);
149:   SNESSetFromOptions(snes);
150:   FormCoordinates(da, &user);

152:   DMCreateGlobalVector(da, &x);
153:   DMCreateGlobalVector(da, &b);
154:   InitialGuess(da, &user, x);
155:   FormRHS(da, &user, b);

157:   PetscPrintf(comm, "lambda: %f mu: %f\n", (double)user.lambda, (double)user.mu);

159:   /* show a cross-section of the initial state */
160:   if (viewline) DisplayLine(snes, x);

162:   /* get the loaded configuration */
163:   SNESSolve(snes, b, x);

165:   SNESGetIterationNumber(snes, &its);
166:   PetscPrintf(comm, "Number of SNES iterations = %" PetscInt_FMT "\n", its);
167:   SNESGetSolution(snes, &X);
168:   /* show a cross-section of the final state */
169:   if (viewline) DisplayLine(snes, X);

171:   if (view) {
172:     PetscViewer viewer;
173:     Vec         coords;
174:     PetscViewerVTKOpen(PETSC_COMM_WORLD, filename, FILE_MODE_WRITE, &viewer);
175:     VecView(x, viewer);
176:     PetscViewerDestroy(&viewer);
177:     DMGetCoordinates(da, &coords);
178:     VecAXPY(coords, 1.0, x);
179:     PetscViewerVTKOpen(PETSC_COMM_WORLD, filename_def, FILE_MODE_WRITE, &viewer);
180:     VecView(x, viewer);
181:     PetscViewerDestroy(&viewer);
182:   }

184:   VecDestroy(&x);
185:   VecDestroy(&b);
186:   DMDestroy(&da);
187:   SNESDestroy(&snes);
188:   PetscFinalize();
189:   return 0;
190: }

192: PetscInt OnBoundary(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz)
193: {
194:   if ((i == 0 || i == mx - 1) && j == my - 1) return 1;
195:   return 0;
196: }

198: void BoundaryValue(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz, PetscScalar *val, AppCtx *user)
199: {
200:   /* reference coordinates */
201:   PetscReal p_x = ((PetscReal)i) / (((PetscReal)(mx - 1)));
202:   PetscReal p_y = ((PetscReal)j) / (((PetscReal)(my - 1)));
203:   PetscReal p_z = ((PetscReal)k) / (((PetscReal)(mz - 1)));
204:   PetscReal o_x = p_x;
205:   PetscReal o_y = p_y;
206:   PetscReal o_z = p_z;
207:   val[0]        = o_x - p_x;
208:   val[1]        = o_y - p_y;
209:   val[2]        = o_z - p_z;
210: }

212: void InvertTensor(PetscScalar *t, PetscScalar *ti, PetscReal *dett)
213: {
214:   const PetscScalar a   = t[0];
215:   const PetscScalar b   = t[1];
216:   const PetscScalar c   = t[2];
217:   const PetscScalar d   = t[3];
218:   const PetscScalar e   = t[4];
219:   const PetscScalar f   = t[5];
220:   const PetscScalar g   = t[6];
221:   const PetscScalar h   = t[7];
222:   const PetscScalar i   = t[8];
223:   const PetscReal   det = PetscRealPart(a * (e * i - f * h) - b * (i * d - f * g) + c * (d * h - e * g));
224:   const PetscReal   di  = 1. / det;
225:   if (ti) {
226:     const PetscScalar A  = (e * i - f * h);
227:     const PetscScalar B  = -(d * i - f * g);
228:     const PetscScalar C  = (d * h - e * g);
229:     const PetscScalar D  = -(b * i - c * h);
230:     const PetscScalar E  = (a * i - c * g);
231:     const PetscScalar F  = -(a * h - b * g);
232:     const PetscScalar G  = (b * f - c * e);
233:     const PetscScalar H  = -(a * f - c * d);
234:     const PetscScalar II = (a * e - b * d);
235:     ti[0]                = di * A;
236:     ti[1]                = di * D;
237:     ti[2]                = di * G;
238:     ti[3]                = di * B;
239:     ti[4]                = di * E;
240:     ti[5]                = di * H;
241:     ti[6]                = di * C;
242:     ti[7]                = di * F;
243:     ti[8]                = di * II;
244:   }
245:   if (dett) *dett = det;
246: }

248: void TensorTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
249: {
250:   PetscInt i, j, m;
251:   for (i = 0; i < 3; i++) {
252:     for (j = 0; j < 3; j++) {
253:       c[i + 3 * j] = 0;
254:       for (m = 0; m < 3; m++) c[i + 3 * j] += a[m + 3 * j] * b[i + 3 * m];
255:     }
256:   }
257: }

259: void TensorTransposeTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
260: {
261:   PetscInt i, j, m;
262:   for (i = 0; i < 3; i++) {
263:     for (j = 0; j < 3; j++) {
264:       c[i + 3 * j] = 0;
265:       for (m = 0; m < 3; m++) c[i + 3 * j] += a[3 * m + j] * b[i + 3 * m];
266:     }
267:   }
268: }

270: void TensorVector(PetscScalar *rot, PetscScalar *vec, PetscScalar *tvec)
271: {
272:   tvec[0] = rot[0] * vec[0] + rot[1] * vec[1] + rot[2] * vec[2];
273:   tvec[1] = rot[3] * vec[0] + rot[4] * vec[1] + rot[5] * vec[2];
274:   tvec[2] = rot[6] * vec[0] + rot[7] * vec[1] + rot[8] * vec[2];
275: }

277: void DeformationGradient(Field *ex, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *invJ, PetscScalar *F)
278: {
279:   PetscInt ii, jj, kk, l;
280:   for (l = 0; l < 9; l++) F[l] = 0.;
281:   F[0] = 1.;
282:   F[4] = 1.;
283:   F[8] = 1.;
284:   /* form the deformation gradient at this basis function -- loop over element unknowns */
285:   for (kk = 0; kk < NB; kk++) {
286:     for (jj = 0; jj < NB; jj++) {
287:       for (ii = 0; ii < NB; ii++) {
288:         PetscInt    idx  = ii + jj * NB + kk * NB * NB;
289:         PetscInt    bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
290:         PetscScalar lgrad[3];
291:         TensorVector(invJ, &grad[3 * bidx], lgrad);
292:         F[0] += lgrad[0] * ex[idx][0];
293:         F[1] += lgrad[1] * ex[idx][0];
294:         F[2] += lgrad[2] * ex[idx][0];
295:         F[3] += lgrad[0] * ex[idx][1];
296:         F[4] += lgrad[1] * ex[idx][1];
297:         F[5] += lgrad[2] * ex[idx][1];
298:         F[6] += lgrad[0] * ex[idx][2];
299:         F[7] += lgrad[1] * ex[idx][2];
300:         F[8] += lgrad[2] * ex[idx][2];
301:       }
302:     }
303:   }
304: }

306: void DeformationGradientJacobian(PetscInt qi, PetscInt qj, PetscInt qk, PetscInt ii, PetscInt jj, PetscInt kk, PetscInt fld, PetscScalar *invJ, PetscScalar *dF)
307: {
308:   PetscInt    l;
309:   PetscScalar lgrad[3];
310:   PetscInt    idx  = ii + jj * NB + kk * NB * NB;
311:   PetscInt    bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
312:   for (l = 0; l < 9; l++) dF[l] = 0.;
313:   /* form the deformation gradient at this basis function -- loop over element unknowns */
314:   TensorVector(invJ, &grad[3 * bidx], lgrad);
315:   dF[3 * fld]     = lgrad[0];
316:   dF[3 * fld + 1] = lgrad[1];
317:   dF[3 * fld + 2] = lgrad[2];
318: }

320: void LagrangeGreenStrain(PetscScalar *F, PetscScalar *E)
321: {
322:   PetscInt i, j, m;
323:   for (i = 0; i < 3; i++) {
324:     for (j = 0; j < 3; j++) {
325:       E[i + 3 * j] = 0;
326:       for (m = 0; m < 3; m++) E[i + 3 * j] += 0.5 * F[3 * m + j] * F[i + 3 * m];
327:     }
328:   }
329:   for (i = 0; i < 3; i++) E[i + 3 * i] -= 0.5;
330: }

332: void SaintVenantKirchoff(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *S)
333: {
334:   PetscInt    i, j;
335:   PetscScalar E[9];
336:   PetscScalar trE = 0;
337:   LagrangeGreenStrain(F, E);
338:   for (i = 0; i < 3; i++) trE += E[i + 3 * i];
339:   for (i = 0; i < 3; i++) {
340:     for (j = 0; j < 3; j++) {
341:       S[i + 3 * j] = 2. * mu * E[i + 3 * j];
342:       if (i == j) S[i + 3 * j] += trE * lambda;
343:     }
344:   }
345: }

347: void SaintVenantKirchoffJacobian(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *dF, PetscScalar *dS)
348: {
349:   PetscScalar FtdF[9], dE[9];
350:   PetscInt    i, j;
351:   PetscScalar dtrE = 0.;
352:   TensorTransposeTensor(dF, F, dE);
353:   TensorTransposeTensor(F, dF, FtdF);
354:   for (i = 0; i < 9; i++) dE[i] += FtdF[i];
355:   for (i = 0; i < 9; i++) dE[i] *= 0.5;
356:   for (i = 0; i < 3; i++) dtrE += dE[i + 3 * i];
357:   for (i = 0; i < 3; i++) {
358:     for (j = 0; j < 3; j++) {
359:       dS[i + 3 * j] = 2. * mu * dE[i + 3 * j];
360:       if (i == j) dS[i + 3 * j] += dtrE * lambda;
361:     }
362:   }
363: }

365: PetscErrorCode FormElements()
366: {
367:   PetscInt  i, j, k, ii, jj, kk;
368:   PetscReal bx, by, bz, dbx, dby, dbz;

370:   /* construct the basis function values and derivatives */
371:   for (k = 0; k < NB; k++) {
372:     for (j = 0; j < NB; j++) {
373:       for (i = 0; i < NB; i++) {
374:         /* loop over the quadrature points */
375:         for (kk = 0; kk < NQ; kk++) {
376:           for (jj = 0; jj < NQ; jj++) {
377:             for (ii = 0; ii < NQ; ii++) {
378:               PetscInt idx = ii + NQ * jj + NQ * NQ * kk + NEQ * i + NEQ * NB * j + NEQ * NB * NB * k;
379:               bx           = pts[ii];
380:               by           = pts[jj];
381:               bz           = pts[kk];
382:               dbx          = 1.;
383:               dby          = 1.;
384:               dbz          = 1.;
385:               if (i == 0) {
386:                 bx  = 1. - bx;
387:                 dbx = -1;
388:               }
389:               if (j == 0) {
390:                 by  = 1. - by;
391:                 dby = -1;
392:               }
393:               if (k == 0) {
394:                 bz  = 1. - bz;
395:                 dbz = -1;
396:               }
397:               vals[idx]         = bx * by * bz;
398:               grad[3 * idx + 0] = dbx * by * bz;
399:               grad[3 * idx + 1] = dby * bx * bz;
400:               grad[3 * idx + 2] = dbz * bx * by;
401:             }
402:           }
403:         }
404:       }
405:     }
406:   }
407:   return 0;
408: }

410: void GatherElementData(PetscInt mx, PetscInt my, PetscInt mz, Field ***x, CoordField ***c, PetscInt i, PetscInt j, PetscInt k, Field *ex, CoordField *ec, AppCtx *user)
411: {
412:   PetscInt m;
413:   PetscInt ii, jj, kk;
414:   /* gather the data -- loop over element unknowns */
415:   for (kk = 0; kk < NB; kk++) {
416:     for (jj = 0; jj < NB; jj++) {
417:       for (ii = 0; ii < NB; ii++) {
418:         PetscInt idx = ii + jj * NB + kk * NB * NB;
419:         /* decouple the boundary nodes for the displacement variables */
420:         if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
421:           BoundaryValue(i + ii, j + jj, k + kk, mx, my, mz, ex[idx], user);
422:         } else {
423:           for (m = 0; m < 3; m++) ex[idx][m] = x[k + kk][j + jj][i + ii][m];
424:         }
425:         for (m = 0; m < 3; m++) ec[idx][m] = c[k + kk][j + jj][i + ii][m];
426:       }
427:     }
428:   }
429: }

431: void QuadraturePointGeometricJacobian(CoordField *ec, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *J)
432: {
433:   PetscInt ii, jj, kk;
434:   /* construct the gradient at the given quadrature point named by i,j,k */
435:   for (ii = 0; ii < 9; ii++) J[ii] = 0;
436:   for (kk = 0; kk < NB; kk++) {
437:     for (jj = 0; jj < NB; jj++) {
438:       for (ii = 0; ii < NB; ii++) {
439:         PetscInt idx  = ii + jj * NB + kk * NB * NB;
440:         PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
441:         J[0] += grad[3 * bidx + 0] * ec[idx][0];
442:         J[1] += grad[3 * bidx + 1] * ec[idx][0];
443:         J[2] += grad[3 * bidx + 2] * ec[idx][0];
444:         J[3] += grad[3 * bidx + 0] * ec[idx][1];
445:         J[4] += grad[3 * bidx + 1] * ec[idx][1];
446:         J[5] += grad[3 * bidx + 2] * ec[idx][1];
447:         J[6] += grad[3 * bidx + 0] * ec[idx][2];
448:         J[7] += grad[3 * bidx + 1] * ec[idx][2];
449:         J[8] += grad[3 * bidx + 2] * ec[idx][2];
450:       }
451:     }
452:   }
453: }

455: void FormElementJacobian(Field *ex, CoordField *ec, Field *ef, PetscScalar *ej, AppCtx *user)
456: {
457:   PetscReal   vol;
458:   PetscScalar J[9];
459:   PetscScalar invJ[9];
460:   PetscScalar F[9], S[9], dF[9], dS[9], dFS[9], FdS[9], FS[9];
461:   PetscReal   scl;
462:   PetscInt    i, j, k, l, ii, jj, kk, ll, qi, qj, qk, m;

464:   if (ej)
465:     for (i = 0; i < NPB * NPB; i++) ej[i] = 0.;
466:   if (ef)
467:     for (i = 0; i < NEB; i++) {
468:       ef[i][0] = 0.;
469:       ef[i][1] = 0.;
470:       ef[i][2] = 0.;
471:     }
472:   /* loop over quadrature */
473:   for (qk = 0; qk < NQ; qk++) {
474:     for (qj = 0; qj < NQ; qj++) {
475:       for (qi = 0; qi < NQ; qi++) {
476:         QuadraturePointGeometricJacobian(ec, qi, qj, qk, J);
477:         InvertTensor(J, invJ, &vol);
478:         scl = vol * wts[qi] * wts[qj] * wts[qk];
479:         DeformationGradient(ex, qi, qj, qk, invJ, F);
480:         SaintVenantKirchoff(user->lambda, user->mu, F, S);
481:         /* form the function */
482:         if (ef) {
483:           TensorTensor(F, S, FS);
484:           for (kk = 0; kk < NB; kk++) {
485:             for (jj = 0; jj < NB; jj++) {
486:               for (ii = 0; ii < NB; ii++) {
487:                 PetscInt    idx  = ii + jj * NB + kk * NB * NB;
488:                 PetscInt    bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
489:                 PetscScalar lgrad[3];
490:                 TensorVector(invJ, &grad[3 * bidx], lgrad);
491:                 /* mu*F : grad phi_{u,v,w} */
492:                 for (m = 0; m < 3; m++) ef[idx][m] += scl * (lgrad[0] * FS[3 * m + 0] + lgrad[1] * FS[3 * m + 1] + lgrad[2] * FS[3 * m + 2]);
493:                 ef[idx][1] -= scl * user->loading * vals[bidx];
494:               }
495:             }
496:           }
497:         }
498:         /* form the jacobian */
499:         if (ej) {
500:           /* loop over trialfunctions */
501:           for (k = 0; k < NB; k++) {
502:             for (j = 0; j < NB; j++) {
503:               for (i = 0; i < NB; i++) {
504:                 for (l = 0; l < 3; l++) {
505:                   PetscInt tridx = l + 3 * (i + j * NB + k * NB * NB);
506:                   DeformationGradientJacobian(qi, qj, qk, i, j, k, l, invJ, dF);
507:                   SaintVenantKirchoffJacobian(user->lambda, user->mu, F, dF, dS);
508:                   TensorTensor(dF, S, dFS);
509:                   TensorTensor(F, dS, FdS);
510:                   for (m = 0; m < 9; m++) dFS[m] += FdS[m];
511:                   /* loop over testfunctions */
512:                   for (kk = 0; kk < NB; kk++) {
513:                     for (jj = 0; jj < NB; jj++) {
514:                       for (ii = 0; ii < NB; ii++) {
515:                         PetscInt    idx  = ii + jj * NB + kk * NB * NB;
516:                         PetscInt    bidx = 8 * idx + qi + NQ * qj + NQ * NQ * qk;
517:                         PetscScalar lgrad[3];
518:                         TensorVector(invJ, &grad[3 * bidx], lgrad);
519:                         for (ll = 0; ll < 3; ll++) {
520:                           PetscInt teidx = ll + 3 * (ii + jj * NB + kk * NB * NB);
521:                           ej[teidx + NPB * tridx] += scl * (lgrad[0] * dFS[3 * ll + 0] + lgrad[1] * dFS[3 * ll + 1] + lgrad[2] * dFS[3 * ll + 2]);
522:                         }
523:                       }
524:                     }
525:                   } /* end of testfunctions */
526:                 }
527:               }
528:             }
529:           } /* end of trialfunctions */
530:         }
531:       }
532:     }
533:   } /* end of quadrature points */
534: }

536: void FormPBJacobian(PetscInt i, PetscInt j, PetscInt k, Field *ex, CoordField *ec, Field *ef, PetscScalar *ej, AppCtx *user)
537: {
538:   PetscReal   vol;
539:   PetscScalar J[9];
540:   PetscScalar invJ[9];
541:   PetscScalar F[9], S[9], dF[9], dS[9], dFS[9], FdS[9], FS[9];
542:   PetscReal   scl;
543:   PetscInt    l, ll, qi, qj, qk, m;
544:   PetscInt    idx = i + j * NB + k * NB * NB;
545:   PetscScalar lgrad[3];

547:   if (ej)
548:     for (l = 0; l < 9; l++) ej[l] = 0.;
549:   if (ef)
550:     for (l = 0; l < 1; l++) {
551:       ef[l][0] = 0.;
552:       ef[l][1] = 0.;
553:       ef[l][2] = 0.;
554:     }
555:   /* loop over quadrature */
556:   for (qk = 0; qk < NQ; qk++) {
557:     for (qj = 0; qj < NQ; qj++) {
558:       for (qi = 0; qi < NQ; qi++) {
559:         PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
560:         QuadraturePointGeometricJacobian(ec, qi, qj, qk, J);
561:         InvertTensor(J, invJ, &vol);
562:         TensorVector(invJ, &grad[3 * bidx], lgrad);
563:         scl = vol * wts[qi] * wts[qj] * wts[qk];
564:         DeformationGradient(ex, qi, qj, qk, invJ, F);
565:         SaintVenantKirchoff(user->lambda, user->mu, F, S);
566:         /* form the function */
567:         if (ef) {
568:           TensorTensor(F, S, FS);
569:           for (m = 0; m < 3; m++) ef[0][m] += scl * (lgrad[0] * FS[3 * m + 0] + lgrad[1] * FS[3 * m + 1] + lgrad[2] * FS[3 * m + 2]);
570:           ef[0][1] -= scl * user->loading * vals[bidx];
571:         }
572:         /* form the jacobian */
573:         if (ej) {
574:           for (l = 0; l < 3; l++) {
575:             DeformationGradientJacobian(qi, qj, qk, i, j, k, l, invJ, dF);
576:             SaintVenantKirchoffJacobian(user->lambda, user->mu, F, dF, dS);
577:             TensorTensor(dF, S, dFS);
578:             TensorTensor(F, dS, FdS);
579:             for (m = 0; m < 9; m++) dFS[m] += FdS[m];
580:             for (ll = 0; ll < 3; ll++) ej[ll + 3 * l] += scl * (lgrad[0] * dFS[3 * ll + 0] + lgrad[1] * dFS[3 * ll + 1] + lgrad[2] * dFS[3 * ll + 2]);
581:           }
582:         }
583:       }
584:     } /* end of quadrature points */
585:   }
586: }

588: void ApplyBCsElement(PetscInt mx, PetscInt my, PetscInt mz, PetscInt i, PetscInt j, PetscInt k, PetscScalar *jacobian)
589: {
590:   PetscInt ii, jj, kk, ll, ei, ej, ek, el;
591:   for (kk = 0; kk < NB; kk++) {
592:     for (jj = 0; jj < NB; jj++) {
593:       for (ii = 0; ii < NB; ii++) {
594:         for (ll = 0; ll < 3; ll++) {
595:           PetscInt tridx = ll + 3 * (ii + jj * NB + kk * NB * NB);
596:           for (ek = 0; ek < NB; ek++) {
597:             for (ej = 0; ej < NB; ej++) {
598:               for (ei = 0; ei < NB; ei++) {
599:                 for (el = 0; el < 3; el++) {
600:                   if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz) || OnBoundary(i + ei, j + ej, k + ek, mx, my, mz)) {
601:                     PetscInt teidx = el + 3 * (ei + ej * NB + ek * NB * NB);
602:                     if (teidx == tridx) {
603:                       jacobian[tridx + NPB * teidx] = 1.;
604:                     } else {
605:                       jacobian[tridx + NPB * teidx] = 0.;
606:                     }
607:                   }
608:                 }
609:               }
610:             }
611:           }
612:         }
613:       }
614:     }
615:   }
616: }

618: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field ***x, Mat jacpre, Mat jac, void *ptr)
619: {
620:   /* values for each basis function at each quadrature point */
621:   AppCtx     *user = (AppCtx *)ptr;
622:   PetscInt    i, j, k, m, l;
623:   PetscInt    ii, jj, kk;
624:   PetscScalar ej[NPB * NPB];
625:   PetscScalar vals[NPB * NPB];
626:   Field       ex[NEB];
627:   CoordField  ec[NEB];

629:   PetscInt      xs = info->xs, ys = info->ys, zs = info->zs;
630:   PetscInt      xm = info->xm, ym = info->ym, zm = info->zm;
631:   PetscInt      xes, yes, zes, xee, yee, zee;
632:   PetscInt      mx = info->mx, my = info->my, mz = info->mz;
633:   DM            cda;
634:   CoordField ***c;
635:   Vec           C;
636:   PetscInt      nrows;
637:   MatStencil    col[NPB], row[NPB];
638:   PetscScalar   v[9];

640:   DMGetCoordinateDM(info->da, &cda);
641:   DMGetCoordinatesLocal(info->da, &C);
642:   DMDAVecGetArray(cda, C, &c);
643:   MatScale(jac, 0.0);

645:   xes = xs;
646:   yes = ys;
647:   zes = zs;
648:   xee = xs + xm;
649:   yee = ys + ym;
650:   zee = zs + zm;
651:   if (xs > 0) xes = xs - 1;
652:   if (ys > 0) yes = ys - 1;
653:   if (zs > 0) zes = zs - 1;
654:   if (xs + xm == mx) xee = xs + xm - 1;
655:   if (ys + ym == my) yee = ys + ym - 1;
656:   if (zs + zm == mz) zee = zs + zm - 1;
657:   for (k = zes; k < zee; k++) {
658:     for (j = yes; j < yee; j++) {
659:       for (i = xes; i < xee; i++) {
660:         GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
661:         FormElementJacobian(ex, ec, NULL, ej, user);
662:         ApplyBCsElement(mx, my, mz, i, j, k, ej);
663:         nrows = 0.;
664:         for (kk = 0; kk < NB; kk++) {
665:           for (jj = 0; jj < NB; jj++) {
666:             for (ii = 0; ii < NB; ii++) {
667:               PetscInt idx = ii + jj * 2 + kk * 4;
668:               for (m = 0; m < 3; m++) {
669:                 col[3 * idx + m].i = i + ii;
670:                 col[3 * idx + m].j = j + jj;
671:                 col[3 * idx + m].k = k + kk;
672:                 col[3 * idx + m].c = m;
673:                 if (i + ii >= xs && i + ii < xm + xs && j + jj >= ys && j + jj < ys + ym && k + kk >= zs && k + kk < zs + zm) {
674:                   row[nrows].i = i + ii;
675:                   row[nrows].j = j + jj;
676:                   row[nrows].k = k + kk;
677:                   row[nrows].c = m;
678:                   for (l = 0; l < NPB; l++) vals[NPB * nrows + l] = ej[NPB * (3 * idx + m) + l];
679:                   nrows++;
680:                 }
681:               }
682:             }
683:           }
684:         }
685:         MatSetValuesStencil(jac, nrows, row, NPB, col, vals, ADD_VALUES);
686:       }
687:     }
688:   }

690:   MatAssemblyBegin(jac, MAT_FLUSH_ASSEMBLY);
691:   MatAssemblyEnd(jac, MAT_FLUSH_ASSEMBLY);

693:   /* set the diagonal */
694:   v[0] = 1.;
695:   v[1] = 0.;
696:   v[2] = 0.;
697:   v[3] = 0.;
698:   v[4] = 1.;
699:   v[5] = 0.;
700:   v[6] = 0.;
701:   v[7] = 0.;
702:   v[8] = 1.;
703:   for (k = zs; k < zs + zm; k++) {
704:     for (j = ys; j < ys + ym; j++) {
705:       for (i = xs; i < xs + xm; i++) {
706:         if (OnBoundary(i, j, k, mx, my, mz)) {
707:           for (m = 0; m < 3; m++) {
708:             col[m].i = i;
709:             col[m].j = j;
710:             col[m].k = k;
711:             col[m].c = m;
712:           }
713:           MatSetValuesStencil(jac, 3, col, 3, col, v, INSERT_VALUES);
714:         }
715:       }
716:     }
717:   }

719:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
720:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);

722:   DMDAVecRestoreArray(cda, C, &c);
723:   return 0;
724: }

726: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field ***x, Field ***f, void *ptr)
727: {
728:   /* values for each basis function at each quadrature point */
729:   AppCtx  *user = (AppCtx *)ptr;
730:   PetscInt i, j, k, l;
731:   PetscInt ii, jj, kk;

733:   Field      ef[NEB];
734:   Field      ex[NEB];
735:   CoordField ec[NEB];

737:   PetscInt      xs = info->xs, ys = info->ys, zs = info->zs;
738:   PetscInt      xm = info->xm, ym = info->ym, zm = info->zm;
739:   PetscInt      xes, yes, zes, xee, yee, zee;
740:   PetscInt      mx = info->mx, my = info->my, mz = info->mz;
741:   DM            cda;
742:   CoordField ***c;
743:   Vec           C;

745:   DMGetCoordinateDM(info->da, &cda);
746:   DMGetCoordinatesLocal(info->da, &C);
747:   DMDAVecGetArray(cda, C, &c);
748:   DMDAGetInfo(info->da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
749:   DMDAGetCorners(info->da, &xs, &ys, &zs, &xm, &ym, &zm);

751:   /* loop over elements */
752:   for (k = zs; k < zs + zm; k++) {
753:     for (j = ys; j < ys + ym; j++) {
754:       for (i = xs; i < xs + xm; i++) {
755:         for (l = 0; l < 3; l++) f[k][j][i][l] = 0.;
756:       }
757:     }
758:   }
759:   /* element starts and ends */
760:   xes = xs;
761:   yes = ys;
762:   zes = zs;
763:   xee = xs + xm;
764:   yee = ys + ym;
765:   zee = zs + zm;
766:   if (xs > 0) xes = xs - 1;
767:   if (ys > 0) yes = ys - 1;
768:   if (zs > 0) zes = zs - 1;
769:   if (xs + xm == mx) xee = xs + xm - 1;
770:   if (ys + ym == my) yee = ys + ym - 1;
771:   if (zs + zm == mz) zee = zs + zm - 1;
772:   for (k = zes; k < zee; k++) {
773:     for (j = yes; j < yee; j++) {
774:       for (i = xes; i < xee; i++) {
775:         GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
776:         FormElementJacobian(ex, ec, ef, NULL, user);
777:         /* put this element's additions into the residuals */
778:         for (kk = 0; kk < NB; kk++) {
779:           for (jj = 0; jj < NB; jj++) {
780:             for (ii = 0; ii < NB; ii++) {
781:               PetscInt idx = ii + jj * NB + kk * NB * NB;
782:               if (k + kk >= zs && j + jj >= ys && i + ii >= xs && k + kk < zs + zm && j + jj < ys + ym && i + ii < xs + xm) {
783:                 if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
784:                   for (l = 0; l < 3; l++) f[k + kk][j + jj][i + ii][l] = x[k + kk][j + jj][i + ii][l] - ex[idx][l];
785:                 } else {
786:                   for (l = 0; l < 3; l++) f[k + kk][j + jj][i + ii][l] += ef[idx][l];
787:                 }
788:               }
789:             }
790:           }
791:         }
792:       }
793:     }
794:   }
795:   DMDAVecRestoreArray(cda, C, &c);
796:   return 0;
797: }

799: PetscErrorCode NonlinearGS(SNES snes, Vec X, Vec B, void *ptr)
800: {
801:   /* values for each basis function at each quadrature point */
802:   AppCtx       *user = (AppCtx *)ptr;
803:   PetscInt      i, j, k, l, m, n, s;
804:   PetscInt      pi, pj, pk;
805:   Field         ef[1];
806:   Field         ex[8];
807:   PetscScalar   ej[9];
808:   CoordField    ec[8];
809:   PetscScalar   pjac[9], pjinv[9];
810:   PetscScalar   pf[3], py[3];
811:   PetscInt      xs, ys, zs;
812:   PetscInt      xm, ym, zm;
813:   PetscInt      mx, my, mz;
814:   DM            cda;
815:   CoordField ***c;
816:   Vec           C;
817:   DM            da;
818:   Vec           Xl, Bl;
819:   Field      ***x, ***b;
820:   PetscInt      sweeps, its;
821:   PetscReal     atol, rtol, stol;
822:   PetscReal     fnorm0 = 0.0, fnorm, ynorm, xnorm = 0.0;

824:   SNESNGSGetSweeps(snes, &sweeps);
825:   SNESNGSGetTolerances(snes, &atol, &rtol, &stol, &its);

827:   SNESGetDM(snes, &da);
828:   DMGetLocalVector(da, &Xl);
829:   if (B) DMGetLocalVector(da, &Bl);
830:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xl);
831:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xl);
832:   if (B) {
833:     DMGlobalToLocalBegin(da, B, INSERT_VALUES, Bl);
834:     DMGlobalToLocalEnd(da, B, INSERT_VALUES, Bl);
835:   }
836:   DMDAVecGetArray(da, Xl, &x);
837:   if (B) DMDAVecGetArray(da, Bl, &b);

839:   DMGetCoordinateDM(da, &cda);
840:   DMGetCoordinatesLocal(da, &C);
841:   DMDAVecGetArray(cda, C, &c);
842:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
843:   DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);

845:   for (s = 0; s < sweeps; s++) {
846:     for (k = zs; k < zs + zm; k++) {
847:       for (j = ys; j < ys + ym; j++) {
848:         for (i = xs; i < xs + xm; i++) {
849:           if (OnBoundary(i, j, k, mx, my, mz)) {
850:             BoundaryValue(i, j, k, mx, my, mz, x[k][j][i], user);
851:           } else {
852:             for (n = 0; n < its; n++) {
853:               for (m = 0; m < 9; m++) pjac[m] = 0.;
854:               for (m = 0; m < 3; m++) pf[m] = 0.;
855:               /* gather the elements for this point */
856:               for (pk = -1; pk < 1; pk++) {
857:                 for (pj = -1; pj < 1; pj++) {
858:                   for (pi = -1; pi < 1; pi++) {
859:                     /* check that this element exists */
860:                     if (i + pi >= 0 && i + pi < mx - 1 && j + pj >= 0 && j + pj < my - 1 && k + pk >= 0 && k + pk < mz - 1) {
861:                       /* create the element function and jacobian */
862:                       GatherElementData(mx, my, mz, x, c, i + pi, j + pj, k + pk, ex, ec, user);
863:                       FormPBJacobian(-pi, -pj, -pk, ex, ec, ef, ej, user);
864:                       /* extract the point named by i,j,k from the whole element jacobian and function */
865:                       for (l = 0; l < 3; l++) {
866:                         pf[l] += ef[0][l];
867:                         for (m = 0; m < 3; m++) pjac[3 * m + l] += ej[3 * m + l];
868:                       }
869:                     }
870:                   }
871:                 }
872:               }
873:               /* invert */
874:               InvertTensor(pjac, pjinv, NULL);
875:               /* apply */
876:               if (B)
877:                 for (m = 0; m < 3; m++) pf[m] -= b[k][j][i][m];
878:               TensorVector(pjinv, pf, py);
879:               xnorm = 0.;
880:               for (m = 0; m < 3; m++) {
881:                 x[k][j][i][m] -= py[m];
882:                 xnorm += PetscRealPart(x[k][j][i][m] * x[k][j][i][m]);
883:               }
884:               fnorm = PetscRealPart(pf[0] * pf[0] + pf[1] * pf[1] + pf[2] * pf[2]);
885:               if (n == 0) fnorm0 = fnorm;
886:               ynorm = PetscRealPart(py[0] * py[0] + py[1] * py[1] + py[2] * py[2]);
887:               if (fnorm < atol * atol || fnorm < rtol * rtol * fnorm0 || ynorm < stol * stol * xnorm) break;
888:             }
889:           }
890:         }
891:       }
892:     }
893:   }
894:   DMDAVecRestoreArray(da, Xl, &x);
895:   DMLocalToGlobalBegin(da, Xl, INSERT_VALUES, X);
896:   DMLocalToGlobalEnd(da, Xl, INSERT_VALUES, X);
897:   DMRestoreLocalVector(da, &Xl);
898:   if (B) {
899:     DMDAVecRestoreArray(da, Bl, &b);
900:     DMRestoreLocalVector(da, &Bl);
901:   }
902:   DMDAVecRestoreArray(cda, C, &c);
903:   return 0;
904: }

906: PetscErrorCode FormCoordinates(DM da, AppCtx *user)
907: {
908:   Vec           coords;
909:   DM            cda;
910:   PetscInt      mx, my, mz;
911:   PetscInt      i, j, k, xs, ys, zs, xm, ym, zm;
912:   CoordField ***x;

914:   DMGetCoordinateDM(da, &cda);
915:   DMCreateGlobalVector(cda, &coords);
916:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
917:   DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
918:   DMDAVecGetArray(da, coords, &x);
919:   for (k = zs; k < zs + zm; k++) {
920:     for (j = ys; j < ys + ym; j++) {
921:       for (i = xs; i < xs + xm; i++) {
922:         PetscReal cx  = ((PetscReal)i) / (((PetscReal)(mx - 1)));
923:         PetscReal cy  = ((PetscReal)j) / (((PetscReal)(my - 1)));
924:         PetscReal cz  = ((PetscReal)k) / (((PetscReal)(mz - 1)));
925:         PetscReal rad = user->rad + cy * user->height;
926:         PetscReal ang = (cx - 0.5) * user->arc;
927:         x[k][j][i][0] = rad * PetscSinReal(ang);
928:         x[k][j][i][1] = rad * PetscCosReal(ang) - (user->rad + 0.5 * user->height) * PetscCosReal(-0.5 * user->arc);
929:         x[k][j][i][2] = user->width * (cz - 0.5);
930:       }
931:     }
932:   }
933:   DMDAVecRestoreArray(da, coords, &x);
934:   DMSetCoordinates(da, coords);
935:   VecDestroy(&coords);
936:   return 0;
937: }

939: PetscErrorCode InitialGuess(DM da, AppCtx *user, Vec X)
940: {
941:   PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
942:   PetscInt mx, my, mz;
943:   Field ***x;

945:   DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
946:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
947:   DMDAVecGetArray(da, X, &x);

949:   for (k = zs; k < zs + zm; k++) {
950:     for (j = ys; j < ys + ym; j++) {
951:       for (i = xs; i < xs + xm; i++) {
952:         /* reference coordinates */
953:         PetscReal p_x = ((PetscReal)i) / (((PetscReal)(mx - 1)));
954:         PetscReal p_y = ((PetscReal)j) / (((PetscReal)(my - 1)));
955:         PetscReal p_z = ((PetscReal)k) / (((PetscReal)(mz - 1)));
956:         PetscReal o_x = p_x;
957:         PetscReal o_y = p_y;
958:         PetscReal o_z = p_z;
959:         x[k][j][i][0] = o_x - p_x;
960:         x[k][j][i][1] = o_y - p_y;
961:         x[k][j][i][2] = o_z - p_z;
962:       }
963:     }
964:   }
965:   DMDAVecRestoreArray(da, X, &x);
966:   return 0;
967: }

969: PetscErrorCode FormRHS(DM da, AppCtx *user, Vec X)
970: {
971:   PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
972:   PetscInt mx, my, mz;
973:   Field ***x;

975:   DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
976:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
977:   DMDAVecGetArray(da, X, &x);

979:   for (k = zs; k < zs + zm; k++) {
980:     for (j = ys; j < ys + ym; j++) {
981:       for (i = xs; i < xs + xm; i++) {
982:         x[k][j][i][0] = 0.;
983:         x[k][j][i][1] = 0.;
984:         x[k][j][i][2] = 0.;
985:         if (i == (mx - 1) / 2 && j == (my - 1)) x[k][j][i][1] = user->ploading / (mz - 1);
986:       }
987:     }
988:   }
989:   DMDAVecRestoreArray(da, X, &x);
990:   return 0;
991: }

993: PetscErrorCode DisplayLine(SNES snes, Vec X)
994: {
995:   PetscInt      r, i, j = 0, k = 0, xs, xm, ys, ym, zs, zm, mx, my, mz;
996:   Field      ***x;
997:   CoordField ***c;
998:   DM            da, cda;
999:   Vec           C;
1000:   PetscMPIInt   size, rank;

1002:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
1003:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1004:   SNESGetDM(snes, &da);
1005:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
1006:   DMGetCoordinateDM(da, &cda);
1007:   DMGetCoordinates(da, &C);
1008:   j = my / 2;
1009:   k = mz / 2;
1010:   DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
1011:   DMDAVecGetArray(da, X, &x);
1012:   DMDAVecGetArray(cda, C, &c);
1013:   for (r = 0; r < size; r++) {
1014:     if (rank == r) {
1015:       if (j >= ys && j < ys + ym && k >= zs && k < zs + zm) {
1016:         for (i = xs; i < xs + xm; i++) {
1017:           PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %d %d: %f %f %f\n", i, 0, 0, (double)PetscRealPart(c[k][j][i][0] + x[k][j][i][0]), (double)PetscRealPart(c[k][j][i][1] + x[k][j][i][1]), (double)PetscRealPart(c[k][j][i][2] + x[k][j][i][2]));
1018:         }
1019:       }
1020:     }
1021:     MPI_Barrier(PETSC_COMM_WORLD);
1022:   }
1023:   DMDAVecRestoreArray(da, X, &x);
1024:   DMDAVecRestoreArray(cda, C, &c);
1025:   return 0;
1026: }

1028: /*TEST

1030:    test:
1031:       nsize: 2
1032:       args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading 0.0 -loading -1. -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -snes_max_it 7
1033:       requires: !single
1034:       timeoutfactor: 3

1036:    test:
1037:       suffix: 2
1038:       args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading 0.0 -loading -1. -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -npc_snes_type fas -npc_fas_levels_snes_type ncg -npc_fas_levels_snes_max_it 3 -npc_snes_monitor_short -snes_max_it 2
1039:       requires: !single

1041:    test:
1042:       suffix: 3
1043:       args: -da_refine 1 -da_overlap 3 -da_local_subdomains 4 -snes_type aspin -rad 10.0 -young 10. -ploading 0.0 -loading -0.5 -snes_monitor_short -ksp_monitor_short -npc_sub_snes_rtol 1e-2 -ksp_rtol 1e-2 -ksp_max_it 14 -snes_converged_reason -snes_max_linear_solve_fail 100 -snes_max_it 4
1044:       requires: !single

1046: TEST*/