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 ArcLengthScaling(DM, AppCtx *, Vec);
     88: PetscErrorCode FormRHS(DM, AppCtx *, Vec);
     89: PetscErrorCode FormCoordinates(DM, AppCtx *);
     90: PetscErrorCode TangentLoad(SNES, Vec, Vec, void *);

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

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

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

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

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

    149:   PetscCall(DMSetApplicationContext(da, &user));
    150:   PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunctionFn *)FormFunctionLocal, &user));
    151:   PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)FormJacobianLocal, &user));
    152:   PetscCall(SNESSetFromOptions(snes));
    153:   PetscCall(SNESNewtonALSetFunction(snes, TangentLoad, &user));
    154:   PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONAL, &alflg));
    155:   if (alflg) {
    156:     user.load_factor = 0.0;

    158:     if (al_rescale) {
    159:       Vec scaling;

    161:       PetscCall(DMCreateGlobalVector(da, &scaling));
    162:       PetscCall(ArcLengthScaling(da, &user, scaling));
    163:       PetscCall(VecPointwiseMult(scaling, scaling, scaling));
    164:       PetscCall(SNESNewtonALSetDiagonalScaling(snes, scaling));
    165:       PetscCall(VecDestroy(&scaling));
    166:     }
    167:   }

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

    171:   PetscCall(DMCreateGlobalVector(da, &x));
    172:   PetscCall(DMCreateGlobalVector(da, &b));
    173:   PetscCall(InitialGuess(da, &user, x));
    174:   PetscCall(FormRHS(da, &user, b));

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

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

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

    184:   PetscCall(SNESGetIterationNumber(snes, &its));
    185:   PetscCall(PetscPrintf(comm, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
    186:   PetscCall(SNESGetSolution(snes, &X));
    187:   /* show a cross-section of the final state */
    188:   if (viewline) PetscCall(DisplayLine(snes, X));

    190:   if (view) {
    191:     PetscViewer viewer;
    192:     Vec         coords;
    193:     PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename, FILE_MODE_WRITE, &viewer));
    194:     PetscCall(VecView(x, viewer));
    195:     PetscCall(PetscViewerDestroy(&viewer));
    196:     PetscCall(DMGetCoordinates(da, &coords));
    197:     PetscCall(VecAXPY(coords, 1.0, x));
    198:     PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename_def, FILE_MODE_WRITE, &viewer));
    199:     PetscCall(VecView(x, viewer));
    200:     PetscCall(PetscViewerDestroy(&viewer));
    201:   }

    203:   PetscCall(VecDestroy(&x));
    204:   PetscCall(VecDestroy(&b));
    205:   PetscCall(DMDestroy(&da));
    206:   PetscCall(SNESDestroy(&snes));
    207:   PetscCall(PetscFinalize());
    208:   return 0;
    209: }

    211: PetscInt OnBoundary(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz)
    212: {
    213:   if ((i == 0 || i == mx - 1) && j == my - 1) return 1;
    214:   return 0;
    215: }

    217: void BoundaryValue(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz, PetscScalar *val, AppCtx *user)
    218: {
    219:   /* reference coordinates */
    220:   PetscReal p_x = ((PetscReal)i) / ((PetscReal)(mx - 1));
    221:   PetscReal p_y = ((PetscReal)j) / ((PetscReal)(my - 1));
    222:   PetscReal p_z = ((PetscReal)k) / ((PetscReal)(mz - 1));
    223:   PetscReal o_x = p_x;
    224:   PetscReal o_y = p_y;
    225:   PetscReal o_z = p_z;
    226:   val[0]        = o_x - p_x;
    227:   val[1]        = o_y - p_y;
    228:   val[2]        = o_z - p_z;
    229: }

    231: void InvertTensor(PetscScalar *t, PetscScalar *ti, PetscReal *dett)
    232: {
    233:   const PetscScalar a   = t[0];
    234:   const PetscScalar b   = t[1];
    235:   const PetscScalar c   = t[2];
    236:   const PetscScalar d   = t[3];
    237:   const PetscScalar e   = t[4];
    238:   const PetscScalar f   = t[5];
    239:   const PetscScalar g   = t[6];
    240:   const PetscScalar h   = t[7];
    241:   const PetscScalar i   = t[8];
    242:   const PetscReal   det = PetscRealPart(a * (e * i - f * h) - b * (i * d - f * g) + c * (d * h - e * g));
    243:   const PetscReal   di  = 1. / det;
    244:   if (ti) {
    245:     const PetscScalar A  = (e * i - f * h);
    246:     const PetscScalar B  = -(d * i - f * g);
    247:     const PetscScalar C  = (d * h - e * g);
    248:     const PetscScalar D  = -(b * i - c * h);
    249:     const PetscScalar E  = (a * i - c * g);
    250:     const PetscScalar F  = -(a * h - b * g);
    251:     const PetscScalar G  = (b * f - c * e);
    252:     const PetscScalar H  = -(a * f - c * d);
    253:     const PetscScalar II = (a * e - b * d);
    254:     ti[0]                = di * A;
    255:     ti[1]                = di * D;
    256:     ti[2]                = di * G;
    257:     ti[3]                = di * B;
    258:     ti[4]                = di * E;
    259:     ti[5]                = di * H;
    260:     ti[6]                = di * C;
    261:     ti[7]                = di * F;
    262:     ti[8]                = di * II;
    263:   }
    264:   if (dett) *dett = det;
    265: }

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

    278: void TensorTransposeTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
    279: {
    280:   PetscInt i, j, m;
    281:   for (i = 0; i < 3; i++) {
    282:     for (j = 0; j < 3; j++) {
    283:       c[i + 3 * j] = 0;
    284:       for (m = 0; m < 3; m++) c[i + 3 * j] += a[3 * m + j] * b[i + 3 * m];
    285:     }
    286:   }
    287: }

    289: void TensorVector(PetscScalar *rot, PetscScalar *vec, PetscScalar *tvec)
    290: {
    291:   tvec[0] = rot[0] * vec[0] + rot[1] * vec[1] + rot[2] * vec[2];
    292:   tvec[1] = rot[3] * vec[0] + rot[4] * vec[1] + rot[5] * vec[2];
    293:   tvec[2] = rot[6] * vec[0] + rot[7] * vec[1] + rot[8] * vec[2];
    294: }

    296: void DeformationGradient(Field *ex, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *invJ, PetscScalar *F)
    297: {
    298:   PetscInt ii, jj, kk, l;
    299:   for (l = 0; l < 9; l++) F[l] = 0.;
    300:   F[0] = 1.;
    301:   F[4] = 1.;
    302:   F[8] = 1.;
    303:   /* form the deformation gradient at this basis function -- loop over element unknowns */
    304:   for (kk = 0; kk < NB; kk++) {
    305:     for (jj = 0; jj < NB; jj++) {
    306:       for (ii = 0; ii < NB; ii++) {
    307:         PetscInt    idx  = ii + jj * NB + kk * NB * NB;
    308:         PetscInt    bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
    309:         PetscScalar lgrad[3];
    310:         TensorVector(invJ, &grad[3 * bidx], lgrad);
    311:         F[0] += lgrad[0] * ex[idx][0];
    312:         F[1] += lgrad[1] * ex[idx][0];
    313:         F[2] += lgrad[2] * ex[idx][0];
    314:         F[3] += lgrad[0] * ex[idx][1];
    315:         F[4] += lgrad[1] * ex[idx][1];
    316:         F[5] += lgrad[2] * ex[idx][1];
    317:         F[6] += lgrad[0] * ex[idx][2];
    318:         F[7] += lgrad[1] * ex[idx][2];
    319:         F[8] += lgrad[2] * ex[idx][2];
    320:       }
    321:     }
    322:   }
    323: }

    325: void DeformationGradientJacobian(PetscInt qi, PetscInt qj, PetscInt qk, PetscInt ii, PetscInt jj, PetscInt kk, PetscInt fld, PetscScalar *invJ, PetscScalar *dF)
    326: {
    327:   PetscInt    l;
    328:   PetscScalar lgrad[3];
    329:   PetscInt    idx  = ii + jj * NB + kk * NB * NB;
    330:   PetscInt    bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
    331:   for (l = 0; l < 9; l++) dF[l] = 0.;
    332:   /* form the deformation gradient at this basis function -- loop over element unknowns */
    333:   TensorVector(invJ, &grad[3 * bidx], lgrad);
    334:   dF[3 * fld]     = lgrad[0];
    335:   dF[3 * fld + 1] = lgrad[1];
    336:   dF[3 * fld + 2] = lgrad[2];
    337: }

    339: void LagrangeGreenStrain(PetscScalar *F, PetscScalar *E)
    340: {
    341:   PetscInt i, j, m;
    342:   for (i = 0; i < 3; i++) {
    343:     for (j = 0; j < 3; j++) {
    344:       E[i + 3 * j] = 0;
    345:       for (m = 0; m < 3; m++) E[i + 3 * j] += 0.5 * F[3 * m + j] * F[i + 3 * m];
    346:     }
    347:   }
    348:   for (i = 0; i < 3; i++) E[i + 3 * i] -= 0.5;
    349: }

    351: void SaintVenantKirchoff(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *S)
    352: {
    353:   PetscInt    i, j;
    354:   PetscScalar E[9];
    355:   PetscScalar trE = 0;
    356:   LagrangeGreenStrain(F, E);
    357:   for (i = 0; i < 3; i++) trE += E[i + 3 * i];
    358:   for (i = 0; i < 3; i++) {
    359:     for (j = 0; j < 3; j++) {
    360:       S[i + 3 * j] = 2. * mu * E[i + 3 * j];
    361:       if (i == j) S[i + 3 * j] += trE * lambda;
    362:     }
    363:   }
    364: }

    366: void SaintVenantKirchoffJacobian(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *dF, PetscScalar *dS)
    367: {
    368:   PetscScalar FtdF[9], dE[9];
    369:   PetscInt    i, j;
    370:   PetscScalar dtrE = 0.;
    371:   TensorTransposeTensor(dF, F, dE);
    372:   TensorTransposeTensor(F, dF, FtdF);
    373:   for (i = 0; i < 9; i++) dE[i] += FtdF[i];
    374:   for (i = 0; i < 9; i++) dE[i] *= 0.5;
    375:   for (i = 0; i < 3; i++) dtrE += dE[i + 3 * i];
    376:   for (i = 0; i < 3; i++) {
    377:     for (j = 0; j < 3; j++) {
    378:       dS[i + 3 * j] = 2. * mu * dE[i + 3 * j];
    379:       if (i == j) dS[i + 3 * j] += dtrE * lambda;
    380:     }
    381:   }
    382: }

    384: PetscErrorCode FormElements(void)
    385: {
    386:   PetscInt  i, j, k, ii, jj, kk;
    387:   PetscReal bx, by, bz, dbx, dby, dbz;

    389:   PetscFunctionBeginUser;
    390:   /* construct the basis function values and derivatives */
    391:   for (k = 0; k < NB; k++) {
    392:     for (j = 0; j < NB; j++) {
    393:       for (i = 0; i < NB; i++) {
    394:         /* loop over the quadrature points */
    395:         for (kk = 0; kk < NQ; kk++) {
    396:           for (jj = 0; jj < NQ; jj++) {
    397:             for (ii = 0; ii < NQ; ii++) {
    398:               PetscInt idx = ii + NQ * jj + NQ * NQ * kk + NEQ * i + NEQ * NB * j + NEQ * NB * NB * k;
    399:               bx           = pts[ii];
    400:               by           = pts[jj];
    401:               bz           = pts[kk];
    402:               dbx          = 1.;
    403:               dby          = 1.;
    404:               dbz          = 1.;
    405:               if (i == 0) {
    406:                 bx  = 1. - bx;
    407:                 dbx = -1;
    408:               }
    409:               if (j == 0) {
    410:                 by  = 1. - by;
    411:                 dby = -1;
    412:               }
    413:               if (k == 0) {
    414:                 bz  = 1. - bz;
    415:                 dbz = -1;
    416:               }
    417:               vals[idx]         = bx * by * bz;
    418:               grad[3 * idx + 0] = dbx * by * bz;
    419:               grad[3 * idx + 1] = dby * bx * bz;
    420:               grad[3 * idx + 2] = dbz * bx * by;
    421:             }
    422:           }
    423:         }
    424:       }
    425:     }
    426:   }
    427:   PetscFunctionReturn(PETSC_SUCCESS);
    428: }

    430: void GatherElementData(PetscInt mx, PetscInt my, PetscInt mz, Field ***x, CoordField ***c, PetscInt i, PetscInt j, PetscInt k, Field *ex, CoordField *ec, AppCtx *user)
    431: {
    432:   PetscInt m;
    433:   PetscInt ii, jj, kk;
    434:   /* gather the data -- loop over element unknowns */
    435:   for (kk = 0; kk < NB; kk++) {
    436:     for (jj = 0; jj < NB; jj++) {
    437:       for (ii = 0; ii < NB; ii++) {
    438:         PetscInt idx = ii + jj * NB + kk * NB * NB;
    439:         /* decouple the boundary nodes for the displacement variables */
    440:         if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
    441:           BoundaryValue(i + ii, j + jj, k + kk, mx, my, mz, ex[idx], user);
    442:         } else {
    443:           for (m = 0; m < 3; m++) ex[idx][m] = x[k + kk][j + jj][i + ii][m];
    444:         }
    445:         for (m = 0; m < 3; m++) ec[idx][m] = c[k + kk][j + jj][i + ii][m];
    446:       }
    447:     }
    448:   }
    449: }

    451: void QuadraturePointGeometricJacobian(CoordField *ec, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *J)
    452: {
    453:   PetscInt ii, jj, kk;
    454:   /* construct the gradient at the given quadrature point named by i,j,k */
    455:   for (ii = 0; ii < 9; ii++) J[ii] = 0;
    456:   for (kk = 0; kk < NB; kk++) {
    457:     for (jj = 0; jj < NB; jj++) {
    458:       for (ii = 0; ii < NB; ii++) {
    459:         PetscInt idx  = ii + jj * NB + kk * NB * NB;
    460:         PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
    461:         J[0] += grad[3 * bidx + 0] * ec[idx][0];
    462:         J[1] += grad[3 * bidx + 1] * ec[idx][0];
    463:         J[2] += grad[3 * bidx + 2] * ec[idx][0];
    464:         J[3] += grad[3 * bidx + 0] * ec[idx][1];
    465:         J[4] += grad[3 * bidx + 1] * ec[idx][1];
    466:         J[5] += grad[3 * bidx + 2] * ec[idx][1];
    467:         J[6] += grad[3 * bidx + 0] * ec[idx][2];
    468:         J[7] += grad[3 * bidx + 1] * ec[idx][2];
    469:         J[8] += grad[3 * bidx + 2] * ec[idx][2];
    470:       }
    471:     }
    472:   }
    473: }

    475: void FormElementJacobian(Field *ex, CoordField *ec, Field *ef, Field *eq, PetscScalar *ej, AppCtx *user)
    476: {
    477:   PetscReal   vol;
    478:   PetscScalar J[9];
    479:   PetscScalar invJ[9];
    480:   PetscScalar F[9], S[9], dF[9], dS[9], dFS[9], FdS[9], FS[9];
    481:   PetscReal   scl;
    482:   PetscInt    i, j, k, l, ii, jj, kk, ll, qi, qj, qk, m;

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

    574: void ApplyBCsElement(PetscInt mx, PetscInt my, PetscInt mz, PetscInt i, PetscInt j, PetscInt k, PetscScalar *jacobian)
    575: {
    576:   PetscInt ii, jj, kk, ll, ei, ej, ek, el;
    577:   for (kk = 0; kk < NB; kk++) {
    578:     for (jj = 0; jj < NB; jj++) {
    579:       for (ii = 0; ii < NB; ii++) {
    580:         for (ll = 0; ll < 3; ll++) {
    581:           PetscInt tridx = ll + 3 * (ii + jj * NB + kk * NB * NB);
    582:           for (ek = 0; ek < NB; ek++) {
    583:             for (ej = 0; ej < NB; ej++) {
    584:               for (ei = 0; ei < NB; ei++) {
    585:                 for (el = 0; el < 3; el++) {
    586:                   if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz) || OnBoundary(i + ei, j + ej, k + ek, mx, my, mz)) {
    587:                     PetscInt teidx = el + 3 * (ei + ej * NB + ek * NB * NB);
    588:                     if (teidx == tridx) {
    589:                       jacobian[tridx + NPB * teidx] = 1.;
    590:                     } else {
    591:                       jacobian[tridx + NPB * teidx] = 0.;
    592:                     }
    593:                   }
    594:                 }
    595:               }
    596:             }
    597:           }
    598:         }
    599:       }
    600:     }
    601:   }
    602: }

    604: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field ***x, Mat jacpre, Mat jac, void *ptr)
    605: {
    606:   /* values for each basis function at each quadrature point */
    607:   AppCtx     *user = (AppCtx *)ptr;
    608:   PetscInt    i, j, k, m, l;
    609:   PetscInt    ii, jj, kk;
    610:   PetscScalar ej[NPB * NPB];
    611:   PetscScalar vals[NPB * NPB];
    612:   Field       ex[NEB];
    613:   CoordField  ec[NEB];

    615:   PetscInt      xs = info->xs, ys = info->ys, zs = info->zs;
    616:   PetscInt      xm = info->xm, ym = info->ym, zm = info->zm;
    617:   PetscInt      xes, yes, zes, xee, yee, zee;
    618:   PetscInt      mx = info->mx, my = info->my, mz = info->mz;
    619:   DM            cda;
    620:   CoordField ***c;
    621:   Vec           C;
    622:   PetscInt      nrows;
    623:   MatStencil    col[NPB], row[NPB];
    624:   PetscScalar   v[9];

    626:   PetscFunctionBeginUser;
    627:   PetscCall(DMGetCoordinateDM(info->da, &cda));
    628:   PetscCall(DMGetCoordinatesLocal(info->da, &C));
    629:   PetscCall(DMDAVecGetArray(cda, C, &c));
    630:   PetscCall(MatScale(jac, 0.0));

    632:   xes = xs;
    633:   yes = ys;
    634:   zes = zs;
    635:   xee = xs + xm;
    636:   yee = ys + ym;
    637:   zee = zs + zm;
    638:   if (xs > 0) xes = xs - 1;
    639:   if (ys > 0) yes = ys - 1;
    640:   if (zs > 0) zes = zs - 1;
    641:   if (xs + xm == mx) xee = xs + xm - 1;
    642:   if (ys + ym == my) yee = ys + ym - 1;
    643:   if (zs + zm == mz) zee = zs + zm - 1;
    644:   for (k = zes; k < zee; k++) {
    645:     for (j = yes; j < yee; j++) {
    646:       for (i = xes; i < xee; i++) {
    647:         GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
    648:         FormElementJacobian(ex, ec, NULL, NULL, ej, user);
    649:         ApplyBCsElement(mx, my, mz, i, j, k, ej);
    650:         nrows = 0.;
    651:         for (kk = 0; kk < NB; kk++) {
    652:           for (jj = 0; jj < NB; jj++) {
    653:             for (ii = 0; ii < NB; ii++) {
    654:               PetscInt idx = ii + jj * 2 + kk * 4;
    655:               for (m = 0; m < 3; m++) {
    656:                 col[3 * idx + m].i = i + ii;
    657:                 col[3 * idx + m].j = j + jj;
    658:                 col[3 * idx + m].k = k + kk;
    659:                 col[3 * idx + m].c = m;
    660:                 if (i + ii >= xs && i + ii < xm + xs && j + jj >= ys && j + jj < ys + ym && k + kk >= zs && k + kk < zs + zm) {
    661:                   row[nrows].i = i + ii;
    662:                   row[nrows].j = j + jj;
    663:                   row[nrows].k = k + kk;
    664:                   row[nrows].c = m;
    665:                   for (l = 0; l < NPB; l++) vals[NPB * nrows + l] = ej[NPB * (3 * idx + m) + l];
    666:                   nrows++;
    667:                 }
    668:               }
    669:             }
    670:           }
    671:         }
    672:         PetscCall(MatSetValuesStencil(jac, nrows, row, NPB, col, vals, ADD_VALUES));
    673:       }
    674:     }
    675:   }

    677:   PetscCall(MatAssemblyBegin(jac, MAT_FLUSH_ASSEMBLY));
    678:   PetscCall(MatAssemblyEnd(jac, MAT_FLUSH_ASSEMBLY));

    680:   /* set the diagonal */
    681:   v[0] = 1.;
    682:   v[1] = 0.;
    683:   v[2] = 0.;
    684:   v[3] = 0.;
    685:   v[4] = 1.;
    686:   v[5] = 0.;
    687:   v[6] = 0.;
    688:   v[7] = 0.;
    689:   v[8] = 1.;
    690:   for (k = zs; k < zs + zm; k++) {
    691:     for (j = ys; j < ys + ym; j++) {
    692:       for (i = xs; i < xs + xm; i++) {
    693:         if (OnBoundary(i, j, k, mx, my, mz)) {
    694:           for (m = 0; m < 3; m++) {
    695:             col[m].i = i;
    696:             col[m].j = j;
    697:             col[m].k = k;
    698:             col[m].c = m;
    699:           }
    700:           PetscCall(MatSetValuesStencil(jac, 3, col, 3, col, v, INSERT_VALUES));
    701:         }
    702:       }
    703:     }
    704:   }

    706:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
    707:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));

    709:   PetscCall(DMDAVecRestoreArray(cda, C, &c));
    710:   PetscFunctionReturn(PETSC_SUCCESS);
    711: }

    713: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field ***x, Field ***f, void *ptr)
    714: {
    715:   /* values for each basis function at each quadrature point */
    716:   AppCtx  *user = (AppCtx *)ptr;
    717:   PetscInt i, j, k, l;
    718:   PetscInt ii, jj, kk;

    720:   Field      ef[NEB];
    721:   Field      ex[NEB];
    722:   CoordField ec[NEB];

    724:   PetscInt      xs = info->xs, ys = info->ys, zs = info->zs;
    725:   PetscInt      xm = info->xm, ym = info->ym, zm = info->zm;
    726:   PetscInt      xes, yes, zes, xee, yee, zee;
    727:   PetscInt      mx = info->mx, my = info->my, mz = info->mz;
    728:   DM            cda;
    729:   CoordField ***c;
    730:   Vec           C;

    732:   PetscFunctionBeginUser;
    733:   PetscCall(DMGetCoordinateDM(info->da, &cda));
    734:   PetscCall(DMGetCoordinatesLocal(info->da, &C));
    735:   PetscCall(DMDAVecGetArray(cda, C, &c));
    736:   PetscCall(DMDAGetInfo(info->da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    737:   PetscCall(DMDAGetCorners(info->da, &xs, &ys, &zs, &xm, &ym, &zm));

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

    787: PetscErrorCode TangentLoad(SNES snes, Vec X, Vec Q, void *ptr)
    788: {
    789:   /* values for each basis function at each quadrature point */
    790:   AppCtx  *user = (AppCtx *)ptr;
    791:   PetscInt xs, ys, zs;
    792:   PetscInt xm, ym, zm;
    793:   PetscInt mx, my, mz;
    794:   DM       da;
    795:   Vec      Xl, Ql;
    796:   Field ***x, ***q;
    797:   PetscInt i, j, k, l;
    798:   PetscInt ii, jj, kk;

    800:   Field      eq[NEB];
    801:   Field      ex[NEB];
    802:   CoordField ec[NEB];

    804:   PetscInt      xes, yes, zes, xee, yee, zee;
    805:   DM            cda;
    806:   CoordField ***c;
    807:   Vec           C;

    809:   PetscFunctionBeginUser;
    810:   /* update user context with current load parameter */
    811:   PetscCall(SNESNewtonALGetLoadParameter(snes, &user->load_factor));

    813:   PetscCall(SNESGetDM(snes, &da));
    814:   PetscCall(DMGetLocalVector(da, &Xl));
    815:   PetscCall(DMGetLocalVector(da, &Ql));
    816:   PetscCall(DMGlobalToLocal(da, X, INSERT_VALUES, Xl));

    818:   PetscCall(DMDAVecGetArray(da, Xl, &x));
    819:   PetscCall(DMDAVecGetArray(da, Ql, &q));

    821:   PetscCall(DMGetCoordinateDM(da, &cda));
    822:   PetscCall(DMGetCoordinatesLocal(da, &C));
    823:   PetscCall(DMDAVecGetArray(cda, C, &c));
    824:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    825:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));

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

    879: PetscErrorCode FormCoordinates(DM da, AppCtx *user)
    880: {
    881:   Vec           coords;
    882:   DM            cda;
    883:   PetscInt      mx, my, mz;
    884:   PetscInt      i, j, k, xs, ys, zs, xm, ym, zm;
    885:   CoordField ***x;

    887:   PetscFunctionBeginUser;
    888:   PetscCall(DMGetCoordinateDM(da, &cda));
    889:   PetscCall(DMCreateGlobalVector(cda, &coords));
    890:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    891:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
    892:   PetscCall(DMDAVecGetArray(da, coords, &x));
    893:   for (k = zs; k < zs + zm; k++) {
    894:     for (j = ys; j < ys + ym; j++) {
    895:       for (i = xs; i < xs + xm; i++) {
    896:         PetscReal cx  = ((PetscReal)i) / ((PetscReal)(mx - 1));
    897:         PetscReal cy  = ((PetscReal)j) / ((PetscReal)(my - 1));
    898:         PetscReal cz  = ((PetscReal)k) / ((PetscReal)(mz - 1));
    899:         PetscReal rad = user->rad + cy * user->height;
    900:         PetscReal ang = (cx - 0.5) * user->arc;
    901:         x[k][j][i][0] = rad * PetscSinReal(ang);
    902:         x[k][j][i][1] = rad * PetscCosReal(ang) - (user->rad + 0.5 * user->height) * PetscCosReal(-0.5 * user->arc);
    903:         x[k][j][i][2] = user->width * (cz - 0.5);
    904:       }
    905:     }
    906:   }
    907:   PetscCall(DMDAVecRestoreArray(da, coords, &x));
    908:   PetscCall(DMSetCoordinates(da, coords));
    909:   PetscCall(VecDestroy(&coords));
    910:   PetscFunctionReturn(PETSC_SUCCESS);
    911: }

    913: PetscErrorCode InitialGuess(DM da, AppCtx *user, Vec X)
    914: {
    915:   PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
    916:   PetscInt mx, my, mz;
    917:   Field ***x;

    919:   PetscFunctionBeginUser;
    920:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
    921:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    922:   PetscCall(DMDAVecGetArray(da, X, &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:         /* reference coordinates */
    928:         PetscReal p_x = ((PetscReal)i) / ((PetscReal)(mx - 1));
    929:         PetscReal p_y = ((PetscReal)j) / ((PetscReal)(my - 1));
    930:         PetscReal p_z = ((PetscReal)k) / ((PetscReal)(mz - 1));
    931:         PetscReal o_x = p_x;
    932:         PetscReal o_y = p_y;
    933:         PetscReal o_z = p_z;
    934:         x[k][j][i][0] = o_x - p_x;
    935:         x[k][j][i][1] = o_y - p_y;
    936:         x[k][j][i][2] = o_z - p_z;
    937:       }
    938:     }
    939:   }
    940:   PetscCall(DMDAVecRestoreArray(da, X, &x));
    941:   PetscFunctionReturn(PETSC_SUCCESS);
    942: }

    944: PetscErrorCode ArcLengthScaling(DM da, AppCtx *user, Vec V)
    945: {
    946:   PetscInt  i, j, k, xs, ys, zs, xm, ym, zm;
    947:   PetscInt  mx, my, mz;
    948:   Field  ***v;
    949:   PetscReal rad = user->rad + user->height;

    951:   PetscFunctionBeginUser;
    952:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
    953:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    954:   PetscCall(DMDAVecGetArray(da, V, &v));

    956:   for (k = zs; k < zs + zm; k++) {
    957:     for (j = ys; j < ys + ym; j++) {
    958:       for (i = xs; i < xs + xm; i++) {
    959:         v[k][j][i][0] = 2. / rad;
    960:         v[k][j][i][1] = 2. / rad;
    961:         v[k][j][i][2] = 2. / user->width;
    962:       }
    963:     }
    964:   }
    965:   PetscCall(DMDAVecRestoreArray(da, V, &v));
    966:   PetscFunctionReturn(PETSC_SUCCESS);
    967: }

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

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

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

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

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

    1030: /*TEST

    1032:    test:
    1033:       nsize: 2
    1034:       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
    1035:       requires: !single
    1036:       timeoutfactor: 3

    1038:    test:
    1039:       suffix: 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 -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
    1041:       requires: !single

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

    1048:    test:
    1049:       suffix: 4
    1050:       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
    1051:       requires: !single defined(PETSC_USE_INFO)
    1052:       filter: grep -h -e "KSP Residual norm" -e "SNES Function norm" -e "Number of SNES iterations" -e "mu:"

    1054:    test:
    1055:       suffix: 5
    1056:       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
    1057:       requires: !single

    1059:    test:
    1060:       suffix: 6
    1061:       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
    1062:       requires: !single

    1064:    test:
    1065:       nsize: {{1 2}}
    1066:       suffix: 7
    1067:       args: -dm_vec_type kokkos -dm_mat_type aijkokkos -da_refine 1 -rad 10.0 -young 10. -ploading -1. -loading -1. -rescale -pc_type mg -mg_levels_ksp_max_it 2 -snes_monitor_short -snes_type newtonal -snes_newtonal_step_size 10 -snes_newtonal_correction_type exact -ksp_rtol 1e-4
    1068:       requires: kokkos_kernels !complex !single
    1069:       timeoutfactor: 3

    1071:    test:
    1072:       nsize: {{1 2}}
    1073:       suffix: 8
    1074:       args: -dm_vec_type cuda -dm_mat_type aijcusparse -da_refine 1 -rad 10.0 -young 10. -ploading -1. -loading -1. -rescale -pc_type mg -mg_levels_ksp_max_it 2 -snes_monitor_short -snes_type newtonal -snes_newtonal_step_size 10 -snes_newtonal_correction_type exact -ksp_rtol 1e-4
    1075:       requires: cuda !complex !single
    1076:       timeoutfactor: 3

    1078:    test:
    1079:       nsize: {{1 2}}
    1080:       suffix: 9
    1081:       args: -da_refine 2 -rad 10.0 -young 10. -ploading -1. -loading -1. -rescale -pc_type mg -mg_levels_ksp_max_it 2 -snes_monitor_short -snes_type newtonal -snes_newtonal_step_size 10 -snes_newtonal_correction_type normal -ksp_rtol 1e-4
    1082:       requires: !single

    1084: TEST*/