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
      9: */

     11: #include <petscdmplex.h>
     12: #include <petscsnes.h>
     13: #include <petscds.h>
     14: #include <petscconvest.h>

     16: typedef enum {SOL_VLAP_QUADRATIC, SOL_ELAS_QUADRATIC, SOL_VLAP_TRIG, SOL_ELAS_TRIG, SOL_ELAS_AXIAL_DISP, SOL_ELAS_UNIFORM_STRAIN, NUM_SOLUTION_TYPES} SolutionType;
     17: const char *solutionTypes[NUM_SOLUTION_TYPES+1] = {"vlap_quad", "elas_quad", "vlap_trig", "elas_trig", "elas_axial_disp", "elas_uniform_strain", "unknown"};

     19: typedef struct {
     20:   /* Domain and mesh definition */
     21:   char         dmType[256]; /* DM type for the solve */
     22:   PetscInt     dim;         /* The topological mesh dimension */
     23:   PetscBool    simplex;     /* Simplicial mesh */
     24:   PetscInt     cells[3];    /* The initial domain division */
     25:   PetscBool    shear;       /* Shear the domain */
     26:   /* Problem definition */
     27:   SolutionType solType;     /* Type of exact solution */
     28:   /* Solver definition */
     29:   PetscBool    useNearNullspace; /* Use the rigid body modes as a near nullspace for AMG */
     30: } AppCtx;

     32: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
     33: {
     34:   PetscInt d;
     35:   for (d = 0; d < dim; ++d) u[d] = 0.0;
     36:   return 0;
     37: }

     39: static PetscErrorCode quadratic_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
     40: {
     41:   u[0] = x[0]*x[0];
     42:   u[1] = x[1]*x[1] - 2.0*x[0]*x[1];
     43:   return 0;
     44: }

     46: static PetscErrorCode quadratic_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
     47: {
     48:   u[0] = x[0]*x[0];
     49:   u[1] = x[1]*x[1] - 2.0*x[0]*x[1];
     50:   u[2] = x[2]*x[2] - 2.0*x[1]*x[2];
     51:   return 0;
     52: }

     54: /*
     55:   u = x^2
     56:   v = y^2 - 2xy
     57:   Delta <u,v> - f = <2, 2> - <2, 2>

     59:   u = x^2
     60:   v = y^2 - 2xy
     61:   w = z^2 - 2yz
     62:   Delta <u,v,w> - f = <2, 2, 2> - <2, 2, 2>
     63: */
     64: static void f0_vlap_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
     65:                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
     66:                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
     67:                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
     68: {
     69:   PetscInt d;
     70:   for (d = 0; d < dim; ++d) f0[d] += 2.0;
     71: }

     73: /*
     74:   u = x^2
     75:   v = y^2 - 2xy
     76:   \varepsilon = / 2x     -y    \
     77:                 \ -y   2y - 2x /
     78:   Tr(\varepsilon) = div u = 2y
     79:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
     80:     = \lambda \partial_j (2y) + 2\mu < 2-1, 2 >
     81:     = \lambda < 0, 2 > + \mu < 2, 4 >

     83:   u = x^2
     84:   v = y^2 - 2xy
     85:   w = z^2 - 2yz
     86:   \varepsilon = / 2x     -y       0   \
     87:                 | -y   2y - 2x   -z   |
     88:                 \  0     -z    2z - 2y/
     89:   Tr(\varepsilon) = div u = 2z
     90:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
     91:     = \lambda \partial_j (2z) + 2\mu < 2-1, 2-1, 2 >
     92:     = \lambda < 0, 0, 2 > + \mu < 2, 2, 4 >
     93: */
     94: static void f0_elas_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
     95:                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
     96:                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
     97:                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
     98: {
     99:   const PetscReal mu     = 1.0;
    100:   const PetscReal lambda = 1.0;
    101:   PetscInt        d;

    103:   for (d = 0; d < dim-1; ++d) f0[d] += 2.0*mu;
    104:   f0[dim-1] += 2.0*lambda + 4.0*mu;
    105: }

    107: static PetscErrorCode trig_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
    108: {
    109:   u[0] = PetscSinReal(2.0*PETSC_PI*x[0]);
    110:   u[1] = PetscSinReal(2.0*PETSC_PI*x[1]) - 2.0*x[0]*x[1];
    111:   return 0;
    112: }

    114: static PetscErrorCode trig_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
    115: {
    116:   u[0] = PetscSinReal(2.0*PETSC_PI*x[0]);
    117:   u[1] = PetscSinReal(2.0*PETSC_PI*x[1]) - 2.0*x[0]*x[1];
    118:   u[2] = PetscSinReal(2.0*PETSC_PI*x[2]) - 2.0*x[1]*x[2];
    119:   return 0;
    120: }

    122: /*
    123:   u = sin(2 pi x)
    124:   v = sin(2 pi y) - 2xy
    125:   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)>

    127:   u = sin(2 pi x)
    128:   v = sin(2 pi y) - 2xy
    129:   w = sin(2 pi z) - 2yz
    130:   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)>
    131: */
    132: static void f0_vlap_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    133:                            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    134:                            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    135:                            PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    136: {
    137:   PetscInt d;
    138:   for (d = 0; d < dim; ++d) f0[d] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]);
    139: }

    141: /*
    142:   u = sin(2 pi x)
    143:   v = sin(2 pi y) - 2xy
    144:   \varepsilon = / 2 pi cos(2 pi x)             -y        \
    145:                 \      -y          2 pi cos(2 pi y) - 2x /
    146:   Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
    147:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
    148:     = \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) >
    149:     = \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) >

    151:   u = sin(2 pi x)
    152:   v = sin(2 pi y) - 2xy
    153:   w = sin(2 pi z) - 2yz
    154:   \varepsilon = / 2 pi cos(2 pi x)            -y                     0         \
    155:                 |         -y       2 pi cos(2 pi y) - 2x            -z         |
    156:                 \          0                  -z         2 pi cos(2 pi z) - 2y /
    157:   Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y
    158:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
    159:     = \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) >
    160:     = \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) >
    161: */
    162: static void f0_elas_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    163:                            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    164:                            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    165:                            PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    166: {
    167:   const PetscReal mu     = 1.0;
    168:   const PetscReal lambda = 1.0;
    169:   const PetscReal fact   = 4.0*PetscSqr(PETSC_PI);
    170:   PetscInt        d;

    172:   for (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);
    173: }

    175: static PetscErrorCode axial_disp_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
    176: {
    177:   const PetscReal mu     = 1.0;
    178:   const PetscReal lambda = 1.0;
    179:   const PetscReal N      = 1.0;
    180:   PetscInt d;

    182:   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];
    183:   u[1] = -0.25*lambda/mu/(lambda+mu)*N*x[1];
    184:   for (d = 2; d < dim; ++d) u[d] = 0.0;
    185:   return 0;
    186: }

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

    192:      n_i \sigma_{ij} = t_i

    194:   u = (1/(2\mu) - 1) x
    195:   v = -y
    196:   f = 0
    197:   t = <4\mu/\lambda (\lambda + \mu), 0>
    198:   \varepsilon = / 1/(2\mu) - 1   0 \
    199:                 \ 0             -1 /
    200:   Tr(\varepsilon) = div u = 1/(2\mu) - 2
    201:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
    202:     = \lambda \partial_j (1/(2\mu) - 2) + 2\mu < 0, 0 >
    203:     = \lambda < 0, 0 > + \mu < 0, 0 > = 0
    204:   NBC =  <1,0> . <4\mu/\lambda (\lambda + \mu), 0> = 4\mu/\lambda (\lambda + \mu)

    206:   u = x - 1/2
    207:   v = 0
    208:   w = 0
    209:   \varepsilon = / x  0  0 \
    210:                 | 0  0  0 |
    211:                 \ 0  0  0 /
    212:   Tr(\varepsilon) = div u = x
    213:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
    214:     = \lambda \partial_j x + 2\mu < 1, 0, 0 >
    215:     = \lambda < 1, 0, 0 > + \mu < 2, 0, 0 >
    216: */
    217: static void f0_elas_axial_disp_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    218:                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    219:                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    220:                                     PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
    221: {
    222:   const PetscReal N = -1.0;

    224:   f0[0] = N;
    225: }

    227: static PetscErrorCode uniform_strain_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
    228: {
    229:   const PetscReal eps_xx = 0.1;
    230:   const PetscReal eps_xy = 0.3;
    231:   const PetscReal eps_yy = 0.25;
    232:   PetscInt d;

    234:   u[0] = eps_xx*x[0] + eps_xy*x[1];
    235:   u[1] = eps_xy*x[0] + eps_yy*x[1];
    236:   for (d = 2; d < dim; ++d) u[d] = 0.0;
    237:   return 0;
    238: }

    240: static void f1_vlap_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    241:                       const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    242:                       const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    243:                       PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
    244: {
    245:   const PetscInt Nc = dim;
    246:   PetscInt       c, d;

    248:   for (c = 0; c < Nc; ++c) for (d = 0; d < dim; ++d) f1[c*dim+d] += u_x[c*dim+d];
    249: }

    251: static void f1_elas_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    252:                       const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    253:                       const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    254:                       PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
    255: {
    256:   const PetscInt  Nc     = dim;
    257:   const PetscReal mu     = 1.0;
    258:   const PetscReal lambda = 1.0;
    259:   PetscInt        c, d;

    261:   for (c = 0; c < Nc; ++c) {
    262:     for (d = 0; d < dim; ++d) {
    263:       f1[c*dim+d] += mu*(u_x[c*dim+d] + u_x[d*dim+c]);
    264:       f1[c*dim+c] += lambda*u_x[d*dim+d];
    265:     }
    266:   }
    267: }

    269: static void g3_vlap_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    270:                        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    271:                        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    272:                        PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
    273: {
    274:   const PetscInt Nc = dim;
    275:   PetscInt       c, d;

    277:   for (c = 0; c < Nc; ++c) {
    278:     for (d = 0; d < dim; ++d) {
    279:       g3[((c*Nc + c)*dim + d)*dim + d] = 1.0;
    280:     }
    281:   }
    282: }

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

    287:   \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg}
    288:   = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc
    289: */
    290: static void g3_elas_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    291:                        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    292:                        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    293:                        PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
    294: {
    295:   const PetscInt  Nc     = dim;
    296:   const PetscReal mu     = 1.0;
    297:   const PetscReal lambda = 1.0;
    298:   PetscInt        c, d;

    300:   for (c = 0; c < Nc; ++c) {
    301:     for (d = 0; d < dim; ++d) {
    302:       g3[((c*Nc + c)*dim + d)*dim + d] += mu;
    303:       g3[((c*Nc + d)*dim + d)*dim + c] += mu;
    304:       g3[((c*Nc + d)*dim + c)*dim + d] += lambda;
    305:     }
    306:   }
    307: }

    309: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
    310: {
    311:   PetscInt       n = 3, sol;

    315:   options->dim              = 2;
    316:   options->cells[0]         = 1;
    317:   options->cells[1]         = 1;
    318:   options->cells[2]         = 1;
    319:   options->simplex          = PETSC_TRUE;
    320:   options->shear            = PETSC_FALSE;
    321:   options->solType          = SOL_VLAP_QUADRATIC;
    322:   options->useNearNullspace = PETSC_TRUE;
    323:   PetscStrncpy(options->dmType, DMPLEX, 256);

    325:   PetscOptionsBegin(comm, "", "Linear Elasticity Problem Options", "DMPLEX");
    326:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex17.c", options->dim, &options->dim, NULL);
    327:   PetscOptionsIntArray("-cells", "The initial mesh division", "ex17.c", options->cells, &n, NULL);
    328:   PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex17.c", options->simplex, &options->simplex, NULL);
    329:   PetscOptionsBool("-shear", "Shear the domain", "ex17.c", options->shear, &options->shear, NULL);
    330:   sol  = options->solType;
    331:   PetscOptionsEList("-sol_type", "Type of exact solution", "ex17.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL);
    332:   options->solType = (SolutionType) sol;
    333:   PetscOptionsBool("-near_nullspace", "Use the rigid body modes as an AMG near nullspace", "ex17.c", options->useNearNullspace, &options->useNearNullspace, NULL);
    334:   PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex17.c", DMList, options->dmType, options->dmType, 256, NULL);
    335:   PetscOptionsEnd();
    336:   return(0);
    337: }

    339: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
    340: {
    341:   PetscBool      flg;

    345:   DMPlexCreateBoxMesh(comm, user->dim, user->simplex, user->cells, NULL, NULL, NULL, PETSC_TRUE, dm);
    346:   {
    347:     DM               pdm = NULL;
    348:     PetscPartitioner part;

    350:     DMPlexGetPartitioner(*dm, &part);
    351:     PetscPartitionerSetFromOptions(part);
    352:     DMPlexDistribute(*dm, 0, NULL, &pdm);
    353:     if (pdm) {
    354:       DMDestroy(dm);
    355:       *dm  = pdm;
    356:     }
    357:   }
    358:   PetscStrcmp(user->dmType, DMPLEX, &flg);
    359:   if (flg) {
    360:     DM ndm;

    362:     DMConvert(*dm, user->dmType, &ndm);
    363:     if (ndm) {
    364:       DMDestroy(dm);
    365:       *dm  = ndm;
    366:     }
    367:   }
    368:   if (user->shear) {DMPlexShearGeometry(*dm, DM_X, NULL);}
    369:   DMLocalizeCoordinates(*dm);

    371:   PetscObjectSetName((PetscObject) *dm, "Mesh");
    372:   DMSetApplicationContext(*dm, user);
    373:   DMSetFromOptions(*dm);
    374:   DMViewFromOptions(*dm, NULL, "-dm_view");
    375:   return(0);
    376: }

    378: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
    379: {
    380:   PetscErrorCode (*exact)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
    381:   PetscDS        prob;
    382:   PetscInt       id;
    383:   PetscInt       dim;

    387:   DMGetDS(dm, &prob);
    388:   PetscDSGetSpatialDimension(prob, &dim);
    389:   switch (user->solType) {
    390:   case SOL_VLAP_QUADRATIC:
    391:     PetscDSSetResidual(prob, 0, f0_vlap_quadratic_u, f1_vlap_u);
    392:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_vlap_uu);
    393:     switch (dim) {
    394:     case 2: exact = quadratic_2d_u;break;
    395:     case 3: exact = quadratic_3d_u;break;
    396:     default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
    397:     }
    398:     break;
    399:   case SOL_ELAS_QUADRATIC:
    400:     PetscDSSetResidual(prob, 0, f0_elas_quadratic_u, f1_elas_u);
    401:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
    402:     switch (dim) {
    403:     case 2: exact = quadratic_2d_u;break;
    404:     case 3: exact = quadratic_3d_u;break;
    405:     default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
    406:     }
    407:     break;
    408:   case SOL_VLAP_TRIG:
    409:     PetscDSSetResidual(prob, 0, f0_vlap_trig_u, f1_vlap_u);
    410:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_vlap_uu);
    411:     switch (dim) {
    412:     case 2: exact = trig_2d_u;break;
    413:     case 3: exact = trig_3d_u;break;
    414:     default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
    415:     }
    416:     break;
    417:   case SOL_ELAS_TRIG:
    418:     PetscDSSetResidual(prob, 0, f0_elas_trig_u, f1_elas_u);
    419:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
    420:     switch (dim) {
    421:     case 2: exact = trig_2d_u;break;
    422:     case 3: exact = trig_3d_u;break;
    423:     default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
    424:     }
    425:     break;
    426:   case SOL_ELAS_AXIAL_DISP:
    427:     PetscDSSetResidual(prob, 0, NULL, f1_elas_u);
    428:     PetscDSSetBdResidual(prob, 0, f0_elas_axial_disp_bd_u, NULL);
    429:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
    430:     exact = axial_disp_u;
    431:     break;
    432:   case SOL_ELAS_UNIFORM_STRAIN:
    433:     PetscDSSetResidual(prob, 0, NULL, f1_elas_u);
    434:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
    435:     exact = uniform_strain_u;
    436:     break;
    437:   default: SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType);
    438:   }
    439:   PetscDSSetExactSolution(prob, 0, exact, user);
    440:   if (user->solType == SOL_ELAS_AXIAL_DISP) {
    441:     PetscInt cmp;

    443:     id   = dim == 3 ? 5 : 2;
    444:     DMAddBoundary(dm,   DM_BC_NATURAL,   "right",  "marker", 0, 0, NULL, (void (*)(void)) zero, NULL, 1, &id, user);
    445:     id   = dim == 3 ? 6 : 4;
    446:     cmp  = 0;
    447:     DMAddBoundary(dm,   DM_BC_ESSENTIAL, "left",   "marker", 0, 1, &cmp, (void (*)(void)) zero, NULL, 1, &id, user);
    448:     cmp  = dim == 3 ? 2 : 1;
    449:     id   = dim == 3 ? 1 : 1;
    450:     DMAddBoundary(dm,   DM_BC_ESSENTIAL, "bottom", "marker", 0, 1, &cmp, (void (*)(void)) zero, NULL, 1, &id, user);
    451:     if (dim == 3) {
    452:       cmp  = 1;
    453:       id   = 3;
    454:       DMAddBoundary(dm, DM_BC_ESSENTIAL, "front",  "marker", 0, 1, &cmp, (void (*)(void)) zero, NULL, 1, &id, user);
    455:     }
    456:   } else {
    457:     id = 1;
    458:     DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) exact, NULL, 1, &id, user);
    459:   }
    460:   return(0);
    461: }

    463: static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
    464: {

    468:   DMPlexCreateRigidBody(dm, origField, nullspace);
    469:   return(0);
    470: }

    472: PetscErrorCode SetupFE(DM dm, PetscInt Nc, PetscBool simplex, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx)
    473: {
    474:   AppCtx        *user = (AppCtx *) ctx;
    475:   DM             cdm  = dm;
    476:   PetscFE        fe;
    477:   char           prefix[PETSC_MAX_PATH_LEN];
    478:   PetscInt       dim;

    482:   /* Create finite element */
    483:   DMGetDimension(dm, &dim);
    484:   PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);
    485:   PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, Nc, simplex, name ? prefix : NULL, -1, &fe);
    486:   PetscObjectSetName((PetscObject) fe, name);
    487:   /* Set discretization and boundary conditions for each mesh */
    488:   DMSetField(dm, 0, NULL, (PetscObject) fe);
    489:   DMCreateDS(dm);
    490:   (*setup)(dm, user);
    491:   while (cdm) {
    492:     DMCopyDisc(dm, cdm);
    493:     if (user->useNearNullspace) {DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace);}
    494:     /* TODO: Check whether the boundary of coarse meshes is marked */
    495:     DMGetCoarseDM(cdm, &cdm);
    496:   }
    497:   PetscFEDestroy(&fe);
    498:   return(0);
    499: }

    501: int main(int argc, char **argv)
    502: {
    503:   DM             dm;   /* Problem specification */
    504:   SNES           snes; /* Nonlinear solver */
    505:   Vec            u;    /* Solutions */
    506:   AppCtx         user; /* User-defined work context */

    509:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
    510:   ProcessOptions(PETSC_COMM_WORLD, &user);
    511:   /* Primal system */
    512:   SNESCreate(PETSC_COMM_WORLD, &snes);
    513:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
    514:   SNESSetDM(snes, dm);
    515:   SetupFE(dm, user.dim, user.simplex, "displacement", SetupPrimalProblem, &user);
    516:   DMCreateGlobalVector(dm, &u);
    517:   VecSet(u, 0.0);
    518:   PetscObjectSetName((PetscObject) u, "displacement");
    519:   DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
    520:   SNESSetFromOptions(snes);
    521:   DMSNESCheckFromOptions(snes, u);
    522:   SNESSolve(snes, NULL, u);
    523:   SNESGetSolution(snes, &u);
    524:   VecViewFromOptions(u, NULL, "-displacement_view");
    525:   /* Cleanup */
    526:   VecDestroy(&u);
    527:   SNESDestroy(&snes);
    528:   DMDestroy(&dm);
    529:   PetscFinalize();
    530:   return ierr;
    531: }

    533: /*TEST

    535:   test:
    536:     suffix: 2d_p1_quad_vlap
    537:     requires: triangle
    538:     args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
    539:   test:
    540:     suffix: 2d_p2_quad_vlap
    541:     requires: triangle
    542:     args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
    543:   test:
    544:     suffix: 2d_p3_quad_vlap
    545:     requires: triangle
    546:     args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
    547:   test:
    548:     suffix: 2d_q1_quad_vlap
    549:     args: -simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
    550:   test:
    551:     suffix: 2d_q2_quad_vlap
    552:     args: -simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
    553:   test:
    554:     suffix: 2d_q3_quad_vlap
    555:     requires: !single
    556:     args: -simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
    557:   test:
    558:     suffix: 2d_p1_quad_elas
    559:     requires: triangle
    560:     args: -sol_type elas_quad -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
    561:   test:
    562:     suffix: 2d_p2_quad_elas
    563:     requires: triangle
    564:     args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001
    565:   test:
    566:     suffix: 2d_p3_quad_elas
    567:     requires: triangle
    568:     args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001
    569:   test:
    570:     suffix: 2d_q1_quad_elas
    571:     args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    572:   test:
    573:     suffix: 2d_q1_quad_elas_shear
    574:     args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    575:   test:
    576:     suffix: 2d_q2_quad_elas
    577:     args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001
    578:   test:
    579:     suffix: 2d_q2_quad_elas_shear
    580:     args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 2 -dmsnes_check
    581:   test:
    582:     suffix: 2d_q3_quad_elas
    583:     args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001
    584:   test:
    585:     suffix: 2d_q3_quad_elas_shear
    586:     requires: !single
    587:     args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 3 -dmsnes_check

    589:   test:
    590:     suffix: 3d_p1_quad_vlap
    591:     requires: ctetgen
    592:     args: -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    593:   test:
    594:     suffix: 3d_p2_quad_vlap
    595:     requires: ctetgen
    596:     args: -dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
    597:   test:
    598:     suffix: 3d_p3_quad_vlap
    599:     requires: ctetgen
    600:     args: -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
    601:   test:
    602:     suffix: 3d_q1_quad_vlap
    603:     args: -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    604:   test:
    605:     suffix: 3d_q2_quad_vlap
    606:     args: -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
    607:   test:
    608:     suffix: 3d_q3_quad_vlap
    609:     args: -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
    610:   test:
    611:     suffix: 3d_p1_quad_elas
    612:     requires: ctetgen
    613:     args: -sol_type elas_quad -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    614:   test:
    615:     suffix: 3d_p2_quad_elas
    616:     requires: ctetgen
    617:     args: -sol_type elas_quad -dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
    618:   test:
    619:     suffix: 3d_p3_quad_elas
    620:     requires: ctetgen
    621:     args: -sol_type elas_quad -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
    622:   test:
    623:     suffix: 3d_q1_quad_elas
    624:     args: -sol_type elas_quad -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    625:   test:
    626:     suffix: 3d_q2_quad_elas
    627:     args: -sol_type elas_quad -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
    628:   test:
    629:     suffix: 3d_q3_quad_elas
    630:     requires: !single
    631:     args: -sol_type elas_quad -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001

    633:   test:
    634:     suffix: 2d_p1_trig_vlap
    635:     requires: triangle
    636:     args: -sol_type vlap_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    637:   test:
    638:     suffix: 2d_p2_trig_vlap
    639:     requires: triangle
    640:     args: -sol_type vlap_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    641:   test:
    642:     suffix: 2d_p3_trig_vlap
    643:     requires: triangle
    644:     args: -sol_type vlap_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    645:   test:
    646:     suffix: 2d_q1_trig_vlap
    647:     args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    648:   test:
    649:     suffix: 2d_q2_trig_vlap
    650:     args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    651:   test:
    652:     suffix: 2d_q3_trig_vlap
    653:     args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    654:   test:
    655:     suffix: 2d_p1_trig_elas
    656:     requires: triangle
    657:     args: -sol_type elas_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    658:   test:
    659:     suffix: 2d_p2_trig_elas
    660:     requires: triangle
    661:     args: -sol_type elas_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    662:   test:
    663:     suffix: 2d_p3_trig_elas
    664:     requires: triangle
    665:     args: -sol_type elas_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    666:   test:
    667:     suffix: 2d_q1_trig_elas
    668:     args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    669:   test:
    670:     suffix: 2d_q1_trig_elas_shear
    671:     args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    672:   test:
    673:     suffix: 2d_q2_trig_elas
    674:     args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    675:   test:
    676:     suffix: 2d_q2_trig_elas_shear
    677:     args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    678:   test:
    679:     suffix: 2d_q3_trig_elas
    680:     args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
    681:   test:
    682:     suffix: 2d_q3_trig_elas_shear
    683:     args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate

    685:   test:
    686:     suffix: 3d_p1_trig_vlap
    687:     requires: ctetgen
    688:     args: -sol_type vlap_trig -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    689:   test:
    690:     suffix: 3d_p2_trig_vlap
    691:     requires: ctetgen
    692:     args: -sol_type vlap_trig -dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    693:   test:
    694:     suffix: 3d_p3_trig_vlap
    695:     requires: ctetgen
    696:     args: -sol_type vlap_trig -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    697:   test:
    698:     suffix: 3d_q1_trig_vlap
    699:     args: -sol_type vlap_trig -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    700:   test:
    701:     suffix: 3d_q2_trig_vlap
    702:     args: -sol_type vlap_trig -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    703:   test:
    704:     suffix: 3d_q3_trig_vlap
    705:     requires: !__float128
    706:     args: -sol_type vlap_trig -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    707:   test:
    708:     suffix: 3d_p1_trig_elas
    709:     requires: ctetgen
    710:     args: -sol_type elas_trig -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    711:   test:
    712:     suffix: 3d_p2_trig_elas
    713:     requires: ctetgen
    714:     args: -sol_type elas_trig -dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    715:   test:
    716:     suffix: 3d_p3_trig_elas
    717:     requires: ctetgen
    718:     args: -sol_type elas_trig -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    719:   test:
    720:     suffix: 3d_q1_trig_elas
    721:     args: -sol_type elas_trig -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
    722:   test:
    723:     suffix: 3d_q2_trig_elas
    724:     args: -sol_type elas_trig -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
    725:   test:
    726:     suffix: 3d_q3_trig_elas
    727:     requires: !__float128
    728:     args: -sol_type elas_trig -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate

    730:   test:
    731:     suffix: 2d_p1_axial_elas
    732:     requires: triangle
    733:     args: -sol_type elas_axial_disp -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    734:   test:
    735:     suffix: 2d_p2_axial_elas
    736:     requires: triangle
    737:     args: -sol_type elas_axial_disp -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
    738:   test:
    739:     suffix: 2d_p3_axial_elas
    740:     requires: triangle
    741:     args: -sol_type elas_axial_disp -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
    742:   test:
    743:     suffix: 2d_q1_axial_elas
    744:     args: -sol_type elas_axial_disp -simplex 0 -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 1 -dmsnes_check .0001 -pc_type lu
    745:   test:
    746:     suffix: 2d_q2_axial_elas
    747:     args: -sol_type elas_axial_disp -simplex 0 -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
    748:   test:
    749:     suffix: 2d_q3_axial_elas
    750:     args: -sol_type elas_axial_disp -simplex 0 -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu

    752:   test:
    753:     suffix: 2d_p1_uniform_elas
    754:     requires: triangle
    755:     args: -sol_type elas_uniform_strain -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    756:   test:
    757:     suffix: 2d_p2_uniform_elas
    758:     requires: triangle
    759:     args: -sol_type elas_uniform_strain -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    760:   test:
    761:     suffix: 2d_p3_uniform_elas
    762:     requires: triangle
    763:     args: -sol_type elas_uniform_strain -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    764:   test:
    765:     suffix: 2d_q1_uniform_elas
    766:     args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    767:   test:
    768:     suffix: 2d_q2_uniform_elas
    769:     requires: !single
    770:     args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
    771:   test:
    772:     suffix: 2d_q3_uniform_elas
    773:     requires: !single
    774:     args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu

    776: TEST*/