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 {SOL_VLAP_QUADRATIC, SOL_ELAS_QUADRATIC, SOL_VLAP_TRIG, SOL_ELAS_TRIG, SOL_ELAS_AXIAL_DISP, SOL_ELAS_UNIFORM_STRAIN, SOL_ELAS_GE, SOL_MASS_QUADRATIC, NUM_SOLUTION_TYPES} SolutionType;
 22: 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"};

 24: typedef enum {DEFORM_NONE, DEFORM_SHEAR, DEFORM_STEP, NUM_DEFORM_TYPES} DeformType;
 25: const char *deformTypes[NUM_DEFORM_TYPES+1] = {"none", "shear", "step", "unknown"};

 27: typedef struct {
 28:   PetscScalar mu;     /* shear modulus */
 29:   PetscScalar lambda; /* Lame's first parameter */
 30: } Parameter;

 32: typedef struct {
 33:   /* Domain and mesh definition */
 34:   char         dmType[256]; /* DM type for the solve */
 35:   DeformType   deform;      /* Domain deformation type */
 36:   /* Problem definition */
 37:   SolutionType solType;     /* Type of exact solution */
 38:   PetscBag     bag;         /* Problem parameters */
 39:   /* Solver definition */
 40:   PetscBool    useNearNullspace; /* Use the rigid body modes as a near nullspace for AMG */
 41: } AppCtx;

 43: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 44: {
 45:   PetscInt d;
 46:   for (d = 0; d < dim; ++d) u[d] = 0.0;
 47:   return 0;
 48: }

 50: static PetscErrorCode ge_shift(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 51: {
 52:   PetscInt d;
 53:   u[0] = 0.1;
 54:   for (d = 1; d < dim; ++d) u[d] = 0.0;
 55:   return 0;
 56: }

 58: static PetscErrorCode quadratic_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 59: {
 60:   u[0] = x[0]*x[0];
 61:   u[1] = x[1]*x[1] - 2.0*x[0]*x[1];
 62:   return 0;
 63: }

 65: static PetscErrorCode quadratic_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 66: {
 67:   u[0] = x[0]*x[0];
 68:   u[1] = x[1]*x[1] - 2.0*x[0]*x[1];
 69:   u[2] = x[2]*x[2] - 2.0*x[1]*x[2];
 70:   return 0;
 71: }

 73: /*
 74:   u = x^2
 75:   v = y^2 - 2xy
 76:   Delta <u,v> - f = <2, 2> - <2, 2>

 78:   u = x^2
 79:   v = y^2 - 2xy
 80:   w = z^2 - 2yz
 81:   Delta <u,v,w> - f = <2, 2, 2> - <2, 2, 2>
 82: */
 83: static void f0_vlap_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 84:                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 85:                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 86:                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 87: {
 88:   PetscInt d;
 89:   for (d = 0; d < dim; ++d) f0[d] += 2.0;
 90: }

 92: /*
 93:   u = x^2
 94:   v = y^2 - 2xy
 95:   \varepsilon = / 2x     -y    \
 96:                 \ -y   2y - 2x /
 97:   Tr(\varepsilon) = div u = 2y
 98:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
 99:     = \lambda \partial_j (2y) + 2\mu < 2-1, 2 >
100:     = \lambda < 0, 2 > + \mu < 2, 4 >

102:   u = x^2
103:   v = y^2 - 2xy
104:   w = z^2 - 2yz
105:   \varepsilon = / 2x     -y       0   \
106:                 | -y   2y - 2x   -z   |
107:                 \  0     -z    2z - 2y/
108:   Tr(\varepsilon) = div u = 2z
109:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
110:     = \lambda \partial_j (2z) + 2\mu < 2-1, 2-1, 2 >
111:     = \lambda < 0, 0, 2 > + \mu < 2, 2, 4 >
112: */
113: static void f0_elas_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
114:                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
115:                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
116:                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
117: {
118:   const PetscReal mu     = 1.0;
119:   const PetscReal lambda = 1.0;
120:   PetscInt        d;

122:   for (d = 0; d < dim-1; ++d) f0[d] += 2.0*mu;
123:   f0[dim-1] += 2.0*lambda + 4.0*mu;
124: }

126: static void f0_mass_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
127:                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
128:                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
129:                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
130: {
131:   if (dim == 2) {
132:     f0[0] -= x[0]*x[0];
133:     f0[1] -= x[1]*x[1] - 2.0*x[0]*x[1];
134:   } else {
135:     f0[0] -= x[0]*x[0];
136:     f0[1] -= x[1]*x[1] - 2.0*x[0]*x[1];
137:     f0[2] -= x[2]*x[2] - 2.0*x[1]*x[2];
138:   }
139: }

141: static PetscErrorCode trig_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
142: {
143:   u[0] = PetscSinReal(2.0*PETSC_PI*x[0]);
144:   u[1] = PetscSinReal(2.0*PETSC_PI*x[1]) - 2.0*x[0]*x[1];
145:   return 0;
146: }

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

156: /*
157:   u = sin(2 pi x)
158:   v = sin(2 pi y) - 2xy
159:   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)>

161:   u = sin(2 pi x)
162:   v = sin(2 pi y) - 2xy
163:   w = sin(2 pi z) - 2yz
164:   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)>
165: */
166: static void f0_vlap_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
167:                            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
168:                            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
169:                            PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
170: {
171:   PetscInt d;
172:   for (d = 0; d < dim; ++d) f0[d] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]);
173: }

