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 {
22: SOL_VLAP_QUADRATIC,
23: SOL_ELAS_QUADRATIC,
24: SOL_VLAP_TRIG,
25: SOL_ELAS_TRIG,
26: SOL_ELAS_AXIAL_DISP,
27: SOL_ELAS_UNIFORM_STRAIN,
28: SOL_ELAS_GE,
29: SOL_MASS_QUADRATIC,
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 **)¶m));
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(PetscOptionsCreateViewer(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(PetscViewerDestroy(&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 **)¶m));
466: switch (user->solType) {
467: case SOL_MASS_QUADRATIC:
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;
482: case SOL_VLAP_QUADRATIC:
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;
496: case SOL_ELAS_QUADRATIC:
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:
679: suffix: 2d_p1_quad_vlap
680: requires: triangle
681: args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
682: test:
683: suffix: 2d_p2_quad_vlap
684: requires: triangle
685: args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
686: test:
687: suffix: 2d_p3_quad_vlap
688: requires: triangle
689: args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
690: test:
691: suffix: 2d_q1_quad_vlap
692: args: -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
693: test:
694: suffix: 2d_q2_quad_vlap
695: args: -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
696: test:
697: suffix: 2d_q3_quad_vlap
698: requires: !single
699: args: -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
700: test:
701: suffix: 2d_p1_quad_elas
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:
705: suffix: 2d_p2_quad_elas
706: requires: triangle
707: args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001
708: test:
709: suffix: 2d_p3_quad_elas
710: requires: triangle
711: args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001
712: test:
713: suffix: 2d_q1_quad_elas
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:
716: suffix: 2d_q1_quad_elas_shear
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:
719: suffix: 2d_q2_quad_elas
720: args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001
721: test:
722: suffix: 2d_q2_quad_elas_shear
723: args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 2 -dmsnes_check
724: test:
725: suffix: 2d_q3_quad_elas
726: args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001
727: test:
728: suffix: 2d_q3_quad_elas_shear
729: requires: !single
730: args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 3 -dmsnes_check
732: test:
733: suffix: 3d_p1_quad_vlap
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:
737: suffix: 3d_p2_quad_vlap
738: requires: ctetgen
739: args: -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
740: test:
741: suffix: 3d_p3_quad_vlap
742: requires: ctetgen
743: args: -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
744: test:
745: suffix: 3d_q1_quad_vlap
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:
748: suffix: 3d_q2_quad_vlap
749: args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
750: test:
751: suffix: 3d_q3_quad_vlap
752: args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
753: test:
754: suffix: 3d_p1_quad_elas
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:
758: suffix: 3d_p2_quad_elas
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:
762: suffix: 3d_p3_quad_elas
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:
766: suffix: 3d_q1_quad_elas
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:
769: suffix: 3d_q2_quad_elas
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:
772: suffix: 3d_q3_quad_elas
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:
930: suffix: 2d_q1_quad_vlap_step
931: args:
932: test:
933: suffix: 2d_q2_quad_vlap_step
934: args: -displacement_petscspace_degree 2
935: test:
936: suffix: 2d_q1_quad_mass_step
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*/