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:     Homogeneous 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";

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

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

114:   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));
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:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-arc", &user.arc, NULL));
125:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, &muflg));
126:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-lambda", &user.lambda, &lambdaflg));
127:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-rad", &user.rad, NULL));
128:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-height", &user.height, NULL));
129:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-width", &user.width, NULL));
130:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-loading", &user.loading, NULL));
131:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-ploading", &user.ploading, NULL));
132:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-poisson", &poisson, &poissonflg));
133:   PetscCall(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:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view", &view, NULL));
140:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_line", &viewline, NULL));

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

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

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

157:   PetscCall(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) PetscCall(DisplayLine(snes, x));

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

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

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

184:   PetscCall(VecDestroy(&x));
185:   PetscCall(VecDestroy(&b));
186:   PetscCall(DMDestroy(&da));
187:   PetscCall(SNESDestroy(&snes));
188:   PetscCall(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:   PetscFunctionBegin;
371:   /* construct the basis function values and derivatives */
372:   for (k = 0; k < NB; k++) {
373:     for (j = 0; j < NB; j++) {
374:       for (i = 0; i < NB; i++) {
375:         /* loop over the quadrature points */
376:         for (kk = 0; kk < NQ; kk++) {
377:           for (jj = 0; jj < NQ; jj++) {
378:             for (ii = 0; ii < NQ; ii++) {
379:               PetscInt idx = ii + NQ * jj + NQ * NQ * kk + NEQ * i + NEQ * NB * j + NEQ * NB * NB * k;
380:               bx           = pts[ii];
381:               by           = pts[jj];
382:               bz           = pts[kk];
383:               dbx          = 1.;
384:               dby          = 1.;
385:               dbz          = 1.;
386:               if (i == 0) {
387:                 bx  = 1. - bx;
388:                 dbx = -1;
389:               }
390:               if (j == 0) {
391:                 by  = 1. - by;
392:                 dby = -1;
393:               }
394:               if (k == 0) {
395:                 bz  = 1. - bz;
396:                 dbz = -1;
397:               }
398:               vals[idx]         = bx * by * bz;
399:               grad[3 * idx + 0] = dbx * by * bz;
400:               grad[3 * idx + 1] = dby * bx * bz;
401:               grad[3 * idx + 2] = dbz * bx * by;
402:             }
403:           }
404:         }
405:       }
406:     }
407:   }
408:   PetscFunctionReturn(PETSC_SUCCESS);
409: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1036: /*TEST

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

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

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

1054: TEST*/