175: /*
176:   u = sin(2 pi x)
177:   v = sin(2 pi y) - 2xy
178:   \varepsilon = / 2 pi cos(2 pi x)             -y        \
179:                 \      -y          2 pi cos(2 pi y) - 2x /
180:   Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
181:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
182:     = \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) >
183:     = \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) >

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

206:   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);
207: }

209: static PetscErrorCode axial_disp_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
210: {
211:   const PetscReal mu     = 1.0;
212:   const PetscReal lambda = 1.0;
213:   const PetscReal N      = 1.0;
214:   PetscInt d;

216:   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];
217:   u[1] = -0.25*lambda/mu/(lambda+mu)*N*x[1];
218:   for (d = 2; d < dim; ++d) u[d] = 0.0;
219:   return 0;
220: }

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

226:      n_i \sigma_{ij} = t_i

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

240:   u = x - 1/2
241:   v = 0
242:   w = 0
243:   \varepsilon = / x  0  0 \
244:                 | 0  0  0 |
245:                 \ 0  0  0 /
246:   Tr(\varepsilon) = div u = x
247:   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
248:     = \lambda \partial_j x + 2\mu < 1, 0, 0 >
249:     = \lambda < 1, 0, 0 > + \mu < 2, 0, 0 >
250: */
251: static void f0_elas_axial_disp_bd_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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
255: {
256:   const PetscReal N = -1.0;

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

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

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

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

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

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

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

296: static void f1_elas_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
297:                       const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
298:                       const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
299:                       PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
300: {
301:   const PetscInt  Nc     = dim;
302:   const PetscReal mu     = 1.0;
303:   const PetscReal lambda = 1.0;
304:   PetscInt        c, d;

306:   for (c = 0; c < Nc; ++c) {
307:     for (d = 0; d < dim; ++d) {
308:       f1[c*dim+d] += mu*(u_x[c*dim+d] + u_x[d*dim+c]);
309:       f1[c*dim+c] += lambda*u_x[d*dim+d];
310:     }
311:   }
312: }

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

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

325: static void g3_vlap_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
326:                        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
327:                        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
328:                        PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
329: {
330:   const PetscInt Nc = dim;
331:   PetscInt       c, d;

333:   for (c = 0; c < Nc; ++c) {
334:     for (d = 0; d < dim; ++d) {
335:       g3[((c*Nc + c)*dim + d)*dim + d] = 1.0;
336:     }
337:   }
338: }

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

343:   \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg}
344:   = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc
345: */
346: static void g3_elas_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
347:                        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
348:                        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
349:                        PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
350: {
351:   const PetscInt  Nc     = dim;
352:   const PetscReal mu     = 1.0;
353:   const PetscReal lambda = 1.0;
354:   PetscInt        c, d;

356:   for (c = 0; c < Nc; ++c) {
357:     for (d = 0; d < dim; ++d) {
358:       g3[((c*Nc + c)*dim + d)*dim + d] += mu;
359:       g3[((c*Nc + d)*dim + d)*dim + c] += mu;
360:       g3[((c*Nc + d)*dim + c)*dim + d] += lambda;
361:     }
362:   }
363: }

