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 {
24:   SOL_VLAP_TRIG,
25:   SOL_ELAS_TRIG,
26:   SOL_ELAS_AXIAL_DISP,
27:   SOL_ELAS_UNIFORM_STRAIN,
28:   SOL_ELAS_GE,
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(PetscOptionsGetViewer(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(PetscOptionsRestoreViewer(&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) {
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;
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;
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:
680:       requires: triangle
681:       args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
682:     test:
684:       requires: triangle
685:       args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
686:     test:
688:       requires: triangle
689:       args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
690:     test:
692:       args: -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
693:     test:
695:       args: -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
696:     test:
698:       requires: !single
699:       args: -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
700:     test:
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:
706:       requires: triangle
707:       args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001
708:     test:
710:       requires: triangle
711:       args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001
712:     test:
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:
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:
720:       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001
721:     test:
723:       args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 2 -dmsnes_check
724:     test:
726:       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001
727:     test:
729:       requires: !single
730:       args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 3 -dmsnes_check

732:     test:
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:
738:       requires: ctetgen
739:       args: -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
740:     test:
742:       requires: ctetgen
743:       args: -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
744:     test:
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:
749:       args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
750:     test:
752:       args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
753:     test:
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:
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:
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:
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:
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:
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:
931:       args:
932:     test:
934:       args: -displacement_petscspace_degree 2
935:     test:
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*/