Actual source code: ex16.c

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

  3: /*F-----------------------------------------------------------------------

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

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

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

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

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

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

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

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

 41:   ------------------------------------------------------------------------F*/

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

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

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

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

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

 65: typedef PetscScalar JacField[9];

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

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

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

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

101:   PetscFunctionBeginUser;
102:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
103:   PetscCall(FormElements());
104:   comm = PETSC_COMM_WORLD;
105:   PetscCall(SNESCreate(comm, &snes));
106:   PetscCall(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));
107:   PetscCall(DMSetFromOptions(da));
108:   PetscCall(DMSetUp(da));
109:   PetscCall(SNESSetDM(snes, (DM)da));

111:   PetscCall(SNESSetNGS(snes, NonlinearGS, &user));

113:   PetscCall(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));
114:   user.loading  = 0.0;
115:   user.arc      = PETSC_PI / 3.;
116:   user.mu       = 4.0;
117:   user.lambda   = 1.0;
118:   user.rad      = 100.0;
119:   user.height   = 3.;
120:   user.width    = 1.;
121:   user.ploading = -5e3;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

369:   PetscFunctionBegin;
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:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBegin;
641:   PetscCall(DMGetCoordinateDM(info->da, &cda));
642:   PetscCall(DMGetCoordinatesLocal(info->da, &C));
643:   PetscCall(DMDAVecGetArray(cda, C, &c));
644:   PetscCall(MatScale(jac, 0.0));

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

691:   PetscCall(MatAssemblyBegin(jac, MAT_FLUSH_ASSEMBLY));
692:   PetscCall(MatAssemblyEnd(jac, MAT_FLUSH_ASSEMBLY));

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

720:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
721:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));

723:   PetscCall(DMDAVecRestoreArray(cda, C, &c));
724:   PetscFunctionReturn(PETSC_SUCCESS);
725: }

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

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

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

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

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

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

826:   PetscFunctionBegin;
827:   PetscCall(SNESNGSGetSweeps(snes, &sweeps));
828:   PetscCall(SNESNGSGetTolerances(snes, &atol, &rtol, &stol, &its));

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

842:   PetscCall(DMGetCoordinateDM(da, &cda));
843:   PetscCall(DMGetCoordinatesLocal(da, &C));
844:   PetscCall(DMDAVecGetArray(cda, C, &c));
845:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
846:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));

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

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

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

943: PetscErrorCode InitialGuess(DM da, AppCtx *user, Vec X)
944: {
945:   PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
946:   PetscInt mx, my, mz;
947:   Field ***x;

949:   PetscFunctionBegin;
950:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
951:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
952:   PetscCall(DMDAVecGetArray(da, X, &x));

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

974: PetscErrorCode FormRHS(DM da, AppCtx *user, Vec X)
975: {
976:   PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
977:   PetscInt mx, my, mz;
978:   Field ***x;

980:   PetscFunctionBegin;
981:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
982:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
983:   PetscCall(DMDAVecGetArray(da, X, &x));

985:   for (k = zs; k < zs + zm; k++) {
986:     for (j = ys; j < ys + ym; j++) {
987:       for (i = xs; i < xs + xm; i++) {
988:         x[k][j][i][0] = 0.;
989:         x[k][j][i][1] = 0.;
990:         x[k][j][i][2] = 0.;
991:         if (i == (mx - 1) / 2 && j == (my - 1)) x[k][j][i][1] = user->ploading / (mz - 1);
992:       }
993:     }
994:   }
995:   PetscCall(DMDAVecRestoreArray(da, X, &x));
996:   PetscFunctionReturn(PETSC_SUCCESS);
997: }

999: PetscErrorCode DisplayLine(SNES snes, Vec X)
1000: {
1001:   PetscInt      r, i, j = 0, k = 0, xs, xm, ys, ym, zs, zm, mx, my, mz;
1002:   Field      ***x;
1003:   CoordField ***c;
1004:   DM            da, cda;
1005:   Vec           C;
1006:   PetscMPIInt   size, rank;

1008:   PetscFunctionBegin;
1009:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1010:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1011:   PetscCall(SNESGetDM(snes, &da));
1012:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1013:   PetscCall(DMGetCoordinateDM(da, &cda));
1014:   PetscCall(DMGetCoordinates(da, &C));
1015:   j = my / 2;
1016:   k = mz / 2;
1017:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
1018:   PetscCall(DMDAVecGetArray(da, X, &x));
1019:   PetscCall(DMDAVecGetArray(cda, C, &c));
1020:   for (r = 0; r < size; r++) {
1021:     if (rank == r) {
1022:       if (j >= ys && j < ys + ym && k >= zs && k < zs + zm) {
1023:         for (i = xs; i < xs + xm; i++) {
1024:           PetscCall(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])));
1025:         }
1026:       }
1027:     }
1028:     PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
1029:   }
1030:   PetscCall(DMDAVecRestoreArray(da, X, &x));
1031:   PetscCall(DMDAVecRestoreArray(cda, C, &c));
1032:   PetscFunctionReturn(PETSC_SUCCESS);
1033: }

1035: /*TEST

1037:    test:
1038:       nsize: 2
1039:       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
1040:       requires: !single
1041:       timeoutfactor: 3

1043:    test:
1044:       suffix: 2
1045:       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
1046:       requires: !single

1048:    test:
1049:       suffix: 3
1050:       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
1051:       requires: !single

1053: TEST*/