365: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
366: {
367:   PetscInt       sol = 0, def = 0;

371:   options->deform           = DEFORM_NONE;
372:   options->solType          = SOL_VLAP_QUADRATIC;
373:   options->useNearNullspace = PETSC_TRUE;
374:   PetscStrncpy(options->dmType, DMPLEX, 256);

376:   PetscOptionsBegin(comm, "", "Linear Elasticity Problem Options", "DMPLEX");
377:   PetscOptionsEList("-deform_type", "Type of domain deformation", "ex17.c", deformTypes, NUM_DEFORM_TYPES, deformTypes[options->deform], &def, NULL);
378:   options->deform = (DeformType) def;
379:   PetscOptionsEList("-sol_type", "Type of exact solution", "ex17.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL);
380:   options->solType = (SolutionType) sol;
381:   PetscOptionsBool("-near_nullspace", "Use the rigid body modes as an AMG near nullspace", "ex17.c", options->useNearNullspace, &options->useNearNullspace, NULL);
382:   PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex17.c", DMList, options->dmType, options->dmType, 256, NULL);
383:   PetscOptionsEnd();
384:   return 0;
385: }

387: static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx)
388: {
389:   PetscBag       bag;
390:   Parameter     *p;

393:   /* setup PETSc parameter bag */
394:   PetscBagGetData(ctx->bag,(void**)&p);
395:   PetscBagSetName(ctx->bag,"par","Elastic Parameters");
396:   bag  = ctx->bag;
397:   PetscBagRegisterScalar(bag, &p->mu,     1.0, "mu",     "Shear Modulus, Pa");
398:   PetscBagRegisterScalar(bag, &p->lambda, 1.0, "lambda", "Lame's first parameter, Pa");
399:   PetscBagSetFromOptions(bag);
400:   {
401:     PetscViewer       viewer;
402:     PetscViewerFormat format;
403:     PetscBool         flg;

405:     PetscOptionsGetViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg);
406:     if (flg) {
407:       PetscViewerPushFormat(viewer, format);
408:       PetscBagView(bag, viewer);
409:       PetscViewerFlush(viewer);
410:       PetscViewerPopFormat(viewer);
411:       PetscViewerDestroy(&viewer);
412:     }
413:   }
414:   return 0;
415: }

417: static PetscErrorCode DMPlexDistortGeometry(DM dm)
418: {
419:   DM             cdm;
420:   DMLabel        label;
421:   Vec            coordinates;
422:   PetscScalar   *coords;
423:   PetscReal      mid = 0.5;
424:   PetscInt       cdim, d, vStart, vEnd, v;

427:   DMGetCoordinateDM(dm, &cdm);
428:   DMGetCoordinateDim(dm, &cdim);
429:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
430:   DMGetLabel(dm, "marker", &label);
431:   DMGetCoordinatesLocal(dm, &coordinates);
432:   VecGetArrayWrite(coordinates, &coords);
433:   for (v = vStart; v < vEnd; ++v) {
434:     PetscScalar *pcoords, shift;
435:     PetscInt     val;

437:     DMLabelGetValue(label, v, &val);
438:     if (val >= 0) continue;
439:     DMPlexPointLocalRef(cdm, v, coords, &pcoords);
440:     shift = 0.2 * PetscAbsScalar(pcoords[0] - mid);
441:     shift = PetscRealPart(pcoords[0]) > mid ? shift : -shift;
442:     for (d = 1; d < cdim; ++d) pcoords[d] += shift;
443:   }
444:   VecRestoreArrayWrite(coordinates, &coords);
445:   return 0;
446: }

448: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
449: {
451:   DMCreate(comm, dm);
452:   DMSetType(*dm, DMPLEX);
453:   DMSetFromOptions(*dm);
454:   switch (user->deform) {
455:     case DEFORM_NONE:  break;
456:     case DEFORM_SHEAR: DMPlexShearGeometry(*dm, DM_X, NULL);break;
457:     case DEFORM_STEP:  DMPlexDistortGeometry(*dm);break;
458:     default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid deformation type: %s (%D)", deformTypes[PetscMin(user->deform, NUM_DEFORM_TYPES)], user->deform);
459:   }
460:   DMSetApplicationContext(*dm, user);
461:   DMViewFromOptions(*dm, NULL, "-dm_view");
462:   return 0;
463: }

