Actual source code: ex17.c

      1: static char help[] = "Linear elasticity in 2d and 3d with finite elements.\n\
      2: We solve the elasticity problem in a rectangular\n\
      3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
      4: This example supports automatic convergence estimation\n\
      5: and eventually adaptivity.\n\n\n";

      7: /*
      8:   https://en.wikipedia.org/wiki/Linear_elasticity

     10:   Converting elastic constants:
     11:     lambda = E nu / ((1 + nu) (1 - 2 nu))
     12:     mu     = E / (2 (1 + nu))
     13: */

     15: #include <petscdmplex.h>
     16: #include <petscsnes.h>
     17: #include <petscds.h>
     18: #include <petscbag.h>
     19: #include <petscconvest.h>

     21: typedef enum {
     22:   SOL_VLAP_QUADRATIC,
     23:   SOL_ELAS_QUADRATIC,
     24:   SOL_VLAP_TRIG,
     25:   SOL_ELAS_TRIG,
     26:   SOL_ELAS_AXIAL_DISP,
     27:   SOL_ELAS_UNIFORM_STRAIN,
     28:   SOL_ELAS_GE,
     29:   SOL_MASS_QUADRATIC,
     30:   NUM_SOLUTION_TYPES
     31: } SolutionType;
     32: const char *solutionTypes[NUM_SOLUTION_TYPES + 1] = {"vlap_quad", "elas_quad", "vlap_trig", "elas_trig", "elas_axial_disp", "elas_uniform_strain", "elas_ge", "mass_quad", "unknown"};

     34: typedef enum {
     35:   DEFORM_NONE,
     36:   DEFORM_SHEAR,
     37:   DEFORM_STEP,
     38:   NUM_DEFORM_TYPES
     39: } DeformType;
     40: const char *deformTypes[NUM_DEFORM_TYPES + 1] = {"none", "shear", "step", "unknown"};

     42: typedef struct {
     43:   PetscScalar mu;     /* shear modulus */
     44:   PetscScalar lambda; /* Lame's first parameter */
     45:   PetscScalar N;      /* Tension force on right wall */
     46: } Parameter;

     48: typedef struct {
     49:   /* Domain and mesh definition */
     50:   char       dmType[256]; /* DM type for the solve */
     51:   DeformType deform;      /* Domain deformation type */
     52:   /* Problem definition */
     53:   SolutionType solType; /* Type of exact solution */
     54:   PetscBag     bag;     /* Problem parameters */
     55:   /* Solver definition */
     56:   PetscBool useNearNullspace; /* Use the rigid body modes as a near nullspace for AMG */
     57: } AppCtx;

     59: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
     60: {
     61:   PetscInt d;
     62:   for (d = 0; d < dim; ++d) u[d] = 0.0;
     63:   return PETSC_SUCCESS;
     64: }

     66: static PetscErrorCode ge_shift(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
     67: {
     68:   PetscInt d;
     69:   u[0] = 0.1;
     70:   for (d = 1; d < dim; ++d) u[d] = 0.0;
     71:   return PETSC_SUCCESS;
     72: }

     74: static PetscErrorCode quadratic_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
     75: {
     76:   u[0] = x[0] * x[0];
     77:   u[1] = x[1] * x[1] - 2.0 * x[0] * x[1];
     78:   return PETSC_SUCCESS;
     79: }

     81: static PetscErrorCode quadratic_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
     82: {
     83:   u[0] = x[0] * x[0];
     84:   u[1] = x[1] * x[1] - 2.0 * x[0] * x[1];
     85:   u[2] = x[2] * x[2] - 2.0 * x[1] * x[2];
     86:   return PETSC_SUCCESS;
     87: }

     89: /*
     90:   u = x^2
     91:   v = y^2 - 2xy
     92:   Delta <u,v> - f = <2, 2> - <2, 2>

     94:   u = x^2
     95:   v = y^2 - 2xy
     96:   w = z^2 - 2yz
     97:   Delta <u,v,w> - f = <2, 2, 2> - <2, 2, 2>
     98: */
     99: static void f0_vlap_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    100: {
    101:   PetscInt d;
    102:   for (d = 0; d < dim; ++d) f0[d] += 2.0;
    103: }

    105: /*
    106:   u = x^2
    107:   v = y^2 - 2xy
    108:   \varepsilon = / 2x     -y    \
    109:                 \ -y   2y - 2x /
    110:   Tr(\varepsilon) = div u = 2y
    111:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
    112:     = \lambda \partial_j (2y) + 2\mu < 2-1, 2 >
    113:     = \lambda < 0, 2 > + \mu < 2, 4 >

    115:   u = x^2
    116:   v = y^2 - 2xy
    117:   w = z^2 - 2yz
    118:   \varepsilon = / 2x     -y       0   \
    119:                 | -y   2y - 2x   -z   |
    120:                 \  0     -z    2z - 2y/
    121:   Tr(\varepsilon) = div u = 2z
    122:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
    123:     = \lambda \partial_j (2z) + 2\mu < 2-1, 2-1, 2 >
    124:     = \lambda < 0, 0, 2 > + \mu < 2, 2, 4 >
    125: */
    126: static void f0_elas_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    127: {
    128:   const PetscReal mu     = PetscRealPart(constants[0]);
    129:   const PetscReal lambda = PetscRealPart(constants[1]);

    131:   for (PetscInt d = 0; d < dim - 1; ++d) f0[d] += 2.0 * mu;
    132:   f0[dim - 1] += 2.0 * lambda + 4.0 * mu;
    133: }

    135: static void f0_mass_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    136: {
    137:   if (dim == 2) {
    138:     f0[0] -= x[0] * x[0];
    139:     f0[1] -= x[1] * x[1] - 2.0 * x[0] * x[1];
    140:   } else {
    141:     f0[0] -= x[0] * x[0];
    142:     f0[1] -= x[1] * x[1] - 2.0 * x[0] * x[1];
    143:     f0[2] -= x[2] * x[2] - 2.0 * x[1] * x[2];
    144:   }
    145: }

    147: static PetscErrorCode trig_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
    148: {
    149:   u[0] = PetscSinReal(2.0 * PETSC_PI * x[0]);
    150:   u[1] = PetscSinReal(2.0 * PETSC_PI * x[1]) - 2.0 * x[0] * x[1];
    151:   return PETSC_SUCCESS;
    152: }

    154: static PetscErrorCode trig_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
    155: {
    156:   u[0] = PetscSinReal(2.0 * PETSC_PI * x[0]);
    157:   u[1] = PetscSinReal(2.0 * PETSC_PI * x[1]) - 2.0 * x[0] * x[1];
    158:   u[2] = PetscSinReal(2.0 * PETSC_PI * x[2]) - 2.0 * x[1] * x[2];
    159:   return PETSC_SUCCESS;
    160: }

    162: /*
    163:   u = sin(2 pi x)
    164:   v = sin(2 pi y) - 2xy
    165:   Delta <u,v> - f = <-4 pi^2 u, -4 pi^2 v> - <-4 pi^2 sin(2 pi x), -4 pi^2 sin(2 pi y)>

    167:   u = sin(2 pi x)
    168:   v = sin(2 pi y) - 2xy
    169:   w = sin(2 pi z) - 2yz
    170:   Delta <u,v,2> - f = <-4 pi^2 u, -4 pi^2 v, -4 pi^2 w> - <-4 pi^2 sin(2 pi x), -4 pi^2 sin(2 pi y), -4 pi^2 sin(2 pi z)>
    171: */
    172: static void f0_vlap_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    173: {
    174:   PetscInt d;
    175:   for (d = 0; d < dim; ++d) f0[d] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
    176: }

    178: /*
    179:   u = sin(2 pi x)
    180:   v = sin(2 pi y) - 2xy
    181:   \varepsilon = / 2 pi cos(2 pi x)             -y        \
    182:                 \      -y          2 pi cos(2 pi y) - 2x /
    183:   Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
    184:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
    185:     = \lambda \partial_j 2 pi (cos(2 pi x) + cos(2 pi y)) + 2\mu < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) >
    186:     = \lambda < -4 pi^2 sin(2 pi x) - 2, -4 pi^2 sin(2 pi y) > + \mu < -8 pi^2 sin(2 pi x) - 2, -8 pi^2 sin(2 pi y) >

    188:   u = sin(2 pi x)
    189:   v = sin(2 pi y) - 2xy
    190:   w = sin(2 pi z) - 2yz
    191:   \varepsilon = / 2 pi cos(2 pi x)            -y                     0         \
    192:                 |         -y       2 pi cos(2 pi y) - 2x            -z         |
    193:                 \          0                  -z         2 pi cos(2 pi z) - 2y /
    194:   Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y
    195:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
    196:     = \lambda \partial_j (2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y) + 2\mu < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) - 1, -4 pi^2 sin(2 pi z) >
    197:     = \lambda < -4 pi^2 sin(2 pi x) - 2, -4 pi^2 sin(2 pi y) - 2, -4 pi^2 sin(2 pi z) > + 2\mu < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) - 1, -4 pi^2 sin(2 pi z) >
    198: */
    199: static void f0_elas_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    200: {
    201:   const PetscReal mu     = PetscRealPart(constants[0]);
    202:   const PetscReal lambda = PetscRealPart(constants[1]);
    203:   const PetscReal fact   = 4.0 * PetscSqr(PETSC_PI);

    205:   for (PetscInt d = 0; d < dim; ++d) f0[d] += -(2.0 * mu + lambda) * fact * PetscSinReal(2.0 * PETSC_PI * x[d]) - (d < dim - 1 ? 2.0 * (mu + lambda) : 0.0);
    206: }

    208: static PetscErrorCode axial_disp_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
    209: {
    210:   AppCtx    *user = (AppCtx *)ctx;
    211:   Parameter *param;

    213:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
    214:   {
    215:     const PetscReal mu     = PetscRealPart(param->mu);
    216:     const PetscReal lambda = PetscRealPart(param->lambda);
    217:     const PetscReal N      = PetscRealPart(param->N);

    219:     u[0] = (3. * lambda * lambda + 8. * lambda * mu + 4 * mu * mu) / (4 * mu * (3 * lambda * lambda + 5. * lambda * mu + 2 * mu * mu)) * N * x[0];
    220:     u[1] = 0.25 * lambda / mu / (lambda + mu) * N * x[1];
    221:     for (PetscInt d = 2; d < dim; ++d) u[d] = 0.0;
    222:   }
    223:   return PETSC_SUCCESS;
    224: }

    226: /*
    227:   We will pull/push on the right side of a block of linearly elastic material. The uniform traction conditions on the
    228:   right side of the box will result in a uniform strain along x and y. The Neumann BC is given by

    230:      n_i \sigma_{ij} = t_i

    232:   u = (1/(2\mu) - 1) x
    233:   v = -y
    234:   f = 0
    235:   t = <4\mu/\lambda (\lambda + \mu), 0>
    236:   \varepsilon = / 1/(2\mu) - 1   0 \
    237:                 \ 0             -1 /
    238:   Tr(\varepsilon) = div u = 1/(2\mu) - 2
    239:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
    240:     = \lambda \partial_j (1/(2\mu) - 2) + 2\mu < 0, 0 >
    241:     = \lambda < 0, 0 > + \mu < 0, 0 > = 0
    242:   NBC =  <1,0> . <4\mu/\lambda (\lambda + \mu), 0> = 4\mu/\lambda (\lambda + \mu)

    244:   u = x - 1/2
    245:   v = 0
    246:   w = 0
    247:   \varepsilon = / x  0  0 \
    248:                 | 0  0  0 |
    249:                 \ 0  0  0 /
    250:   Tr(\varepsilon) = div u = x
    251:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
    252:     = \lambda \partial_j x + 2\mu < 1, 0, 0 >
    253:     = \lambda < 1, 0, 0 > + \mu < 2, 0, 0 >
    254: */
    255: static void f0_elas_axial_disp_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    256: {
    257:   const PetscReal N = PetscRealPart(constants[2]);

    259:   f0[0] = N;
    260: }

    262: static PetscErrorCode uniform_strain_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
    263: {
    264:   const PetscReal eps_xx = 0.1;
    265:   const PetscReal eps_xy = 0.3;
    266:   const PetscReal eps_yy = 0.25;
    267:   PetscInt        d;

    269:   u[0] = eps_xx * x[0] + eps_xy * x[1];
    270:   u[1] = eps_xy * x[0] + eps_yy * x[1];
    271:   for (d = 2; d < dim; ++d) u[d] = 0.0;
    272:   return PETSC_SUCCESS;
    273: }

    275: static void f0_mass_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    276: {
    277:   const PetscInt Nc = dim;
    278:   PetscInt       c;

    280:   for (c = 0; c < Nc; ++c) f0[c] = u[c];
    281: }

    283: static void f1_vlap_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
    284: {
    285:   const PetscInt Nc = dim;
    286:   PetscInt       c, d;

    288:   for (c = 0; c < Nc; ++c)
    289:     for (d = 0; d < dim; ++d) f1[c * dim + d] += u_x[c * dim + d];
    290: }

    292: static void f1_elas_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
    293: {
    294:   const PetscReal mu     = PetscRealPart(constants[0]);
    295:   const PetscReal lambda = PetscRealPart(constants[1]);
    296:   const PetscInt  Nc     = dim;

    298:   for (PetscInt c = 0; c < Nc; ++c) {
    299:     for (PetscInt d = 0; d < dim; ++d) {
    300:       f1[c * dim + d] += mu * (u_x[c * dim + d] + u_x[d * dim + c]);
    301:       f1[c * dim + c] += lambda * u_x[d * dim + d];
    302:     }
    303:   }
    304: }

    306: static void g0_mass_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
    307: {
    308:   const PetscInt Nc = dim;
    309:   PetscInt       c;

    311:   for (c = 0; c < Nc; ++c) g0[c * Nc + c] = 1.0;
    312: }

    314: static void g3_vlap_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
    315: {
    316:   const PetscInt Nc = dim;
    317:   PetscInt       c, d;

    319:   for (c = 0; c < Nc; ++c) {
    320:     for (d = 0; d < dim; ++d) g3[((c * Nc + c) * dim + d) * dim + d] = 1.0;
    321:   }
    322: }

    324: /*
    325:   \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc

    327:   \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg}
    328:   = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc
    329: */
    330: static void g3_elas_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
    331: {
    332:   const PetscReal mu     = PetscRealPart(constants[0]);
    333:   const PetscReal lambda = PetscRealPart(constants[1]);
    334:   const PetscInt  Nc     = dim;

    336:   for (PetscInt c = 0; c < Nc; ++c) {
    337:     for (PetscInt d = 0; d < dim; ++d) {
    338:       g3[((c * Nc + c) * dim + d) * dim + d] += mu;
    339:       g3[((c * Nc + d) * dim + d) * dim + c] += mu;
    340:       g3[((c * Nc + d) * dim + c) * dim + d] += lambda;
    341:     }
    342:   }
    343: }

    345: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
    346: {
    347:   PetscInt sol = 0, def = 0;

    349:   PetscFunctionBeginUser;
    350:   options->deform           = DEFORM_NONE;
    351:   options->solType          = SOL_VLAP_QUADRATIC;
    352:   options->useNearNullspace = PETSC_TRUE;
    353:   PetscCall(PetscStrncpy(options->dmType, DMPLEX, 256));

    355:   PetscOptionsBegin(comm, "", "Linear Elasticity Problem Options", "DMPLEX");
    356:   PetscCall(PetscOptionsEList("-deform_type", "Type of domain deformation", "ex17.c", deformTypes, NUM_DEFORM_TYPES, deformTypes[options->deform], &def, NULL));
    357:   options->deform = (DeformType)def;
    358:   PetscCall(PetscOptionsEList("-sol_type", "Type of exact solution", "ex17.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL));
    359:   options->solType = (SolutionType)sol;
    360:   PetscCall(PetscOptionsBool("-near_nullspace", "Use the rigid body modes as an AMG near nullspace", "ex17.c", options->useNearNullspace, &options->useNearNullspace, NULL));
    361:   PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex17.c", DMList, options->dmType, options->dmType, 256, NULL));
    362:   PetscOptionsEnd();
    363:   PetscFunctionReturn(PETSC_SUCCESS);
    364: }

    366: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
    367: {
    368:   PetscBag   bag;
    369:   Parameter *p;

    371:   PetscFunctionBeginUser;
    372:   /* setup PETSc parameter bag */
    373:   PetscCall(PetscBagGetData(ctx->bag, (void **)&p));
    374:   PetscCall(PetscBagSetName(ctx->bag, "par", "Elastic Parameters"));
    375:   bag = ctx->bag;
    376:   PetscCall(PetscBagRegisterScalar(bag, &p->mu, 1.0, "mu", "Shear Modulus, Pa"));
    377:   PetscCall(PetscBagRegisterScalar(bag, &p->lambda, 1.0, "lambda", "Lame's first parameter, Pa"));
    378:   PetscCall(PetscBagRegisterScalar(bag, &p->N, -1.0, "N", "Tension on right wall, Pa"));
    379:   PetscCall(PetscBagSetFromOptions(bag));
    380:   {
    381:     PetscViewer       viewer;
    382:     PetscViewerFormat format;
    383:     PetscBool         flg;

    385:     PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg));
    386:     if (flg) {
    387:       PetscCall(PetscViewerPushFormat(viewer, format));
    388:       PetscCall(PetscBagView(bag, viewer));
    389:       PetscCall(PetscViewerFlush(viewer));
    390:       PetscCall(PetscViewerPopFormat(viewer));
    391:       PetscCall(PetscViewerDestroy(&viewer));
    392:     }
    393:   }
    394:   PetscFunctionReturn(PETSC_SUCCESS);
    395: }

    397: static PetscErrorCode DMPlexDistortGeometry(DM dm)
    398: {
    399:   DM           cdm;
    400:   DMLabel      label;
    401:   Vec          coordinates;
    402:   PetscScalar *coords;
    403:   PetscReal    mid = 0.5;
    404:   PetscInt     cdim, d, vStart, vEnd, v;

    406:   PetscFunctionBeginUser;
    407:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    408:   PetscCall(DMGetCoordinateDim(dm, &cdim));
    409:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    410:   PetscCall(DMGetLabel(dm, "marker", &label));
    411:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
    412:   PetscCall(VecGetArrayWrite(coordinates, &coords));
    413:   for (v = vStart; v < vEnd; ++v) {
    414:     PetscScalar *pcoords, shift;
    415:     PetscInt     val;

    417:     PetscCall(DMLabelGetValue(label, v, &val));
    418:     if (val >= 0) continue;
    419:     PetscCall(DMPlexPointLocalRef(cdm, v, coords, &pcoords));
    420:     shift = 0.2 * PetscAbsScalar(pcoords[0] - mid);
    421:     shift = PetscRealPart(pcoords[0]) > mid ? shift : -shift;
    422:     for (d = 1; d < cdim; ++d) pcoords[d] += shift;
    423:   }
    424:   PetscCall(VecRestoreArrayWrite(coordinates, &coords));
    425:   PetscFunctionReturn(PETSC_SUCCESS);
    426: }

    428: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
    429: {
    430:   PetscFunctionBeginUser;
    431:   PetscCall(DMCreate(comm, dm));
    432:   PetscCall(DMSetType(*dm, DMPLEX));
    433:   PetscCall(DMSetFromOptions(*dm));
    434:   switch (user->deform) {
    435:   case DEFORM_NONE:
    436:     break;
    437:   case DEFORM_SHEAR:
    438:     PetscCall(DMPlexShearGeometry(*dm, DM_X, NULL));
    439:     break;
    440:   case DEFORM_STEP:
    441:     PetscCall(DMPlexDistortGeometry(*dm));
    442:     break;
    443:   default:
    444:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid deformation type: %s (%d)", deformTypes[PetscMin(user->deform, NUM_DEFORM_TYPES)], user->deform);
    445:   }
    446:   PetscCall(DMSetApplicationContext(*dm, user));
    447:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
    448:   PetscFunctionReturn(PETSC_SUCCESS);
    449: }

    451: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
    452: {
    453:   PetscErrorCode (*exact)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
    454:   Parameter    *param;
    455:   PetscDS       ds;
    456:   PetscWeakForm wf;
    457:   DMLabel       label;
    458:   PetscInt      id, bd;
    459:   PetscInt      dim;

    461:   PetscFunctionBeginUser;
    462:   PetscCall(DMGetDS(dm, &ds));
    463:   PetscCall(PetscDSGetWeakForm(ds, &wf));
    464:   PetscCall(PetscDSGetSpatialDimension(ds, &dim));
    465:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
    466:   switch (user->solType) {
    467:   case SOL_MASS_QUADRATIC:
    468:     PetscCall(PetscDSSetResidual(ds, 0, f0_mass_u, NULL));
    469:     PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_mass_uu, NULL, NULL, NULL));
    470:     PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 1, f0_mass_quadratic_u, 0, NULL));
    471:     switch (dim) {
    472:     case 2:
    473:       exact = quadratic_2d_u;
    474:       break;
    475:     case 3:
    476:       exact = quadratic_3d_u;
    477:       break;
    478:     default:
    479:       SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %" PetscInt_FMT, dim);
    480:     }
    481:     break;
    482:   case SOL_VLAP_QUADRATIC:
    483:     PetscCall(PetscDSSetResidual(ds, 0, f0_vlap_quadratic_u, f1_vlap_u));
    484:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_vlap_uu));
    485:     switch (dim) {
    486:     case 2:
    487:       exact = quadratic_2d_u;
    488:       break;
    489:     case 3:
    490:       exact = quadratic_3d_u;
    491:       break;
    492:     default:
    493:       SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %" PetscInt_FMT, dim);
    494:     }
    495:     break;
    496:   case SOL_ELAS_QUADRATIC:
    497:     PetscCall(PetscDSSetResidual(ds, 0, f0_elas_quadratic_u, f1_elas_u));
    498:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu));
    499:     switch (dim) {
    500:     case 2:
    501:       exact = quadratic_2d_u;
    502:       break;
    503:     case 3:
    504:       exact = quadratic_3d_u;
    505:       break;
    506:     default:
    507:       SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %" PetscInt_FMT, dim);
    508:     }
    509:     break;
    510:   case SOL_VLAP_TRIG:
    511:     PetscCall(PetscDSSetResidual(ds, 0, f0_vlap_trig_u, f1_vlap_u));
    512:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_vlap_uu));
    513:     switch (dim) {
    514:     case 2:
    515:       exact = trig_2d_u;
    516:       break;
    517:     case 3:
    518:       exact = trig_3d_u;
    519:       break;
    520:     default:
    521:       SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %" PetscInt_FMT, dim);
    522:     }
    523:     break;
    524:   case SOL_ELAS_TRIG:
    525:     PetscCall(PetscDSSetResidual(ds, 0, f0_elas_trig_u, f1_elas_u));
    526:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu));
    527:     switch (dim) {
    528:     case 2:
    529:       exact = trig_2d_u;
    530:       break;
    531:     case 3:
    532:       exact = trig_3d_u;
    533:       break;
    534:     default:
    535:       SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %" PetscInt_FMT, dim);
    536:     }
    537:     break;
    538:   case SOL_ELAS_AXIAL_DISP:
    539:     PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_elas_u));
    540:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu));
    541:     id = dim == 3 ? 5 : 2;
    542:     PetscCall(DMGetLabel(dm, "marker", &label));
    543:     PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right", label, 1, &id, 0, 0, NULL, (void (*)(void))NULL, NULL, user, &bd));
    544:     PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    545:     PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_elas_axial_disp_bd_u, 0, NULL));
    546:     exact = axial_disp_u;
    547:     break;
    548:   case SOL_ELAS_UNIFORM_STRAIN:
    549:     PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_elas_u));
    550:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu));
    551:     exact = uniform_strain_u;
    552:     break;
    553:   case SOL_ELAS_GE:
    554:     PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_elas_u));
    555:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu));
    556:     exact = zero; /* No exact solution available */
    557:     break;
    558:   default:
    559:     SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType);
    560:   }
    561:   PetscCall(PetscDSSetExactSolution(ds, 0, exact, user));
    562:   PetscCall(DMGetLabel(dm, "marker", &label));
    563:   if (user->solType == SOL_ELAS_AXIAL_DISP) {
    564:     PetscInt cmp;

    566:     id  = dim == 3 ? 6 : 4;
    567:     cmp = 0;
    568:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left", label, 1, &id, 0, 1, &cmp, (void (*)(void))zero, NULL, user, NULL));
    569:     cmp = dim == 3 ? 2 : 1;
    570:     id  = dim == 3 ? 1 : 1;
    571:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom", label, 1, &id, 0, 1, &cmp, (void (*)(void))zero, NULL, user, NULL));
    572:     if (dim == 3) {
    573:       cmp = 1;
    574:       id  = 3;
    575:       PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "front", label, 1, &id, 0, 1, &cmp, (void (*)(void))zero, NULL, user, NULL));
    576:     }
    577:   } else if (user->solType == SOL_ELAS_GE) {
    578:     PetscInt cmp;

    580:     id = dim == 3 ? 6 : 4;
    581:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left", label, 1, &id, 0, 0, NULL, (void (*)(void))zero, NULL, user, NULL));
    582:     id  = dim == 3 ? 5 : 2;
    583:     cmp = 0;
    584:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right", label, 1, &id, 0, 1, &cmp, (void (*)(void))ge_shift, NULL, user, NULL));
    585:   } else {
    586:     id = 1;
    587:     PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))exact, NULL, user, NULL));
    588:   }
    589:   /* Setup constants */
    590:   {
    591:     PetscScalar constants[3];

    593:     constants[0] = param->mu;     /* shear modulus, Pa */
    594:     constants[1] = param->lambda; /* Lame's first parameter, Pa */
    595:     constants[2] = param->N;      /* Tension on right wall, Pa */
    596:     PetscCall(PetscDSSetConstants(ds, 3, constants));
    597:   }
    598:   PetscFunctionReturn(PETSC_SUCCESS);
    599: }

    601: static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
    602: {
    603:   PetscFunctionBegin;
    604:   PetscCall(DMPlexCreateRigidBody(dm, origField, nullspace));
    605:   PetscFunctionReturn(PETSC_SUCCESS);
    606: }

    608: PetscErrorCode SetupFE(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx)
    609: {
    610:   AppCtx        *user = (AppCtx *)ctx;
    611:   DM             cdm  = dm;
    612:   PetscFE        fe;
    613:   char           prefix[PETSC_MAX_PATH_LEN];
    614:   DMPolytopeType ct;
    615:   PetscInt       dim, cStart;

    617:   PetscFunctionBegin;
    618:   /* Create finite element */
    619:   PetscCall(DMGetDimension(dm, &dim));
    620:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
    621:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
    622:   PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
    623:   PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, dim, ct, name ? prefix : NULL, -1, &fe));
    624:   PetscCall(PetscObjectSetName((PetscObject)fe, name));
    625:   /* Set discretization and boundary conditions for each mesh */
    626:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
    627:   PetscCall(DMCreateDS(dm));
    628:   PetscCall((*setup)(dm, user));
    629:   while (cdm) {
    630:     PetscCall(DMCopyDisc(dm, cdm));
    631:     if (user->useNearNullspace) PetscCall(DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace));
    632:     PetscCall(DMGetCoarseDM(cdm, &cdm));
    633:   }
    634:   PetscCall(PetscFEDestroy(&fe));
    635:   PetscFunctionReturn(PETSC_SUCCESS);
    636: }

    638: int main(int argc, char **argv)
    639: {
    640:   DM     dm;   /* Problem specification */
    641:   SNES   snes; /* Nonlinear solver */
    642:   Vec    u;    /* Solutions */
    643:   AppCtx user; /* User-defined work context */

    645:   PetscFunctionBeginUser;
    646:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    647:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
    648:   PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &user.bag));
    649:   PetscCall(SetupParameters(PETSC_COMM_WORLD, &user));
    650:   /* Primal system */
    651:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
    652:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
    653:   PetscCall(SNESSetDM(snes, dm));
    654:   PetscCall(SetupFE(dm, "displacement", SetupPrimalProblem, &user));
    655:   PetscCall(DMCreateGlobalVector(dm, &u));
    656:   PetscCall(VecSet(u, 0.0));
    657:   PetscCall(PetscObjectSetName((PetscObject)u, "displacement"));
    658:   PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
    659:   PetscCall(SNESSetFromOptions(snes));
    660:   PetscCall(DMSNESCheckFromOptions(snes, u));
    661:   PetscCall(SNESSolve(snes, NULL, u));
    662:   PetscCall(SNESGetSolution(snes, &u));
    663:   PetscCall(VecViewFromOptions(u, NULL, "-displacement_view"));
    664:   /* Cleanup */
    665:   PetscCall(VecDestroy(&u));
    666:   PetscCall(SNESDestroy(&snes));
    667:   PetscCall(DMDestroy(&dm));
    668:   PetscCall(PetscBagDestroy(&user.bag));
    669:   PetscCall(PetscFinalize());
    670:   return 0;
    671: }

    673: /*TEST

    675:   testset:
    676:     args: -dm_plex_box_faces 1,1,1

    678:     test:
    679:       suffix: 2d_p1_quad_vlap
    680:       requires: triangle
    681:       args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
    682:     test:
    683:       suffix: 2d_p2_quad_vlap
    684:       requires: triangle
    685:       args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
    686:     test:
    687:       suffix: 2d_p3_quad_vlap
    688:       requires: triangle
    689:       args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
    690:     test:
    691:       suffix: 2d_q1_quad_vlap
    692:       args: -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
    693:     test:
    694:       suffix: 2d_q2_quad_vlap
    695:       args: -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
    696:     test:
    697:       suffix: 2d_q3_quad_vlap
    698:       requires: !single
    699:       args: -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
    700:     test:
    701:       suffix: 2d_p1_quad_elas
    702:       requires: triangle
    703:       args: -sol_type elas_quad -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
    704:     test:
    705:       suffix: 2d_p2_quad_elas
    706:       requires: triangle
    707:       args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001
    708:     test:
    709:       suffix: 2d_p3_quad_elas
    710:       requires: triangle
    711:       args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001
    712:     test:
    713:       suffix: 2d_q1_quad_elas
    714:       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    715:     test:
    716:       suffix: 2d_q1_quad_elas_shear
    717:       args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    718:     test:
    719:       suffix: 2d_q2_quad_elas
    720:       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001
    721:     test:
    722:       suffix: 2d_q2_quad_elas_shear
    723:       args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 2 -dmsnes_check
    724:     test:
    725:       suffix: 2d_q3_quad_elas
    726:       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001
    727:     test:
    728:       suffix: 2d_q3_quad_elas_shear
    729:       requires: !single
    730:       args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 3 -dmsnes_check

    732:     test:
    733:       suffix: 3d_p1_quad_vlap
    734:       requires: ctetgen
    735:       args: -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    736:     test:
    737:       suffix: 3d_p2_quad_vlap
    738:       requires: ctetgen
    739:       args: -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
    740:     test:
    741:       suffix: 3d_p3_quad_vlap
    742:       requires: ctetgen
    743:       args: -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
    744:     test:
    745:       suffix: 3d_q1_quad_vlap
    746:       args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    747:     test:
    748:       suffix: 3d_q2_quad_vlap
    749:       args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
    750:     test:
    751:       suffix: 3d_q3_quad_vlap
    752:       args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
    753:     test:
    754:       suffix: 3d_p1_quad_elas
    755:       requires: ctetgen
    756:       args: -sol_type elas_quad -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    757:     test:
    758:       suffix: 3d_p2_quad_elas
    759:       requires: ctetgen
    760:       args: -sol_type elas_quad -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
    761:     test:
    762:       suffix: 3d_p3_quad_elas
    763:       requires: ctetgen
    764:       args: -sol_type elas_quad -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
    765:     test:
    766:       suffix: 3d_q1_quad_elas
    767:       args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    768:     test:
    769:       suffix: 3d_q2_quad_elas
    770:       args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
    771:     test:
    772:       suffix: 3d_q3_quad_elas
    773:       requires: !single
    774:       args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001

    776:     test:
    777:       suffix: 2d_p1_trig_vlap
    778:       requires: triangle
    779:       args: -sol_type vlap_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    780:     test:
    781:       suffix: 2d_p2_trig_vlap
    782:       requires: triangle
    783:       args: -sol_type vlap_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    784:     test:
    785:       suffix: 2d_p3_trig_vlap
    786:       requires: triangle
    787:       args: -sol_type vlap_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    788:     test:
    789:       suffix: 2d_q1_trig_vlap
    790:       args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    791:     test:
    792:       suffix: 2d_q2_trig_vlap
    793:       args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    794:     test:
    795:       suffix: 2d_q3_trig_vlap
    796:       args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    797:     test:
    798:       suffix: 2d_p1_trig_elas
    799:       requires: triangle
    800:       args: -sol_type elas_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    801:     test:
    802:       suffix: 2d_p2_trig_elas
    803:       requires: triangle
    804:       args: -sol_type elas_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    805:     test:
    806:       suffix: 2d_p3_trig_elas
    807:       requires: triangle
    808:       args: -sol_type elas_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    809:     test:
    810:       suffix: 2d_q1_trig_elas
    811:       args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    812:     test:
    813:       suffix: 2d_q1_trig_elas_shear
    814:       args: -sol_type elas_trig -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    815:     test:
    816:       suffix: 2d_q2_trig_elas
    817:       args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    818:     test:
    819:       suffix: 2d_q2_trig_elas_shear
    820:       args: -sol_type elas_trig -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    821:     test:
    822:       suffix: 2d_q3_trig_elas
    823:       args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    824:     test:
    825:       suffix: 2d_q3_trig_elas_shear
    826:       args: -sol_type elas_trig -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate

    828:     test:
    829:       suffix: 3d_p1_trig_vlap
    830:       requires: ctetgen
    831:       args: -sol_type vlap_trig -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    832:     test:
    833:       suffix: 3d_p2_trig_vlap
    834:       requires: ctetgen
    835:       args: -sol_type vlap_trig -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    836:     test:
    837:       suffix: 3d_p3_trig_vlap
    838:       requires: ctetgen
    839:       args: -sol_type vlap_trig -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    840:     test:
    841:       suffix: 3d_q1_trig_vlap
    842:       args: -sol_type vlap_trig -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    843:     test:
    844:       suffix: 3d_q2_trig_vlap
    845:       args: -sol_type vlap_trig -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    846:     test:
    847:       suffix: 3d_q3_trig_vlap
    848:       requires: !__float128
    849:       args: -sol_type vlap_trig -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    850:     test:
    851:       suffix: 3d_p1_trig_elas
    852:       requires: ctetgen
    853:       args: -sol_type elas_trig -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    854:     test:
    855:       suffix: 3d_p2_trig_elas
    856:       requires: ctetgen
    857:       args: -sol_type elas_trig -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    858:     test:
    859:       suffix: 3d_p3_trig_elas
    860:       requires: ctetgen
    861:       args: -sol_type elas_trig -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    862:     test:
    863:       suffix: 3d_q1_trig_elas
    864:       args: -sol_type elas_trig -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    865:     test:
    866:       suffix: 3d_q2_trig_elas
    867:       args: -sol_type elas_trig -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    868:     test:
    869:       suffix: 3d_q3_trig_elas
    870:       requires: !__float128
    871:       args: -sol_type elas_trig -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate

    873:     test:
    874:       suffix: 2d_p1_axial_elas
    875:       requires: triangle
    876:       args: -sol_type elas_axial_disp -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    877:     test:
    878:       suffix: 2d_p2_axial_elas
    879:       requires: triangle
    880:       args: -sol_type elas_axial_disp -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
    881:     test:
    882:       suffix: 2d_p3_axial_elas
    883:       requires: triangle
    884:       args: -sol_type elas_axial_disp -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
    885:     test:
    886:       suffix: 2d_q1_axial_elas
    887:       args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 1 -dmsnes_check .0001 -pc_type lu
    888:     test:
    889:       suffix: 2d_q2_axial_elas
    890:       args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
    891:     test:
    892:       suffix: 2d_q3_axial_elas
    893:       args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu

    895:     test:
    896:       suffix: 2d_p1_uniform_elas
    897:       requires: triangle
    898:       args: -sol_type elas_uniform_strain -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    899:     test:
    900:       suffix: 2d_p2_uniform_elas
    901:       requires: triangle
    902:       args: -sol_type elas_uniform_strain -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    903:     test:
    904:       suffix: 2d_p3_uniform_elas
    905:       requires: triangle
    906:       args: -sol_type elas_uniform_strain -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    907:     test:
    908:       suffix: 2d_q1_uniform_elas
    909:       args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    910:     test:
    911:       suffix: 2d_q2_uniform_elas
    912:       requires: !single
    913:       args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    914:     test:
    915:       suffix: 2d_q3_uniform_elas
    916:       requires: !single
    917:       args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    918:     test:
    919:       suffix: 2d_p1_uniform_elas_step
    920:       requires: triangle
    921:       args: -sol_type elas_uniform_strain -deform_type step -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu

    923:   testset:
    924:     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -deform_type step -displacement_petscspace_degree 1 -dmsnes_check .0001 -pc_type lu

    926:     test:
    927:       suffix: 2d_q1_uniform_elas_step
    928:       args: -sol_type elas_uniform_strain -dm_refine 2
    929:     test:
    930:       suffix: 2d_q1_quad_vlap_step
    931:       args:
    932:     test:
    933:       suffix: 2d_q2_quad_vlap_step
    934:       args: -displacement_petscspace_degree 2
    935:     test:
    936:       suffix: 2d_q1_quad_mass_step
    937:       args: -sol_type mass_quad

    939:   testset:
    940:     filter: grep -v "variant HERMITIAN"
    941:     args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_lower -5,-5,-0.25 -dm_plex_box_upper 5,5,0.25 \
    942:           -dm_plex_box_faces 5,5,2 -dm_plex_separate_marker -dm_refine 0 -petscpartitioner_type simple \
    943:           -sol_type elas_ge

    945:     test:
    946:       suffix: ge_q1_0
    947:       args: -displacement_petscspace_degree 1 \
    948:             -snes_max_it 2 -snes_rtol 1.e-10 \
    949:             -ksp_type cg -ksp_rtol 1.e-10 -ksp_max_it 100 -ksp_norm_type unpreconditioned \
    950:             -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \
    951:               -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true \
    952:               -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 \
    953:               -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.1 -mg_levels_pc_type jacobi \
    954:               -matptap_via scalable
    955:     test:
    956:       suffix: ge_q1_gmg
    957:       args: -displacement_petscspace_degree 1 \
    958:             -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
    959:             -snes_max_it 2 -snes_rtol 1.e-10 \
    960:             -ksp_type cg -ksp_rtol 1.e-10 -ksp_max_it 100 -ksp_norm_type unpreconditioned \
    961:             -pc_type mg -pc_mg_type full \
    962:               -mg_levels_ksp_max_it 4 -mg_levels_esteig_ksp_type cg \
    963:               -mg_levels_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \
    964:               -mg_levels_pc_type jacobi
    965:     test:
    966:       nsize: 5
    967:       suffix: ge_q1_gdsw
    968:       args: -snes_max_it 1 -ksp_type cg -ksp_norm_type natural -displacement_petscspace_degree 1 -snes_monitor_short -ksp_monitor_short -pc_type mg -pc_mg_adapt_interp_coarse_space gdsw -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type bjacobi -mg_levels_esteig_ksp_type cg -mg_levels_sub_pc_type icc -mg_coarse_redundant_pc_type cholesky -ksp_view

    970: TEST*/