xref: /petsc/src/snes/tutorials/ex17.c (revision 34ad990445da331e7a633c30db32e380550a8b3a)
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";
6 
7 /*
8   https://en.wikipedia.org/wiki/Linear_elasticity
9 */
10 
11 #include <petscdmplex.h>
12 #include <petscsnes.h>
13 #include <petscds.h>
14 #include <petscconvest.h>
15 
16 typedef enum {SOL_VLAP_QUADRATIC, SOL_ELAS_QUADRATIC, SOL_VLAP_TRIG, SOL_ELAS_TRIG, SOL_ELAS_AXIAL_DISP, SOL_ELAS_UNIFORM_STRAIN, SOL_MASS_QUADRATIC, 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", "mass_quad", "unknown"};
18 
19 typedef enum {DEFORM_NONE, DEFORM_SHEAR, DEFORM_STEP, NUM_DEFORM_TYPES} DeformType;
20 const char *deformTypes[NUM_DEFORM_TYPES+1] = {"none", "shear", "step", "unknown"};
21 
22 typedef struct {
23   /* Domain and mesh definition */
24   char         dmType[256]; /* DM type for the solve */
25   DeformType   deform;      /* Domain deformation type */
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;
31 
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 }
38 
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 }
45 
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 }
53 
54 /*
55   u = x^2
56   v = y^2 - 2xy
57   Delta <u,v> - f = <2, 2> - <2, 2>
58 
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 }
72 
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 >
82 
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;
102 
103   for (d = 0; d < dim-1; ++d) f0[d] += 2.0*mu;
104   f0[dim-1] += 2.0*lambda + 4.0*mu;
105 }
106 
107 static void f0_mass_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
108                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
109                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
110                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
111 {
112   if (dim == 2) {
113     f0[0] -= x[0]*x[0];
114     f0[1] -= x[1]*x[1] - 2.0*x[0]*x[1];
115   } else {
116     f0[0] -= x[0]*x[0];
117     f0[1] -= x[1]*x[1] - 2.0*x[0]*x[1];
118     f0[2] -= x[2]*x[2] - 2.0*x[1]*x[2];
119   }
120 }
121 
122 static PetscErrorCode trig_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
123 {
124   u[0] = PetscSinReal(2.0*PETSC_PI*x[0]);
125   u[1] = PetscSinReal(2.0*PETSC_PI*x[1]) - 2.0*x[0]*x[1];
126   return 0;
127 }
128 
129 static PetscErrorCode trig_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
130 {
131   u[0] = PetscSinReal(2.0*PETSC_PI*x[0]);
132   u[1] = PetscSinReal(2.0*PETSC_PI*x[1]) - 2.0*x[0]*x[1];
133   u[2] = PetscSinReal(2.0*PETSC_PI*x[2]) - 2.0*x[1]*x[2];
134   return 0;
135 }
136 
137 /*
138   u = sin(2 pi x)
139   v = sin(2 pi y) - 2xy
140   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)>
141 
142   u = sin(2 pi x)
143   v = sin(2 pi y) - 2xy
144   w = sin(2 pi z) - 2yz
145   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)>
146 */
147 static void f0_vlap_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
148                            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
149                            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
150                            PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
151 {
152   PetscInt d;
153   for (d = 0; d < dim; ++d) f0[d] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]);
154 }
155 
156 /*
157   u = sin(2 pi x)
158   v = sin(2 pi y) - 2xy
159   \varepsilon = / 2 pi cos(2 pi x)             -y        \
160                 \      -y          2 pi cos(2 pi y) - 2x /
161   Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
162   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
163     = \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) >
164     = \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) >
165 
166   u = sin(2 pi x)
167   v = sin(2 pi y) - 2xy
168   w = sin(2 pi z) - 2yz
169   \varepsilon = / 2 pi cos(2 pi x)            -y                     0         \
170                 |         -y       2 pi cos(2 pi y) - 2x            -z         |
171                 \          0                  -z         2 pi cos(2 pi z) - 2y /
172   Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y
173   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
174     = \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) >
175     = \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) >
176 */
177 static void f0_elas_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
178                            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
179                            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
180                            PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
181 {
182   const PetscReal mu     = 1.0;
183   const PetscReal lambda = 1.0;
184   const PetscReal fact   = 4.0*PetscSqr(PETSC_PI);
185   PetscInt        d;
186 
187   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);
188 }
189 
190 static PetscErrorCode axial_disp_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
191 {
192   const PetscReal mu     = 1.0;
193   const PetscReal lambda = 1.0;
194   const PetscReal N      = 1.0;
195   PetscInt d;
196 
197   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];
198   u[1] = -0.25*lambda/mu/(lambda+mu)*N*x[1];
199   for (d = 2; d < dim; ++d) u[d] = 0.0;
200   return 0;
201 }
202 
203 /*
204   We will pull/push on the right side of a block of linearly elastic material. The uniform traction conditions on the
205   right side of the box will result in a uniform strain along x and y. The Neumann BC is given by
206 
207      n_i \sigma_{ij} = t_i
208 
209   u = (1/(2\mu) - 1) x
210   v = -y
211   f = 0
212   t = <4\mu/\lambda (\lambda + \mu), 0>
213   \varepsilon = / 1/(2\mu) - 1   0 \
214                 \ 0             -1 /
215   Tr(\varepsilon) = div u = 1/(2\mu) - 2
216   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
217     = \lambda \partial_j (1/(2\mu) - 2) + 2\mu < 0, 0 >
218     = \lambda < 0, 0 > + \mu < 0, 0 > = 0
219   NBC =  <1,0> . <4\mu/\lambda (\lambda + \mu), 0> = 4\mu/\lambda (\lambda + \mu)
220 
221   u = x - 1/2
222   v = 0
223   w = 0
224   \varepsilon = / x  0  0 \
225                 | 0  0  0 |
226                 \ 0  0  0 /
227   Tr(\varepsilon) = div u = x
228   div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
229     = \lambda \partial_j x + 2\mu < 1, 0, 0 >
230     = \lambda < 1, 0, 0 > + \mu < 2, 0, 0 >
231 */
232 static void f0_elas_axial_disp_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
233                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
234                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
235                                     PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
236 {
237   const PetscReal N = -1.0;
238 
239   f0[0] = N;
240 }
241 
242 static PetscErrorCode uniform_strain_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
243 {
244   const PetscReal eps_xx = 0.1;
245   const PetscReal eps_xy = 0.3;
246   const PetscReal eps_yy = 0.25;
247   PetscInt d;
248 
249   u[0] = eps_xx*x[0] + eps_xy*x[1];
250   u[1] = eps_xy*x[0] + eps_yy*x[1];
251   for (d = 2; d < dim; ++d) u[d] = 0.0;
252   return 0;
253 }
254 
255 static void f0_mass_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
256                       const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
257                       const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
258                       PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
259 {
260   const PetscInt Nc = dim;
261   PetscInt       c;
262 
263   for (c = 0; c < Nc; ++c) f0[c] = u[c];
264 }
265 
266 static void f1_vlap_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
267                       const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
268                       const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
269                       PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
270 {
271   const PetscInt Nc = dim;
272   PetscInt       c, d;
273 
274   for (c = 0; c < Nc; ++c) for (d = 0; d < dim; ++d) f1[c*dim+d] += u_x[c*dim+d];
275 }
276 
277 static void f1_elas_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
278                       const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
279                       const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
280                       PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
281 {
282   const PetscInt  Nc     = dim;
283   const PetscReal mu     = 1.0;
284   const PetscReal lambda = 1.0;
285   PetscInt        c, d;
286 
287   for (c = 0; c < Nc; ++c) {
288     for (d = 0; d < dim; ++d) {
289       f1[c*dim+d] += mu*(u_x[c*dim+d] + u_x[d*dim+c]);
290       f1[c*dim+c] += lambda*u_x[d*dim+d];
291     }
292   }
293 }
294 
295 static void g0_mass_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
296                        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
297                        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
298                        PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
299 {
300   const PetscInt Nc = dim;
301   PetscInt       c;
302 
303   for (c = 0; c < Nc; ++c) g0[c*Nc + c] = 1.0;
304 }
305 
306 static void g3_vlap_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
307                        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
308                        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
309                        PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
310 {
311   const PetscInt Nc = dim;
312   PetscInt       c, d;
313 
314   for (c = 0; c < Nc; ++c) {
315     for (d = 0; d < dim; ++d) {
316       g3[((c*Nc + c)*dim + d)*dim + d] = 1.0;
317     }
318   }
319 }
320 
321 /*
322   \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc
323 
324   \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg}
325   = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc
326 */
327 static void g3_elas_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
328                        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
329                        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
330                        PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
331 {
332   const PetscInt  Nc     = dim;
333   const PetscReal mu     = 1.0;
334   const PetscReal lambda = 1.0;
335   PetscInt        c, d;
336 
337   for (c = 0; c < Nc; ++c) {
338     for (d = 0; d < dim; ++d) {
339       g3[((c*Nc + c)*dim + d)*dim + d] += mu;
340       g3[((c*Nc + d)*dim + d)*dim + c] += mu;
341       g3[((c*Nc + d)*dim + c)*dim + d] += lambda;
342     }
343   }
344 }
345 
346 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
347 {
348   PetscInt       sol = 0, def = 0;
349   PetscErrorCode ierr;
350 
351   PetscFunctionBeginUser;
352   options->deform           = DEFORM_NONE;
353   options->solType          = SOL_VLAP_QUADRATIC;
354   options->useNearNullspace = PETSC_TRUE;
355   ierr = PetscStrncpy(options->dmType, DMPLEX, 256);CHKERRQ(ierr);
356 
357   ierr = PetscOptionsBegin(comm, "", "Linear Elasticity Problem Options", "DMPLEX");CHKERRQ(ierr);
358   ierr = PetscOptionsEList("-deform_type", "Type of domain deformation", "ex17.c", deformTypes, NUM_DEFORM_TYPES, deformTypes[options->deform], &def, NULL);CHKERRQ(ierr);
359   options->deform = (DeformType) def;
360   ierr = PetscOptionsEList("-sol_type", "Type of exact solution", "ex17.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL);CHKERRQ(ierr);
361   options->solType = (SolutionType) sol;
362   ierr = PetscOptionsBool("-near_nullspace", "Use the rigid body modes as an AMG near nullspace", "ex17.c", options->useNearNullspace, &options->useNearNullspace, NULL);CHKERRQ(ierr);
363   ierr = PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex17.c", DMList, options->dmType, options->dmType, 256, NULL);CHKERRQ(ierr);
364   ierr = PetscOptionsEnd();CHKERRQ(ierr);
365   PetscFunctionReturn(0);
366 }
367 
368 static PetscErrorCode DMPlexDistortGeometry(DM dm)
369 {
370   DM             cdm;
371   DMLabel        label;
372   Vec            coordinates;
373   PetscScalar   *coords;
374   PetscReal      mid = 0.5;
375   PetscInt       cdim, d, vStart, vEnd, v;
376   PetscErrorCode ierr;
377 
378   PetscFunctionBeginUser;
379   ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
380   ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr);
381   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
382   ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr);
383   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
384   ierr = VecGetArrayWrite(coordinates, &coords);CHKERRQ(ierr);
385   for (v = vStart; v < vEnd; ++v) {
386     PetscScalar *pcoords, shift;
387     PetscInt     val;
388 
389     ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr);
390     if (val >= 0) continue;
391     ierr = DMPlexPointLocalRef(cdm, v, coords, &pcoords);CHKERRQ(ierr);
392     shift = 0.2 * PetscAbsScalar(pcoords[0] - mid);
393     shift = PetscRealPart(pcoords[0]) > mid ? shift : -shift;
394     for (d = 1; d < cdim; ++d) pcoords[d] += shift;
395   }
396   ierr = VecRestoreArrayWrite(coordinates, &coords);CHKERRQ(ierr);
397   PetscFunctionReturn(0);
398 }
399 
400 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
401 {
402   PetscErrorCode ierr;
403 
404   PetscFunctionBeginUser;
405   ierr = DMCreate(comm, dm);CHKERRQ(ierr);
406   ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
407   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
408   switch (user->deform) {
409     case DEFORM_NONE:  break;
410     case DEFORM_SHEAR: ierr = DMPlexShearGeometry(*dm, DM_X, NULL);CHKERRQ(ierr);break;
411     case DEFORM_STEP:  ierr = DMPlexDistortGeometry(*dm);CHKERRQ(ierr);break;
412     default: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid deformation type: %s (%D)", deformTypes[PetscMin(user->deform, NUM_DEFORM_TYPES)], user->deform);
413   }
414   ierr = DMSetApplicationContext(*dm, user);CHKERRQ(ierr);
415   ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
416   PetscFunctionReturn(0);
417 }
418 
419 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
420 {
421   PetscErrorCode (*exact)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
422   PetscDS        ds;
423   PetscWeakForm  wf;
424   DMLabel        label;
425   PetscInt       id, bd;
426   PetscInt       dim;
427   PetscErrorCode ierr;
428 
429   PetscFunctionBeginUser;
430   ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);
431   ierr = PetscDSGetWeakForm(ds, &wf);CHKERRQ(ierr);
432   ierr = PetscDSGetSpatialDimension(ds, &dim);CHKERRQ(ierr);
433   switch (user->solType) {
434   case SOL_MASS_QUADRATIC:
435     ierr = PetscDSSetResidual(ds, 0, f0_mass_u, NULL);CHKERRQ(ierr);
436     ierr = PetscDSSetJacobian(ds, 0, 0, g0_mass_uu, NULL, NULL, NULL);CHKERRQ(ierr);
437     ierr = PetscWeakFormSetIndexResidual(wf, NULL, 0, 0,  0, 1, f0_mass_quadratic_u, 0, NULL);CHKERRQ(ierr);
438     switch (dim) {
439     case 2: exact = quadratic_2d_u;break;
440     case 3: exact = quadratic_3d_u;break;
441     default: SETERRQ1(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
442     }
443     break;
444   case SOL_VLAP_QUADRATIC:
445     ierr = PetscDSSetResidual(ds, 0, f0_vlap_quadratic_u, f1_vlap_u);CHKERRQ(ierr);
446     ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_vlap_uu);CHKERRQ(ierr);
447     switch (dim) {
448     case 2: exact = quadratic_2d_u;break;
449     case 3: exact = quadratic_3d_u;break;
450     default: SETERRQ1(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
451     }
452     break;
453   case SOL_ELAS_QUADRATIC:
454     ierr = PetscDSSetResidual(ds, 0, f0_elas_quadratic_u, f1_elas_u);CHKERRQ(ierr);
455     ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr);
456     switch (dim) {
457     case 2: exact = quadratic_2d_u;break;
458     case 3: exact = quadratic_3d_u;break;
459     default: SETERRQ1(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
460     }
461     break;
462   case SOL_VLAP_TRIG:
463     ierr = PetscDSSetResidual(ds, 0, f0_vlap_trig_u, f1_vlap_u);CHKERRQ(ierr);
464     ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_vlap_uu);CHKERRQ(ierr);
465     switch (dim) {
466     case 2: exact = trig_2d_u;break;
467     case 3: exact = trig_3d_u;break;
468     default: SETERRQ1(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
469     }
470     break;
471   case SOL_ELAS_TRIG:
472     ierr = PetscDSSetResidual(ds, 0, f0_elas_trig_u, f1_elas_u);CHKERRQ(ierr);
473     ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr);
474     switch (dim) {
475     case 2: exact = trig_2d_u;break;
476     case 3: exact = trig_3d_u;break;
477     default: SETERRQ1(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
478     }
479     break;
480   case SOL_ELAS_AXIAL_DISP:
481     ierr = PetscDSSetResidual(ds, 0, NULL, f1_elas_u);CHKERRQ(ierr);
482     ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr);
483     id   = dim == 3 ? 5 : 2;
484     ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr);
485     ierr = DMAddBoundary(dm, DM_BC_NATURAL, "right", label, 1, &id, 0, 0, NULL, (void (*)(void)) NULL, NULL, user, &bd);CHKERRQ(ierr);
486     ierr = PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
487     ierr = PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_elas_axial_disp_bd_u, 0, NULL);CHKERRQ(ierr);
488     exact = axial_disp_u;
489     break;
490   case SOL_ELAS_UNIFORM_STRAIN:
491     ierr = PetscDSSetResidual(ds, 0, NULL, f1_elas_u);CHKERRQ(ierr);
492     ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr);
493     exact = uniform_strain_u;
494     break;
495   default: SETERRQ2(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType);
496   }
497   ierr = PetscDSSetExactSolution(ds, 0, exact, user);CHKERRQ(ierr);
498   ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr);
499   if (user->solType == SOL_ELAS_AXIAL_DISP) {
500     PetscInt cmp;
501 
502     id   = dim == 3 ? 6 : 4;
503     cmp  = 0;
504     ierr = DMAddBoundary(dm,   DM_BC_ESSENTIAL, "left",   label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);CHKERRQ(ierr);
505     cmp  = dim == 3 ? 2 : 1;
506     id   = dim == 3 ? 1 : 1;
507     ierr = DMAddBoundary(dm,   DM_BC_ESSENTIAL, "bottom", label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);CHKERRQ(ierr);
508     if (dim == 3) {
509       cmp  = 1;
510       id   = 3;
511       ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "front",  label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);CHKERRQ(ierr);
512     }
513   } else {
514     id = 1;
515     ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) exact, NULL, user, NULL);CHKERRQ(ierr);
516   }
517   PetscFunctionReturn(0);
518 }
519 
520 static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
521 {
522   PetscErrorCode ierr;
523 
524   PetscFunctionBegin;
525   ierr = DMPlexCreateRigidBody(dm, origField, nullspace);CHKERRQ(ierr);
526   PetscFunctionReturn(0);
527 }
528 
529 PetscErrorCode SetupFE(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx)
530 {
531   AppCtx        *user = (AppCtx *) ctx;
532   DM             cdm  = dm;
533   PetscFE        fe;
534   char           prefix[PETSC_MAX_PATH_LEN];
535   DMPolytopeType ct;
536   PetscBool      simplex;
537   PetscInt       dim, cStart;
538   PetscErrorCode ierr;
539 
540   PetscFunctionBegin;
541   /* Create finite element */
542   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
543   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr);
544   ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr);
545   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
546   ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr);
547   ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, dim, simplex, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr);
548   ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr);
549   /* Set discretization and boundary conditions for each mesh */
550   ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr);
551   ierr = DMCreateDS(dm);CHKERRQ(ierr);
552   ierr = (*setup)(dm, user);CHKERRQ(ierr);
553   while (cdm) {
554     ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr);
555     if (user->useNearNullspace) {ierr = DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace);CHKERRQ(ierr);}
556     /* TODO: Check whether the boundary of coarse meshes is marked */
557     ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr);
558   }
559   ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
560   PetscFunctionReturn(0);
561 }
562 
563 int main(int argc, char **argv)
564 {
565   DM             dm;   /* Problem specification */
566   SNES           snes; /* Nonlinear solver */
567   Vec            u;    /* Solutions */
568   AppCtx         user; /* User-defined work context */
569   PetscErrorCode ierr;
570 
571   ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
572   ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
573   /* Primal system */
574   ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);
575   ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
576   ierr = SNESSetDM(snes, dm);CHKERRQ(ierr);
577   ierr = SetupFE(dm, "displacement", SetupPrimalProblem, &user);CHKERRQ(ierr);
578   ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
579   ierr = VecSet(u, 0.0);CHKERRQ(ierr);
580   ierr = PetscObjectSetName((PetscObject) u, "displacement");CHKERRQ(ierr);
581   ierr = DMPlexSetSNESLocalFEM(dm, &user, &user, &user);CHKERRQ(ierr);
582   ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
583   ierr = DMSNESCheckFromOptions(snes, u);CHKERRQ(ierr);
584   ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr);
585   ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr);
586   ierr = VecViewFromOptions(u, NULL, "-displacement_view");CHKERRQ(ierr);
587   /* Cleanup */
588   ierr = VecDestroy(&u);CHKERRQ(ierr);
589   ierr = SNESDestroy(&snes);CHKERRQ(ierr);
590   ierr = DMDestroy(&dm);CHKERRQ(ierr);
591   ierr = PetscFinalize();
592   return ierr;
593 }
594 
595 /*TEST
596 
597   testset:
598     args: -dm_plex_box_faces 1,1,1
599 
600     test:
601       suffix: 2d_p1_quad_vlap
602       requires: triangle
603       args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
604     test:
605       suffix: 2d_p2_quad_vlap
606       requires: triangle
607       args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
608     test:
609       suffix: 2d_p3_quad_vlap
610       requires: triangle
611       args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
612     test:
613       suffix: 2d_q1_quad_vlap
614       args: -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
615     test:
616       suffix: 2d_q2_quad_vlap
617       args: -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
618     test:
619       suffix: 2d_q3_quad_vlap
620       requires: !single
621       args: -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
622     test:
623       suffix: 2d_p1_quad_elas
624       requires: triangle
625       args: -sol_type elas_quad -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
626     test:
627       suffix: 2d_p2_quad_elas
628       requires: triangle
629       args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001
630     test:
631       suffix: 2d_p3_quad_elas
632       requires: triangle
633       args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001
634     test:
635       suffix: 2d_q1_quad_elas
636       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
637     test:
638       suffix: 2d_q1_quad_elas_shear
639       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
640     test:
641       suffix: 2d_q2_quad_elas
642       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001
643     test:
644       suffix: 2d_q2_quad_elas_shear
645       args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 2 -dmsnes_check
646     test:
647       suffix: 2d_q3_quad_elas
648       args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001
649     test:
650       suffix: 2d_q3_quad_elas_shear
651       requires: !single
652       args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 3 -dmsnes_check
653 
654     test:
655       suffix: 3d_p1_quad_vlap
656       requires: ctetgen
657       args: -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
658     test:
659       suffix: 3d_p2_quad_vlap
660       requires: ctetgen
661       args: -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
662     test:
663       suffix: 3d_p3_quad_vlap
664       requires: ctetgen
665       args: -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
666     test:
667       suffix: 3d_q1_quad_vlap
668       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
669     test:
670       suffix: 3d_q2_quad_vlap
671       args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
672     test:
673       suffix: 3d_q3_quad_vlap
674       args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
675     test:
676       suffix: 3d_p1_quad_elas
677       requires: ctetgen
678       args: -sol_type elas_quad -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
679     test:
680       suffix: 3d_p2_quad_elas
681       requires: ctetgen
682       args: -sol_type elas_quad -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
683     test:
684       suffix: 3d_p3_quad_elas
685       requires: ctetgen
686       args: -sol_type elas_quad -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
687     test:
688       suffix: 3d_q1_quad_elas
689       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
690     test:
691       suffix: 3d_q2_quad_elas
692       args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
693     test:
694       suffix: 3d_q3_quad_elas
695       requires: !single
696       args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
697 
698     test:
699       suffix: 2d_p1_trig_vlap
700       requires: triangle
701       args: -sol_type vlap_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
702     test:
703       suffix: 2d_p2_trig_vlap
704       requires: triangle
705       args: -sol_type vlap_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
706     test:
707       suffix: 2d_p3_trig_vlap
708       requires: triangle
709       args: -sol_type vlap_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
710     test:
711       suffix: 2d_q1_trig_vlap
712       args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
713     test:
714       suffix: 2d_q2_trig_vlap
715       args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
716     test:
717       suffix: 2d_q3_trig_vlap
718       args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
719     test:
720       suffix: 2d_p1_trig_elas
721       requires: triangle
722       args: -sol_type elas_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
723     test:
724       suffix: 2d_p2_trig_elas
725       requires: triangle
726       args: -sol_type elas_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
727     test:
728       suffix: 2d_p3_trig_elas
729       requires: triangle
730       args: -sol_type elas_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
731     test:
732       suffix: 2d_q1_trig_elas
733       args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
734     test:
735       suffix: 2d_q1_trig_elas_shear
736       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
737     test:
738       suffix: 2d_q2_trig_elas
739       args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
740     test:
741       suffix: 2d_q2_trig_elas_shear
742       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
743     test:
744       suffix: 2d_q3_trig_elas
745       args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
746     test:
747       suffix: 2d_q3_trig_elas_shear
748       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
749 
750     test:
751       suffix: 3d_p1_trig_vlap
752       requires: ctetgen
753       args: -sol_type vlap_trig -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
754     test:
755       suffix: 3d_p2_trig_vlap
756       requires: ctetgen
757       args: -sol_type vlap_trig -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
758     test:
759       suffix: 3d_p3_trig_vlap
760       requires: ctetgen
761       args: -sol_type vlap_trig -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
762     test:
763       suffix: 3d_q1_trig_vlap
764       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
765     test:
766       suffix: 3d_q2_trig_vlap
767       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
768     test:
769       suffix: 3d_q3_trig_vlap
770       requires: !__float128
771       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
772     test:
773       suffix: 3d_p1_trig_elas
774       requires: ctetgen
775       args: -sol_type elas_trig -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
776     test:
777       suffix: 3d_p2_trig_elas
778       requires: ctetgen
779       args: -sol_type elas_trig -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
780     test:
781       suffix: 3d_p3_trig_elas
782       requires: ctetgen
783       args: -sol_type elas_trig -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
784     test:
785       suffix: 3d_q1_trig_elas
786       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
787     test:
788       suffix: 3d_q2_trig_elas
789       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
790     test:
791       suffix: 3d_q3_trig_elas
792       requires: !__float128
793       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
794 
795     test:
796       suffix: 2d_p1_axial_elas
797       requires: triangle
798       args: -sol_type elas_axial_disp -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 2 -dmsnes_check .0001 -pc_type lu
799     test:
800       suffix: 2d_p2_axial_elas
801       requires: triangle
802       args: -sol_type elas_axial_disp -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
803     test:
804       suffix: 2d_p3_axial_elas
805       requires: triangle
806       args: -sol_type elas_axial_disp -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
807     test:
808       suffix: 2d_q1_axial_elas
809       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
810     test:
811       suffix: 2d_q2_axial_elas
812       args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type   lu
813     test:
814       suffix: 2d_q3_axial_elas
815       args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type   lu
816 
817     test:
818       suffix: 2d_p1_uniform_elas
819       requires: triangle
820       args: -sol_type elas_uniform_strain -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
821     test:
822       suffix: 2d_p2_uniform_elas
823       requires: triangle
824       args: -sol_type elas_uniform_strain -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
825     test:
826       suffix: 2d_p3_uniform_elas
827       requires: triangle
828       args: -sol_type elas_uniform_strain -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
829     test:
830       suffix: 2d_q1_uniform_elas
831       args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
832     test:
833       suffix: 2d_q2_uniform_elas
834       requires: !single
835       args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
836     test:
837       suffix: 2d_q3_uniform_elas
838       requires: !single
839       args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
840     test:
841       suffix: 2d_p1_uniform_elas_step
842       requires: triangle
843       args: -sol_type elas_uniform_strain -deform_type step -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
844 
845   testset:
846     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -deform_type step -displacement_petscspace_degree 1 -dmsnes_check .0001 -pc_type lu
847 
848     test:
849       suffix: 2d_q1_uniform_elas_step
850       args: -sol_type elas_uniform_strain -dm_refine 2
851     test:
852       suffix: 2d_q1_quad_vlap_step
853       args:
854     test:
855       suffix: 2d_q2_quad_vlap_step
856       args: -displacement_petscspace_degree 2
857     test:
858       suffix: 2d_q1_quad_mass_step
859       args: -sol_type mass_quad
860 
861 TEST*/
862