465: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
466: {
467:   PetscErrorCode (*exact)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
468:   Parameter       *param;
469:   PetscDS          ds;
470:   PetscWeakForm    wf;
471:   DMLabel          label;
472:   PetscInt         id, bd;
473:   PetscInt         dim;

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

554:     id   = dim == 3 ? 6 : 4;
555:     cmp  = 0;
556:     DMAddBoundary(dm,   DM_BC_ESSENTIAL, "left",   label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);
557:     cmp  = dim == 3 ? 2 : 1;
558:     id   = dim == 3 ? 1 : 1;
559:     DMAddBoundary(dm,   DM_BC_ESSENTIAL, "bottom", label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);
560:     if (dim == 3) {
561:       cmp  = 1;
562:       id   = 3;
563:       DMAddBoundary(dm, DM_BC_ESSENTIAL, "front",  label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);
564:     }
565:   } else if (user->solType == SOL_ELAS_GE) {
566:     PetscInt cmp;

568:     id   = dim == 3 ? 6 : 4;
569:     DMAddBoundary(dm,   DM_BC_ESSENTIAL, "left",   label, 1, &id, 0, 0, NULL, (void (*)(void)) zero, NULL, user, NULL);
570:     id   = dim == 3 ? 5 : 2;
571:     cmp  = 0;
572:     DMAddBoundary(dm,   DM_BC_ESSENTIAL, "right",  label, 1, &id, 0, 1, &cmp, (void (*)(void)) ge_shift, NULL, user, NULL);
573:   } else {
574:     id = 1;
575:     DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) exact, NULL, user, NULL);
576:   }
577:   /* Setup constants */
578:   {
579:     PetscScalar constants[2];

581:     constants[0] = param->mu;     /* shear modulus, Pa */
582:     constants[1] = param->lambda; /* Lame's first parameter, Pa */
583:     PetscDSSetConstants(ds, 2, constants);
584:   }
585:   return 0;
586: }

588: static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
589: {
590:   DMPlexCreateRigidBody(dm, origField, nullspace);
591:   return 0;
592: }

594: PetscErrorCode SetupFE(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx)
595: {
596:   AppCtx        *user = (AppCtx *) ctx;
597:   DM             cdm  = dm;
598:   PetscFE        fe;
599:   char           prefix[PETSC_MAX_PATH_LEN];
600:   DMPolytopeType ct;
601:   PetscBool      simplex;
602:   PetscInt       dim, cStart;

604:   /* Create finite element */
605:   DMGetDimension(dm, &dim);
606:   DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
607:   DMPlexGetCellType(dm, cStart, &ct);
608:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
609:   PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);
610:   PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, dim, simplex, name ? prefix : NULL, -1, &fe);
611:   PetscObjectSetName((PetscObject) fe, name);
612:   /* Set discretization and boundary conditions for each mesh */
613:   DMSetField(dm, 0, NULL, (PetscObject) fe);
614:   DMCreateDS(dm);
615:   (*setup)(dm, user);
616:   while (cdm) {
617:     DMCopyDisc(dm, cdm);
618:     if (user->useNearNullspace) DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace);
619:     /* TODO: Check whether the boundary of coarse meshes is marked */
620:     DMGetCoarseDM(cdm, &cdm);
621:   }
622:   PetscFEDestroy(&fe);
623:   return 0;
624: }

626: int main(int argc, char **argv)
627: {
628:   DM             dm;   /* Problem specification */
629:   SNES           snes; /* Nonlinear solver */
630:   Vec            u;    /* Solutions */
631:   AppCtx         user; /* User-defined work context */

633:   PetscInitialize(&argc, &argv, NULL,help);
634:   ProcessOptions(PETSC_COMM_WORLD, &user);
635:   PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &user.bag);
636:   SetupParameters(PETSC_COMM_WORLD, &user);
637:   /* Primal system */
638:   SNESCreate(PETSC_COMM_WORLD, &snes);
639:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
640:   SNESSetDM(snes, dm);
641:   SetupFE(dm, "displacement", SetupPrimalProblem, &user);
642:   DMCreateGlobalVector(dm, &u);
643:   VecSet(u, 0.0);
644:   PetscObjectSetName((PetscObject) u, "displacement");
645:   DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
646:   SNESSetFromOptions(snes);
647:   DMSNESCheckFromOptions(snes, u);
648:   SNESSolve(snes, NULL, u);
649:   SNESGetSolution(snes, &u);
650:   VecViewFromOptions(u, NULL, "-displacement_view");
651:   /* Cleanup */
652:   VecDestroy(&u);
653:   SNESDestroy(&snes);
654:   DMDestroy(&dm);
655:   PetscBagDestroy(&user.bag);
656:   PetscFinalize();
657:   return 0;
658: }

