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.  This example
     37:     also demonstrates the use of the arc length continuation method NEWTONAL, which avoids the numerical difficulties
     38:     of the snap-through via tracing the equilibrium path through load increments.

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

     43:   ------------------------------------------------------------------------F*/

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

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

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

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

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

     67: typedef PetscScalar JacField[9];

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

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

     86: PetscErrorCode InitialGuess(DM, AppCtx *, Vec);
     87: PetscErrorCode FormRHS(DM, AppCtx *, Vec);
     88: PetscErrorCode FormCoordinates(DM, AppCtx *);
     89: PetscErrorCode TangentLoad(SNES, Vec, Vec, void *);

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

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

    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;
    123:   user.load_factor = 1.0;

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

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

    147:   PetscCall(DMSetApplicationContext(da, &user));
    148:   PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode (*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user));
    149:   PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)FormJacobianLocal, &user));
    150:   PetscCall(SNESSetFromOptions(snes));
    151:   PetscCall(SNESNewtonALSetFunction(snes, TangentLoad, &user));
    152:   PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONAL, &alflg));
    153:   if (alflg) user.load_factor = 0.0;

    155:   PetscCall(FormCoordinates(da, &user));

    157:   PetscCall(DMCreateGlobalVector(da, &x));
    158:   PetscCall(DMCreateGlobalVector(da, &b));
    159:   PetscCall(InitialGuess(da, &user, x));
    160:   PetscCall(FormRHS(da, &user, b));

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

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

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

    170:   PetscCall(SNESGetIterationNumber(snes, &its));
    171:   PetscCall(PetscPrintf(comm, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
    172:   PetscCall(SNESGetSolution(snes, &X));
    173:   /* show a cross-section of the final state */
    174:   if (viewline) PetscCall(DisplayLine(snes, X));

    176:   if (view) {
    177:     PetscViewer viewer;
    178:     Vec         coords;
    179:     PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename, FILE_MODE_WRITE, &viewer));
    180:     PetscCall(VecView(x, viewer));
    181:     PetscCall(PetscViewerDestroy(&viewer));
    182:     PetscCall(DMGetCoordinates(da, &coords));
    183:     PetscCall(VecAXPY(coords, 1.0, x));
    184:     PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename_def, FILE_MODE_WRITE, &viewer));
    185:     PetscCall(VecView(x, viewer));
    186:     PetscCall(PetscViewerDestroy(&viewer));
    187:   }

    189:   PetscCall(VecDestroy(&x));
    190:   PetscCall(VecDestroy(&b));
    191:   PetscCall(DMDestroy(&da));
    192:   PetscCall(SNESDestroy(&snes));
    193:   PetscCall(PetscFinalize());
    194:   return 0;
    195: }

    197: PetscInt OnBoundary(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz)
    198: {
    199:   if ((i == 0 || i == mx - 1) && j == my - 1) return 1;
    200:   return 0;
    201: }

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

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

    253: void TensorTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
    254: {
    255:   PetscInt i, j, m;
    256:   for (i = 0; i < 3; i++) {
    257:     for (j = 0; j < 3; j++) {
    258:       c[i + 3 * j] = 0;
    259:       for (m = 0; m < 3; m++) c[i + 3 * j] += a[m + 3 * j] * b[i + 3 * m];
    260:     }
    261:   }
    262: }

    264: void TensorTransposeTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
    265: {
    266:   PetscInt i, j, m;
    267:   for (i = 0; i < 3; i++) {
    268:     for (j = 0; j < 3; j++) {
    269:       c[i + 3 * j] = 0;
    270:       for (m = 0; m < 3; m++) c[i + 3 * j] += a[3 * m + j] * b[i + 3 * m];
    271:     }
    272:   }
    273: }

    275: void TensorVector(PetscScalar *rot, PetscScalar *vec, PetscScalar *tvec)
    276: {
    277:   tvec[0] = rot[0] * vec[0] + rot[1] * vec[1] + rot[2] * vec[2];
    278:   tvec[1] = rot[3] * vec[0] + rot[4] * vec[1] + rot[5] * vec[2];
    279:   tvec[2] = rot[6] * vec[0] + rot[7] * vec[1] + rot[8] * vec[2];
    280: }

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

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

    325: void LagrangeGreenStrain(PetscScalar *F, PetscScalar *E)
    326: {
    327:   PetscInt i, j, m;
    328:   for (i = 0; i < 3; i++) {
    329:     for (j = 0; j < 3; j++) {
    330:       E[i + 3 * j] = 0;
    331:       for (m = 0; m < 3; m++) E[i + 3 * j] += 0.5 * F[3 * m + j] * F[i + 3 * m];
    332:     }
    333:   }
    334:   for (i = 0; i < 3; i++) E[i + 3 * i] -= 0.5;
    335: }

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

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

    370: PetscErrorCode FormElements()
    371: {
    372:   PetscInt  i, j, k, ii, jj, kk;
    373:   PetscReal bx, by, bz, dbx, dby, dbz;

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

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

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

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

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

    560: void ApplyBCsElement(PetscInt mx, PetscInt my, PetscInt mz, PetscInt i, PetscInt j, PetscInt k, PetscScalar *jacobian)
    561: {
    562:   PetscInt ii, jj, kk, ll, ei, ej, ek, el;
    563:   for (kk = 0; kk < NB; kk++) {
    564:     for (jj = 0; jj < NB; jj++) {
    565:       for (ii = 0; ii < NB; ii++) {
    566:         for (ll = 0; ll < 3; ll++) {
    567:           PetscInt tridx = ll + 3 * (ii + jj * NB + kk * NB * NB);
    568:           for (ek = 0; ek < NB; ek++) {
    569:             for (ej = 0; ej < NB; ej++) {
    570:               for (ei = 0; ei < NB; ei++) {
    571:                 for (el = 0; el < 3; el++) {
    572:                   if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz) || OnBoundary(i + ei, j + ej, k + ek, mx, my, mz)) {
    573:                     PetscInt teidx = el + 3 * (ei + ej * NB + ek * NB * NB);
    574:                     if (teidx == tridx) {
    575:                       jacobian[tridx + NPB * teidx] = 1.;
    576:                     } else {
    577:                       jacobian[tridx + NPB * teidx] = 0.;
    578:                     }
    579:                   }
    580:                 }
    581:               }
    582:             }
    583:           }
    584:         }
    585:       }
    586:     }
    587:   }
    588: }

    590: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field ***x, Mat jacpre, Mat jac, void *ptr)
    591: {
    592:   /* values for each basis function at each quadrature point */
    593:   AppCtx     *user = (AppCtx *)ptr;
    594:   PetscInt    i, j, k, m, l;
    595:   PetscInt    ii, jj, kk;
    596:   PetscScalar ej[NPB * NPB];
    597:   PetscScalar vals[NPB * NPB];
    598:   Field       ex[NEB];
    599:   CoordField  ec[NEB];

    601:   PetscInt      xs = info->xs, ys = info->ys, zs = info->zs;
    602:   PetscInt      xm = info->xm, ym = info->ym, zm = info->zm;
    603:   PetscInt      xes, yes, zes, xee, yee, zee;
    604:   PetscInt      mx = info->mx, my = info->my, mz = info->mz;
    605:   DM            cda;
    606:   CoordField ***c;
    607:   Vec           C;
    608:   PetscInt      nrows;
    609:   MatStencil    col[NPB], row[NPB];
    610:   PetscScalar   v[9];

    612:   PetscFunctionBegin;
    613:   PetscCall(DMGetCoordinateDM(info->da, &cda));
    614:   PetscCall(DMGetCoordinatesLocal(info->da, &C));
    615:   PetscCall(DMDAVecGetArray(cda, C, &c));
    616:   PetscCall(MatScale(jac, 0.0));

    618:   xes = xs;
    619:   yes = ys;
    620:   zes = zs;
    621:   xee = xs + xm;
    622:   yee = ys + ym;
    623:   zee = zs + zm;
    624:   if (xs > 0) xes = xs - 1;
    625:   if (ys > 0) yes = ys - 1;
    626:   if (zs > 0) zes = zs - 1;
    627:   if (xs + xm == mx) xee = xs + xm - 1;
    628:   if (ys + ym == my) yee = ys + ym - 1;
    629:   if (zs + zm == mz) zee = zs + zm - 1;
    630:   for (k = zes; k < zee; k++) {
    631:     for (j = yes; j < yee; j++) {
    632:       for (i = xes; i < xee; i++) {
    633:         GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
    634:         FormElementJacobian(ex, ec, NULL, NULL, ej, user);
    635:         ApplyBCsElement(mx, my, mz, i, j, k, ej);
    636:         nrows = 0.;
    637:         for (kk = 0; kk < NB; kk++) {
    638:           for (jj = 0; jj < NB; jj++) {
    639:             for (ii = 0; ii < NB; ii++) {
    640:               PetscInt idx = ii + jj * 2 + kk * 4;
    641:               for (m = 0; m < 3; m++) {
    642:                 col[3 * idx + m].i = i + ii;
    643:                 col[3 * idx + m].j = j + jj;
    644:                 col[3 * idx + m].k = k + kk;
    645:                 col[3 * idx + m].c = m;
    646:                 if (i + ii >= xs && i + ii < xm + xs && j + jj >= ys && j + jj < ys + ym && k + kk >= zs && k + kk < zs + zm) {
    647:                   row[nrows].i = i + ii;
    648:                   row[nrows].j = j + jj;
    649:                   row[nrows].k = k + kk;
    650:                   row[nrows].c = m;
    651:                   for (l = 0; l < NPB; l++) vals[NPB * nrows + l] = ej[NPB * (3 * idx + m) + l];
    652:                   nrows++;
    653:                 }
    654:               }
    655:             }
    656:           }
    657:         }
    658:         PetscCall(MatSetValuesStencil(jac, nrows, row, NPB, col, vals, ADD_VALUES));
    659:       }
    660:     }
    661:   }

    663:   PetscCall(MatAssemblyBegin(jac, MAT_FLUSH_ASSEMBLY));
    664:   PetscCall(MatAssemblyEnd(jac, MAT_FLUSH_ASSEMBLY));

    666:   /* set the diagonal */
    667:   v[0] = 1.;
    668:   v[1] = 0.;
    669:   v[2] = 0.;
    670:   v[3] = 0.;
    671:   v[4] = 1.;
    672:   v[5] = 0.;
    673:   v[6] = 0.;
    674:   v[7] = 0.;
    675:   v[8] = 1.;
    676:   for (k = zs; k < zs + zm; k++) {
    677:     for (j = ys; j < ys + ym; j++) {
    678:       for (i = xs; i < xs + xm; i++) {
    679:         if (OnBoundary(i, j, k, mx, my, mz)) {
    680:           for (m = 0; m < 3; m++) {
    681:             col[m].i = i;
    682:             col[m].j = j;
    683:             col[m].k = k;
    684:             col[m].c = m;
    685:           }
    686:           PetscCall(MatSetValuesStencil(jac, 3, col, 3, col, v, INSERT_VALUES));
    687:         }
    688:       }
    689:     }
    690:   }

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

    695:   PetscCall(DMDAVecRestoreArray(cda, C, &c));
    696:   PetscFunctionReturn(PETSC_SUCCESS);
    697: }

    699: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field ***x, Field ***f, void *ptr)
    700: {
    701:   /* values for each basis function at each quadrature point */
    702:   AppCtx  *user = (AppCtx *)ptr;
    703:   PetscInt i, j, k, l;
    704:   PetscInt ii, jj, kk;

    706:   Field      ef[NEB];
    707:   Field      ex[NEB];
    708:   CoordField ec[NEB];

    710:   PetscInt      xs = info->xs, ys = info->ys, zs = info->zs;
    711:   PetscInt      xm = info->xm, ym = info->ym, zm = info->zm;
    712:   PetscInt      xes, yes, zes, xee, yee, zee;
    713:   PetscInt      mx = info->mx, my = info->my, mz = info->mz;
    714:   DM            cda;
    715:   CoordField ***c;
    716:   Vec           C;

    718:   PetscFunctionBegin;
    719:   PetscCall(DMGetCoordinateDM(info->da, &cda));
    720:   PetscCall(DMGetCoordinatesLocal(info->da, &C));
    721:   PetscCall(DMDAVecGetArray(cda, C, &c));
    722:   PetscCall(DMDAGetInfo(info->da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    723:   PetscCall(DMDAGetCorners(info->da, &xs, &ys, &zs, &xm, &ym, &zm));

    725:   /* loop over elements */
    726:   for (k = zs; k < zs + zm; k++) {
    727:     for (j = ys; j < ys + ym; j++) {
    728:       for (i = xs; i < xs + xm; i++) {
    729:         for (l = 0; l < 3; l++) f[k][j][i][l] = 0.;
    730:       }
    731:     }
    732:   }
    733:   /* element starts and ends */
    734:   xes = xs;
    735:   yes = ys;
    736:   zes = zs;
    737:   xee = xs + xm;
    738:   yee = ys + ym;
    739:   zee = zs + zm;
    740:   if (xs > 0) xes = xs - 1;
    741:   if (ys > 0) yes = ys - 1;
    742:   if (zs > 0) zes = zs - 1;
    743:   if (xs + xm == mx) xee = xs + xm - 1;
    744:   if (ys + ym == my) yee = ys + ym - 1;
    745:   if (zs + zm == mz) zee = zs + zm - 1;
    746:   for (k = zes; k < zee; k++) {
    747:     for (j = yes; j < yee; j++) {
    748:       for (i = xes; i < xee; i++) {
    749:         GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
    750:         FormElementJacobian(ex, ec, ef, NULL, NULL, user);
    751:         /* put this element's additions into the residuals */
    752:         for (kk = 0; kk < NB; kk++) {
    753:           for (jj = 0; jj < NB; jj++) {
    754:             for (ii = 0; ii < NB; ii++) {
    755:               PetscInt idx = ii + jj * NB + kk * NB * NB;
    756:               if (k + kk >= zs && j + jj >= ys && i + ii >= xs && k + kk < zs + zm && j + jj < ys + ym && i + ii < xs + xm) {
    757:                 if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
    758:                   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];
    759:                 } else {
    760:                   for (l = 0; l < 3; l++) f[k + kk][j + jj][i + ii][l] += ef[idx][l];
    761:                 }
    762:               }
    763:             }
    764:           }
    765:         }
    766:       }
    767:     }
    768:   }
    769:   PetscCall(DMDAVecRestoreArray(cda, C, &c));
    770:   PetscFunctionReturn(PETSC_SUCCESS);
    771: }

    773: PetscErrorCode TangentLoad(SNES snes, Vec X, Vec Q, void *ptr)
    774: {
    775:   /* values for each basis function at each quadrature point */
    776:   AppCtx  *user = (AppCtx *)ptr;
    777:   PetscInt xs, ys, zs;
    778:   PetscInt xm, ym, zm;
    779:   PetscInt mx, my, mz;
    780:   DM       da;
    781:   Vec      Xl, Ql;
    782:   Field ***x, ***q;
    783:   PetscInt i, j, k, l;
    784:   PetscInt ii, jj, kk;

    786:   Field      eq[NEB];
    787:   Field      ex[NEB];
    788:   CoordField ec[NEB];

    790:   PetscInt      xes, yes, zes, xee, yee, zee;
    791:   DM            cda;
    792:   CoordField ***c;
    793:   Vec           C;

    795:   PetscFunctionBegin;
    796:   /* update user context with current load parameter */
    797:   PetscCall(SNESNewtonALGetLoadParameter(snes, &user->load_factor));

    799:   PetscCall(SNESGetDM(snes, &da));
    800:   PetscCall(DMGetLocalVector(da, &Xl));
    801:   PetscCall(DMGetLocalVector(da, &Ql));
    802:   PetscCall(DMGlobalToLocal(da, X, INSERT_VALUES, Xl));

    804:   PetscCall(DMDAVecGetArray(da, Xl, &x));
    805:   PetscCall(DMDAVecGetArray(da, Ql, &q));

    807:   PetscCall(DMGetCoordinateDM(da, &cda));
    808:   PetscCall(DMGetCoordinatesLocal(da, &C));
    809:   PetscCall(DMDAVecGetArray(cda, C, &c));
    810:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    811:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));

    813:   /* loop over elements */
    814:   for (k = zs; k < zs + zm; k++) {
    815:     for (j = ys; j < ys + ym; j++) {
    816:       for (i = xs; i < xs + xm; i++) {
    817:         for (l = 0; l < 3; l++) q[k][j][i][l] = 0.;
    818:       }
    819:     }
    820:   }
    821:   /* element starts and ends */
    822:   xes = xs;
    823:   yes = ys;
    824:   zes = zs;
    825:   xee = xs + xm;
    826:   yee = ys + ym;
    827:   zee = zs + zm;
    828:   if (xs > 0) xes = xs - 1;
    829:   if (ys > 0) yes = ys - 1;
    830:   if (zs > 0) zes = zs - 1;
    831:   if (xs + xm == mx) xee = xs + xm - 1;
    832:   if (ys + ym == my) yee = ys + ym - 1;
    833:   if (zs + zm == mz) zee = zs + zm - 1;
    834:   for (k = zes; k < zee; k++) {
    835:     for (j = yes; j < yee; j++) {
    836:       for (i = xes; i < xee; i++) {
    837:         GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
    838:         FormElementJacobian(ex, ec, NULL, eq, NULL, user);
    839:         /* put this element's additions into the residuals */
    840:         for (kk = 0; kk < NB; kk++) {
    841:           for (jj = 0; jj < NB; jj++) {
    842:             for (ii = 0; ii < NB; ii++) {
    843:               PetscInt idx = ii + jj * NB + kk * NB * NB;
    844:               if (k + kk >= zs && j + jj >= ys && i + ii >= xs && k + kk < zs + zm && j + jj < ys + ym && i + ii < xs + xm) {
    845:                 if (!OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
    846:                   for (l = 0; l < 3; l++) q[k + kk][j + jj][i + ii][l] += eq[idx][l];
    847:                 }
    848:               }
    849:             }
    850:           }
    851:         }
    852:       }
    853:     }
    854:   }
    855:   PetscCall(DMDAVecRestoreArray(da, Xl, &x));
    856:   PetscCall(DMDAVecRestoreArray(da, Ql, &q));
    857:   PetscCall(VecZeroEntries(Q));
    858:   PetscCall(DMLocalToGlobal(da, Ql, INSERT_VALUES, Q));
    859:   PetscCall(DMRestoreLocalVector(da, &Ql));
    860:   PetscCall(DMRestoreLocalVector(da, &Xl));
    861:   PetscCall(DMDAVecRestoreArray(cda, C, &c));
    862:   PetscFunctionReturn(PETSC_SUCCESS);
    863: }

    865: PetscErrorCode FormCoordinates(DM da, AppCtx *user)
    866: {
    867:   Vec           coords;
    868:   DM            cda;
    869:   PetscInt      mx, my, mz;
    870:   PetscInt      i, j, k, xs, ys, zs, xm, ym, zm;
    871:   CoordField ***x;

    873:   PetscFunctionBegin;
    874:   PetscCall(DMGetCoordinateDM(da, &cda));
    875:   PetscCall(DMCreateGlobalVector(cda, &coords));
    876:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    877:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
    878:   PetscCall(DMDAVecGetArray(da, coords, &x));
    879:   for (k = zs; k < zs + zm; k++) {
    880:     for (j = ys; j < ys + ym; j++) {
    881:       for (i = xs; i < xs + xm; i++) {
    882:         PetscReal cx  = ((PetscReal)i) / ((PetscReal)(mx - 1));
    883:         PetscReal cy  = ((PetscReal)j) / ((PetscReal)(my - 1));
    884:         PetscReal cz  = ((PetscReal)k) / ((PetscReal)(mz - 1));
    885:         PetscReal rad = user->rad + cy * user->height;
    886:         PetscReal ang = (cx - 0.5) * user->arc;
    887:         x[k][j][i][0] = rad * PetscSinReal(ang);
    888:         x[k][j][i][1] = rad * PetscCosReal(ang) - (user->rad + 0.5 * user->height) * PetscCosReal(-0.5 * user->arc);
    889:         x[k][j][i][2] = user->width * (cz - 0.5);
    890:       }
    891:     }
    892:   }
    893:   PetscCall(DMDAVecRestoreArray(da, coords, &x));
    894:   PetscCall(DMSetCoordinates(da, coords));
    895:   PetscCall(VecDestroy(&coords));
    896:   PetscFunctionReturn(PETSC_SUCCESS);
    897: }

    899: PetscErrorCode InitialGuess(DM da, AppCtx *user, Vec X)
    900: {
    901:   PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
    902:   PetscInt mx, my, mz;
    903:   Field ***x;

    905:   PetscFunctionBegin;
    906:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
    907:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    908:   PetscCall(DMDAVecGetArray(da, X, &x));

    910:   for (k = zs; k < zs + zm; k++) {
    911:     for (j = ys; j < ys + ym; j++) {
    912:       for (i = xs; i < xs + xm; i++) {
    913:         /* reference coordinates */
    914:         PetscReal p_x = ((PetscReal)i) / ((PetscReal)(mx - 1));
    915:         PetscReal p_y = ((PetscReal)j) / ((PetscReal)(my - 1));
    916:         PetscReal p_z = ((PetscReal)k) / ((PetscReal)(mz - 1));
    917:         PetscReal o_x = p_x;
    918:         PetscReal o_y = p_y;
    919:         PetscReal o_z = p_z;
    920:         x[k][j][i][0] = o_x - p_x;
    921:         x[k][j][i][1] = o_y - p_y;
    922:         x[k][j][i][2] = o_z - p_z;
    923:       }
    924:     }
    925:   }
    926:   PetscCall(DMDAVecRestoreArray(da, X, &x));
    927:   PetscFunctionReturn(PETSC_SUCCESS);
    928: }

    930: PetscErrorCode FormRHS(DM da, AppCtx *user, Vec X)
    931: {
    932:   PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
    933:   PetscInt mx, my, mz;
    934:   Field ***x;

    936:   PetscFunctionBegin;
    937:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
    938:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    939:   PetscCall(DMDAVecGetArray(da, X, &x));

    941:   for (k = zs; k < zs + zm; k++) {
    942:     for (j = ys; j < ys + ym; j++) {
    943:       for (i = xs; i < xs + xm; i++) {
    944:         x[k][j][i][0] = 0.;
    945:         x[k][j][i][1] = 0.;
    946:         x[k][j][i][2] = 0.;
    947:         if (i == (mx - 1) / 2 && j == (my - 1)) x[k][j][i][1] = user->ploading / (mz - 1);
    948:       }
    949:     }
    950:   }
    951:   PetscCall(DMDAVecRestoreArray(da, X, &x));
    952:   PetscFunctionReturn(PETSC_SUCCESS);
    953: }

    955: PetscErrorCode DisplayLine(SNES snes, Vec X)
    956: {
    957:   PetscInt      r, i, j = 0, k = 0, xs, xm, ys, ym, zs, zm, mx, my, mz;
    958:   Field      ***x;
    959:   CoordField ***c;
    960:   DM            da, cda;
    961:   Vec           C;
    962:   PetscMPIInt   size, rank;

    964:   PetscFunctionBegin;
    965:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    966:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    967:   PetscCall(SNESGetDM(snes, &da));
    968:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    969:   PetscCall(DMGetCoordinateDM(da, &cda));
    970:   PetscCall(DMGetCoordinates(da, &C));
    971:   j = my / 2;
    972:   k = mz / 2;
    973:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
    974:   PetscCall(DMDAVecGetArray(da, X, &x));
    975:   PetscCall(DMDAVecGetArray(cda, C, &c));
    976:   for (r = 0; r < size; r++) {
    977:     if (rank == r) {
    978:       if (j >= ys && j < ys + ym && k >= zs && k < zs + zm) {
    979:         for (i = xs; i < xs + xm; i++) {
    980:           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])));
    981:         }
    982:       }
    983:     }
    984:     PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
    985:   }
    986:   PetscCall(DMDAVecRestoreArray(da, X, &x));
    987:   PetscCall(DMDAVecRestoreArray(cda, C, &c));
    988:   PetscFunctionReturn(PETSC_SUCCESS);
    989: }

    991: /*TEST

    993:    test:
    994:       nsize: 2
    995:       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
    996:       requires: !single
    997:       timeoutfactor: 3

    999:    test:
    1000:       suffix: 2
    1001:       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 -npc_fas_levels_snes_linesearch_maxlambda 2.0
    1002:       requires: !single

    1004:    test:
    1005:       suffix: 3
    1006:       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 -npc_sub_ksp_type preonly -npc_sub_pc_type lu
    1007:       requires: !single

    1009:    test:
    1010:       suffix: 4
    1011:       args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading -1. -loading -1. -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -snes_type newtonal -snes_newtonal_step_size 30 -ksp_rtol 1e-4 -info
    1012:       requires: !single defined(PETSC_USE_INFO)
    1013:       filter: grep -h -e "KSP Residual norm" -e "SNES Function norm" -e "Number of SNES iterations" -e "mu:"

    1015:    test:
    1016:       suffix: 5
    1017:       args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading -1. -loading -1. -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -snes_type newtonal -snes_newtonal_step_size 30 -snes_newtonal_correction_type normal -ksp_rtol 1e-4
    1018:       requires: !single

    1020:    test:
    1021:       suffix: 6
    1022:       args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading -0.5 -loading -1 -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -snes_type newtonal -snes_newtonal_step_size 0.5 -snes_newtonal_max_continuation_steps 3 -snes_newtonal_scale_rhs false -snes_newtonal_lambda_min -0.1 -ksp_rtol 1e-4
    1023:       requires: !single

    1025: TEST*/