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 **) ¶m);
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*/