660: /*TEST

662:   testset:
663:     args: -dm_plex_box_faces 1,1,1

665:     test:
666:       suffix: 2d_p1_quad_vlap
667:       requires: triangle
668:       args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
669:     test:
670:       suffix: 2d_p2_quad_vlap
671:       requires: triangle
672:       args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
673:     test:
674:       suffix: 2d_p3_quad_vlap
675:       requires: triangle
676:       args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
677:     test:
678:       suffix: 2d_q1_quad_vlap
679:       args: -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
680:     test:
681:       suffix: 2d_q2_quad_vlap
682:       args: -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
683:     test:
684:       suffix: 2d_q3_quad_vlap
685:       requires: !single
686:       args: -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
687:     test:
688:       suffix: 2d_p1_quad_elas
689:       requires: triangle
690:       args: -sol_type elas_quad -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
691:     test:
692:       suffix: 2d_p2_quad_elas
693:       requires: triangle
694:       args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001
695:     test:
696:       suffix: 2d_p3_quad_elas
697:       requires: triangle
698:       args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001
699:     test:
700:       suffix: 2d_q1_quad_elas
701:       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
702:     test:
703:       suffix: 2d_q1_quad_elas_shear
704:       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
705:     test:
706:       suffix: 2d_q2_quad_elas
707:       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001
708:     test:
709:       suffix: 2d_q2_quad_elas_shear
710:       args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 2 -dmsnes_check
711:     test:
712:       suffix: 2d_q3_quad_elas
713:       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001
714:     test:
715:       suffix: 2d_q3_quad_elas_shear
716:       requires: !single
717:       args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 3 -dmsnes_check

719:     test:
720:       suffix: 3d_p1_quad_vlap
721:       requires: ctetgen
722:       args: -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
723:     test:
724:       suffix: 3d_p2_quad_vlap
725:       requires: ctetgen
726:       args: -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
727:     test:
728:       suffix: 3d_p3_quad_vlap
729:       requires: ctetgen
730:       args: -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
731:     test:
732:       suffix: 3d_q1_quad_vlap
733:       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
734:     test:
735:       suffix: 3d_q2_quad_vlap
736:       args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
737:     test:
738:       suffix: 3d_q3_quad_vlap
739:       args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
740:     test:
741:       suffix: 3d_p1_quad_elas
742:       requires: ctetgen
743:       args: -sol_type elas_quad -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
744:     test:
745:       suffix: 3d_p2_quad_elas
746:       requires: ctetgen
747:       args: -sol_type elas_quad -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
748:     test:
749:       suffix: 3d_p3_quad_elas
750:       requires: ctetgen
751:       args: -sol_type elas_quad -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
752:     test:
753:       suffix: 3d_q1_quad_elas
754:       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
755:     test:
756:       suffix: 3d_q2_quad_elas
757:       args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
758:     test:
759:       suffix: 3d_q3_quad_elas
760:       requires: !single
761:       args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001

763:     test:
764:       suffix: 2d_p1_trig_vlap
765:       requires: triangle
766:       args: -sol_type vlap_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
767:     test:
768:       suffix: 2d_p2_trig_vlap
769:       requires: triangle
770:       args: -sol_type vlap_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
771:     test:
772:       suffix: 2d_p3_trig_vlap
773:       requires: triangle
774:       args: -sol_type vlap_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
775:     test:
776:       suffix: 2d_q1_trig_vlap
777:       args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
778:     test:
779:       suffix: 2d_q2_trig_vlap
780:       args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
781:     test:
782:       suffix: 2d_q3_trig_vlap
783:       args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
784:     test:
785:       suffix: 2d_p1_trig_elas
786:       requires: triangle
787:       args: -sol_type elas_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
788:     test:
789:       suffix: 2d_p2_trig_elas
790:       requires: triangle
791:       args: -sol_type elas_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
792:     test:
793:       suffix: 2d_p3_trig_elas
794:       requires: triangle
795:       args: -sol_type elas_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
796:     test:
797:       suffix: 2d_q1_trig_elas
798:       args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
799:     test:
800:       suffix: 2d_q1_trig_elas_shear
801:       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
802:     test:
803:       suffix: 2d_q2_trig_elas
804:       args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
805:     test:
806:       suffix: 2d_q2_trig_elas_shear
807:       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
808:     test:
809:       suffix: 2d_q3_trig_elas
810:       args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
811:     test:
812:       suffix: 2d_q3_trig_elas_shear
813:       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

815:     test:
816:       suffix: 3d_p1_trig_vlap
817:       requires: ctetgen
818:       args: -sol_type vlap_trig -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
819:     test:
820:       suffix: 3d_p2_trig_vlap
821:       requires: ctetgen
822:       args: -sol_type vlap_trig -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
823:     test:
824:       suffix: 3d_p3_trig_vlap
825:       requires: ctetgen
826:       args: -sol_type vlap_trig -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
827:     test:
828:       suffix: 3d_q1_trig_vlap
829:       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
830:     test:
831:       suffix: 3d_q2_trig_vlap
832:       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
833:     test:
834:       suffix: 3d_q3_trig_vlap
835:       requires: !__float128
836:       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
837:     test:
838:       suffix: 3d_p1_trig_elas
839:       requires: ctetgen
840:       args: -sol_type elas_trig -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
841:     test:
842:       suffix: 3d_p2_trig_elas
843:       requires: ctetgen
844:       args: -sol_type elas_trig -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
845:     test:
846:       suffix: 3d_p3_trig_elas
847:       requires: ctetgen
848:       args: -sol_type elas_trig -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
849:     test:
850:       suffix: 3d_q1_trig_elas
851:       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
852:     test:
853:       suffix: 3d_q2_trig_elas
854:       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
855:     test:
856:       suffix: 3d_q3_trig_elas
857:       requires: !__float128
858:       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

860:     test:
861:       suffix: 2d_p1_axial_elas
862:       requires: triangle
863:       args: -sol_type elas_axial_disp -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 2 -dmsnes_check .0001 -pc_type lu
864:     test:
865:       suffix: 2d_p2_axial_elas
866:       requires: triangle
867:       args: -sol_type elas_axial_disp -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
868:     test:
869:       suffix: 2d_p3_axial_elas
870:       requires: triangle
871:       args: -sol_type elas_axial_disp -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
872:     test:
873:       suffix: 2d_q1_axial_elas
874:       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
875:     test:
876:       suffix: 2d_q2_axial_elas
877:       args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type   lu
878:     test:
879:       suffix: 2d_q3_axial_elas
880:       args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type   lu

882:     test:
883:       suffix: 2d_p1_uniform_elas
884:       requires: triangle
885:       args: -sol_type elas_uniform_strain -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
886:     test:
887:       suffix: 2d_p2_uniform_elas
888:       requires: triangle
889:       args: -sol_type elas_uniform_strain -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
890:     test:
891:       suffix: 2d_p3_uniform_elas
892:       requires: triangle
893:       args: -sol_type elas_uniform_strain -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
894:     test:
895:       suffix: 2d_q1_uniform_elas
896:       args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
897:     test:
898:       suffix: 2d_q2_uniform_elas
899:       requires: !single
900:       args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
901:     test:
902:       suffix: 2d_q3_uniform_elas
903:       requires: !single
904:       args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
905:     test:
906:       suffix: 2d_p1_uniform_elas_step
907:       requires: triangle
908:       args: -sol_type elas_uniform_strain -deform_type step -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu

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

913:     test:
914:       suffix: 2d_q1_uniform_elas_step
915:       args: -sol_type elas_uniform_strain -dm_refine 2
916:     test:
917:       suffix: 2d_q1_quad_vlap_step
918:       args:
919:     test:
920:       suffix: 2d_q2_quad_vlap_step
921:       args: -displacement_petscspace_degree 2
922:     test:
923:       suffix: 2d_q1_quad_mass_step
924:       args: -sol_type mass_quad

926:   testset:
927:     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 \
928:           -dm_plex_box_faces 5,5,2 -dm_plex_separate_marker -dm_refine 0 -petscpartitioner_type simple \
929:           -sol_type elas_ge \
930:           -snes_max_it 2 -snes_rtol 1.e-10 \
931:           -ksp_type cg -ksp_rtol 1.e-10 -ksp_max_it 100 -ksp_norm_type unpreconditioned \
932:           -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \
933:             -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true \
934:             -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 \
935:             -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 \
936:             -matptap_via scalable
937:     test:
938:       suffix: ge_q1_0
939:       args: -displacement_petscspace_degree 1

941: TEST*/