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
9: */
11: #include <petscdmplex.h>
12: #include <petscsnes.h>
13: #include <petscds.h>
14: #include <petscconvest.h>
16: typedef enum {SOL_VLAP_QUADRATIC, SOL_ELAS_QUADRATIC, SOL_VLAP_TRIG, SOL_ELAS_TRIG, SOL_ELAS_AXIAL_DISP, SOL_ELAS_UNIFORM_STRAIN, 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", "unknown"};
19: typedef struct {
20: /* Domain and mesh definition */
21: char dmType[256]; /* DM type for the solve */
22: PetscInt dim; /* The topological mesh dimension */
23: PetscBool simplex; /* Simplicial mesh */
24: PetscInt cells[3]; /* The initial domain division */
25: PetscBool shear; /* Shear the domain */
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;
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: }
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: }
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: }
54: /*
55: u = x^2
56: v = y^2 - 2xy
57: Delta <u,v> - f = <2, 2> - <2, 2>
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: }
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 >
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;
103: for (d = 0; d < dim-1; ++d) f0[d] += 2.0*mu;
104: f0[dim-1] += 2.0*lambda + 4.0*mu;
105: }
107: static PetscErrorCode trig_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
108: {
109: u[0] = PetscSinReal(2.0*PETSC_PI*x[0]);
110: u[1] = PetscSinReal(2.0*PETSC_PI*x[1]) - 2.0*x[0]*x[1];
111: return 0;
112: }
114: static PetscErrorCode trig_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
115: {
116: u[0] = PetscSinReal(2.0*PETSC_PI*x[0]);
117: u[1] = PetscSinReal(2.0*PETSC_PI*x[1]) - 2.0*x[0]*x[1];
118: u[2] = PetscSinReal(2.0*PETSC_PI*x[2]) - 2.0*x[1]*x[2];
119: return 0;
120: }
122: /*
123: u = sin(2 pi x)
124: v = sin(2 pi y) - 2xy
125: 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)>
127: u = sin(2 pi x)
128: v = sin(2 pi y) - 2xy
129: w = sin(2 pi z) - 2yz
130: 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)>
131: */
132: static void f0_vlap_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
133: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
134: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
135: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
136: {
137: PetscInt d;
138: for (d = 0; d < dim; ++d) f0[d] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]);
139: }
141: /*
142: u = sin(2 pi x)
143: v = sin(2 pi y) - 2xy
144: \varepsilon = / 2 pi cos(2 pi x) -y \
145: \ -y 2 pi cos(2 pi y) - 2x /
146: Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x
147: div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
148: = \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) >
149: = \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) >
151: u = sin(2 pi x)
152: v = sin(2 pi y) - 2xy
153: w = sin(2 pi z) - 2yz
154: \varepsilon = / 2 pi cos(2 pi x) -y 0 \
155: | -y 2 pi cos(2 pi y) - 2x -z |
156: \ 0 -z 2 pi cos(2 pi z) - 2y /
157: Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y
158: div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
159: = \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) >
160: = \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) >
161: */
162: static void f0_elas_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
163: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
164: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
165: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
166: {
167: const PetscReal mu = 1.0;
168: const PetscReal lambda = 1.0;
169: const PetscReal fact = 4.0*PetscSqr(PETSC_PI);
170: PetscInt d;
172: 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);
173: }
175: static PetscErrorCode axial_disp_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
176: {
177: const PetscReal mu = 1.0;
178: const PetscReal lambda = 1.0;
179: const PetscReal N = 1.0;
180: PetscInt d;
182: 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];
183: u[1] = -0.25*lambda/mu/(lambda+mu)*N*x[1];
184: for (d = 2; d < dim; ++d) u[d] = 0.0;
185: return 0;
186: }
188: /*
189: We will pull/push on the right side of a block of linearly elastic material. The uniform traction conditions on the
190: right side of the box will result in a uniform strain along x and y. The Neumann BC is given by
192: n_i \sigma_{ij} = t_i
194: u = (1/(2\mu) - 1) x
195: v = -y
196: f = 0
197: t = <4\mu/\lambda (\lambda + \mu), 0>
198: \varepsilon = / 1/(2\mu) - 1 0 \
199: \ 0 -1 /
200: Tr(\varepsilon) = div u = 1/(2\mu) - 2
201: div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
202: = \lambda \partial_j (1/(2\mu) - 2) + 2\mu < 0, 0 >
203: = \lambda < 0, 0 > + \mu < 0, 0 > = 0
204: NBC = <1,0> . <4\mu/\lambda (\lambda + \mu), 0> = 4\mu/\lambda (\lambda + \mu)
206: u = x - 1/2
207: v = 0
208: w = 0
209: \varepsilon = / x 0 0 \
210: | 0 0 0 |
211: \ 0 0 0 /
212: Tr(\varepsilon) = div u = x
213: div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij}
214: = \lambda \partial_j x + 2\mu < 1, 0, 0 >
215: = \lambda < 1, 0, 0 > + \mu < 2, 0, 0 >
216: */
217: static void f0_elas_axial_disp_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
218: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
219: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
220: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
221: {
222: const PetscReal N = -1.0;
224: f0[0] = N;
225: }
227: static PetscErrorCode uniform_strain_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
228: {
229: const PetscReal eps_xx = 0.1;
230: const PetscReal eps_xy = 0.3;
231: const PetscReal eps_yy = 0.25;
232: PetscInt d;
234: u[0] = eps_xx*x[0] + eps_xy*x[1];
235: u[1] = eps_xy*x[0] + eps_yy*x[1];
236: for (d = 2; d < dim; ++d) u[d] = 0.0;
237: return 0;
238: }
240: static void f1_vlap_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
241: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
242: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
243: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
244: {
245: const PetscInt Nc = dim;
246: PetscInt c, d;
248: for (c = 0; c < Nc; ++c) for (d = 0; d < dim; ++d) f1[c*dim+d] += u_x[c*dim+d];
249: }
251: static void f1_elas_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[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
255: {
256: const PetscInt Nc = dim;
257: const PetscReal mu = 1.0;
258: const PetscReal lambda = 1.0;
259: PetscInt c, d;
261: for (c = 0; c < Nc; ++c) {
262: for (d = 0; d < dim; ++d) {
263: f1[c*dim+d] += mu*(u_x[c*dim+d] + u_x[d*dim+c]);
264: f1[c*dim+c] += lambda*u_x[d*dim+d];
265: }
266: }
267: }
269: static void g3_vlap_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
270: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
271: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
272: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
273: {
274: const PetscInt Nc = dim;
275: PetscInt c, d;
277: for (c = 0; c < Nc; ++c) {
278: for (d = 0; d < dim; ++d) {
279: g3[((c*Nc + c)*dim + d)*dim + d] = 1.0;
280: }
281: }
282: }
284: /*
285: \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc
287: \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg}
288: = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc
289: */
290: static void g3_elas_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
291: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
292: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
293: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
294: {
295: const PetscInt Nc = dim;
296: const PetscReal mu = 1.0;
297: const PetscReal lambda = 1.0;
298: PetscInt c, d;
300: for (c = 0; c < Nc; ++c) {
301: for (d = 0; d < dim; ++d) {
302: g3[((c*Nc + c)*dim + d)*dim + d] += mu;
303: g3[((c*Nc + d)*dim + d)*dim + c] += mu;
304: g3[((c*Nc + d)*dim + c)*dim + d] += lambda;
305: }
306: }
307: }
309: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
310: {
311: PetscInt n = 3, sol;
315: options->dim = 2;
316: options->cells[0] = 1;
317: options->cells[1] = 1;
318: options->cells[2] = 1;
319: options->simplex = PETSC_TRUE;
320: options->shear = PETSC_FALSE;
321: options->solType = SOL_VLAP_QUADRATIC;
322: options->useNearNullspace = PETSC_TRUE;
323: PetscStrncpy(options->dmType, DMPLEX, 256);
325: PetscOptionsBegin(comm, "", "Linear Elasticity Problem Options", "DMPLEX");
326: PetscOptionsInt("-dim", "The topological mesh dimension", "ex17.c", options->dim, &options->dim, NULL);
327: PetscOptionsIntArray("-cells", "The initial mesh division", "ex17.c", options->cells, &n, NULL);
328: PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex17.c", options->simplex, &options->simplex, NULL);
329: PetscOptionsBool("-shear", "Shear the domain", "ex17.c", options->shear, &options->shear, NULL);
330: sol = options->solType;
331: PetscOptionsEList("-sol_type", "Type of exact solution", "ex17.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL);
332: options->solType = (SolutionType) sol;
333: PetscOptionsBool("-near_nullspace", "Use the rigid body modes as an AMG near nullspace", "ex17.c", options->useNearNullspace, &options->useNearNullspace, NULL);
334: PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex17.c", DMList, options->dmType, options->dmType, 256, NULL);
335: PetscOptionsEnd();
336: return(0);
337: }
339: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
340: {
341: PetscBool flg;
345: DMPlexCreateBoxMesh(comm, user->dim, user->simplex, user->cells, NULL, NULL, NULL, PETSC_TRUE, dm);
346: {
347: DM pdm = NULL;
348: PetscPartitioner part;
350: DMPlexGetPartitioner(*dm, &part);
351: PetscPartitionerSetFromOptions(part);
352: DMPlexDistribute(*dm, 0, NULL, &pdm);
353: if (pdm) {
354: DMDestroy(dm);
355: *dm = pdm;
356: }
357: }
358: PetscStrcmp(user->dmType, DMPLEX, &flg);
359: if (flg) {
360: DM ndm;
362: DMConvert(*dm, user->dmType, &ndm);
363: if (ndm) {
364: DMDestroy(dm);
365: *dm = ndm;
366: }
367: }
368: if (user->shear) {DMPlexShearGeometry(*dm, DM_X, NULL);}
369: DMLocalizeCoordinates(*dm);
371: PetscObjectSetName((PetscObject) *dm, "Mesh");
372: DMSetApplicationContext(*dm, user);
373: DMSetFromOptions(*dm);
374: DMViewFromOptions(*dm, NULL, "-dm_view");
375: return(0);
376: }
378: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
379: {
380: PetscErrorCode (*exact)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
381: PetscDS prob;
382: PetscInt id;
383: PetscInt dim;
387: DMGetDS(dm, &prob);
388: PetscDSGetSpatialDimension(prob, &dim);
389: switch (user->solType) {
390: case SOL_VLAP_QUADRATIC:
391: PetscDSSetResidual(prob, 0, f0_vlap_quadratic_u, f1_vlap_u);
392: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_vlap_uu);
393: switch (dim) {
394: case 2: exact = quadratic_2d_u;break;
395: case 3: exact = quadratic_3d_u;break;
396: default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
397: }
398: break;
399: case SOL_ELAS_QUADRATIC:
400: PetscDSSetResidual(prob, 0, f0_elas_quadratic_u, f1_elas_u);
401: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
402: switch (dim) {
403: case 2: exact = quadratic_2d_u;break;
404: case 3: exact = quadratic_3d_u;break;
405: default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
406: }
407: break;
408: case SOL_VLAP_TRIG:
409: PetscDSSetResidual(prob, 0, f0_vlap_trig_u, f1_vlap_u);
410: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_vlap_uu);
411: switch (dim) {
412: case 2: exact = trig_2d_u;break;
413: case 3: exact = trig_3d_u;break;
414: default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
415: }
416: break;
417: case SOL_ELAS_TRIG:
418: PetscDSSetResidual(prob, 0, f0_elas_trig_u, f1_elas_u);
419: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
420: switch (dim) {
421: case 2: exact = trig_2d_u;break;
422: case 3: exact = trig_3d_u;break;
423: default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
424: }
425: break;
426: case SOL_ELAS_AXIAL_DISP:
427: PetscDSSetResidual(prob, 0, NULL, f1_elas_u);
428: PetscDSSetBdResidual(prob, 0, f0_elas_axial_disp_bd_u, NULL);
429: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
430: exact = axial_disp_u;
431: break;
432: case SOL_ELAS_UNIFORM_STRAIN:
433: PetscDSSetResidual(prob, 0, NULL, f1_elas_u);
434: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
435: exact = uniform_strain_u;
436: break;
437: default: SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType);
438: }
439: PetscDSSetExactSolution(prob, 0, exact, user);
440: if (user->solType == SOL_ELAS_AXIAL_DISP) {
441: PetscInt cmp;
443: id = dim == 3 ? 5 : 2;
444: DMAddBoundary(dm, DM_BC_NATURAL, "right", "marker", 0, 0, NULL, (void (*)(void)) zero, NULL, 1, &id, user);
445: id = dim == 3 ? 6 : 4;
446: cmp = 0;
447: DMAddBoundary(dm, DM_BC_ESSENTIAL, "left", "marker", 0, 1, &cmp, (void (*)(void)) zero, NULL, 1, &id, user);
448: cmp = dim == 3 ? 2 : 1;
449: id = dim == 3 ? 1 : 1;
450: DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom", "marker", 0, 1, &cmp, (void (*)(void)) zero, NULL, 1, &id, user);
451: if (dim == 3) {
452: cmp = 1;
453: id = 3;
454: DMAddBoundary(dm, DM_BC_ESSENTIAL, "front", "marker", 0, 1, &cmp, (void (*)(void)) zero, NULL, 1, &id, user);
455: }
456: } else {
457: id = 1;
458: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) exact, NULL, 1, &id, user);
459: }
460: return(0);
461: }
463: static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
464: {
468: DMPlexCreateRigidBody(dm, origField, nullspace);
469: return(0);
470: }
472: PetscErrorCode SetupFE(DM dm, PetscInt Nc, PetscBool simplex, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx)
473: {
474: AppCtx *user = (AppCtx *) ctx;
475: DM cdm = dm;
476: PetscFE fe;
477: char prefix[PETSC_MAX_PATH_LEN];
478: PetscInt dim;
482: /* Create finite element */
483: DMGetDimension(dm, &dim);
484: PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);
485: PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, Nc, simplex, name ? prefix : NULL, -1, &fe);
486: PetscObjectSetName((PetscObject) fe, name);
487: /* Set discretization and boundary conditions for each mesh */
488: DMSetField(dm, 0, NULL, (PetscObject) fe);
489: DMCreateDS(dm);
490: (*setup)(dm, user);
491: while (cdm) {
492: DMCopyDisc(dm, cdm);
493: if (user->useNearNullspace) {DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace);}
494: /* TODO: Check whether the boundary of coarse meshes is marked */
495: DMGetCoarseDM(cdm, &cdm);
496: }
497: PetscFEDestroy(&fe);
498: return(0);
499: }
501: int main(int argc, char **argv)
502: {
503: DM dm; /* Problem specification */
504: SNES snes; /* Nonlinear solver */
505: Vec u; /* Solutions */
506: AppCtx user; /* User-defined work context */
509: PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
510: ProcessOptions(PETSC_COMM_WORLD, &user);
511: /* Primal system */
512: SNESCreate(PETSC_COMM_WORLD, &snes);
513: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
514: SNESSetDM(snes, dm);
515: SetupFE(dm, user.dim, user.simplex, "displacement", SetupPrimalProblem, &user);
516: DMCreateGlobalVector(dm, &u);
517: VecSet(u, 0.0);
518: PetscObjectSetName((PetscObject) u, "displacement");
519: DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
520: SNESSetFromOptions(snes);
521: DMSNESCheckFromOptions(snes, u);
522: SNESSolve(snes, NULL, u);
523: SNESGetSolution(snes, &u);
524: VecViewFromOptions(u, NULL, "-displacement_view");
525: /* Cleanup */
526: VecDestroy(&u);
527: SNESDestroy(&snes);
528: DMDestroy(&dm);
529: PetscFinalize();
530: return ierr;
531: }
533: /*TEST
535: test:
536: suffix: 2d_p1_quad_vlap
537: requires: triangle
538: args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
539: test:
540: suffix: 2d_p2_quad_vlap
541: requires: triangle
542: args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
543: test:
544: suffix: 2d_p3_quad_vlap
545: requires: triangle
546: args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
547: test:
548: suffix: 2d_q1_quad_vlap
549: args: -simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
550: test:
551: suffix: 2d_q2_quad_vlap
552: args: -simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
553: test:
554: suffix: 2d_q3_quad_vlap
555: requires: !single
556: args: -simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
557: test:
558: suffix: 2d_p1_quad_elas
559: requires: triangle
560: args: -sol_type elas_quad -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
561: test:
562: suffix: 2d_p2_quad_elas
563: requires: triangle
564: args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001
565: test:
566: suffix: 2d_p3_quad_elas
567: requires: triangle
568: args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001
569: test:
570: suffix: 2d_q1_quad_elas
571: args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
572: test:
573: suffix: 2d_q1_quad_elas_shear
574: args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
575: test:
576: suffix: 2d_q2_quad_elas
577: args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001
578: test:
579: suffix: 2d_q2_quad_elas_shear
580: args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 2 -dmsnes_check
581: test:
582: suffix: 2d_q3_quad_elas
583: args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001
584: test:
585: suffix: 2d_q3_quad_elas_shear
586: requires: !single
587: args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 3 -dmsnes_check
589: test:
590: suffix: 3d_p1_quad_vlap
591: requires: ctetgen
592: args: -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
593: test:
594: suffix: 3d_p2_quad_vlap
595: requires: ctetgen
596: args: -dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
597: test:
598: suffix: 3d_p3_quad_vlap
599: requires: ctetgen
600: args: -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
601: test:
602: suffix: 3d_q1_quad_vlap
603: args: -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
604: test:
605: suffix: 3d_q2_quad_vlap
606: args: -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
607: test:
608: suffix: 3d_q3_quad_vlap
609: args: -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
610: test:
611: suffix: 3d_p1_quad_elas
612: requires: ctetgen
613: args: -sol_type elas_quad -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
614: test:
615: suffix: 3d_p2_quad_elas
616: requires: ctetgen
617: args: -sol_type elas_quad -dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
618: test:
619: suffix: 3d_p3_quad_elas
620: requires: ctetgen
621: args: -sol_type elas_quad -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
622: test:
623: suffix: 3d_q1_quad_elas
624: args: -sol_type elas_quad -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
625: test:
626: suffix: 3d_q2_quad_elas
627: args: -sol_type elas_quad -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
628: test:
629: suffix: 3d_q3_quad_elas
630: requires: !single
631: args: -sol_type elas_quad -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
633: test:
634: suffix: 2d_p1_trig_vlap
635: requires: triangle
636: args: -sol_type vlap_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
637: test:
638: suffix: 2d_p2_trig_vlap
639: requires: triangle
640: args: -sol_type vlap_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
641: test:
642: suffix: 2d_p3_trig_vlap
643: requires: triangle
644: args: -sol_type vlap_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
645: test:
646: suffix: 2d_q1_trig_vlap
647: args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
648: test:
649: suffix: 2d_q2_trig_vlap
650: args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
651: test:
652: suffix: 2d_q3_trig_vlap
653: args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
654: test:
655: suffix: 2d_p1_trig_elas
656: requires: triangle
657: args: -sol_type elas_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
658: test:
659: suffix: 2d_p2_trig_elas
660: requires: triangle
661: args: -sol_type elas_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
662: test:
663: suffix: 2d_p3_trig_elas
664: requires: triangle
665: args: -sol_type elas_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
666: test:
667: suffix: 2d_q1_trig_elas
668: args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
669: test:
670: suffix: 2d_q1_trig_elas_shear
671: args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
672: test:
673: suffix: 2d_q2_trig_elas
674: args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
675: test:
676: suffix: 2d_q2_trig_elas_shear
677: args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
678: test:
679: suffix: 2d_q3_trig_elas
680: args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
681: test:
682: suffix: 2d_q3_trig_elas_shear
683: args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
685: test:
686: suffix: 3d_p1_trig_vlap
687: requires: ctetgen
688: args: -sol_type vlap_trig -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
689: test:
690: suffix: 3d_p2_trig_vlap
691: requires: ctetgen
692: args: -sol_type vlap_trig -dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
693: test:
694: suffix: 3d_p3_trig_vlap
695: requires: ctetgen
696: args: -sol_type vlap_trig -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
697: test:
698: suffix: 3d_q1_trig_vlap
699: args: -sol_type vlap_trig -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
700: test:
701: suffix: 3d_q2_trig_vlap
702: args: -sol_type vlap_trig -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
703: test:
704: suffix: 3d_q3_trig_vlap
705: requires: !__float128
706: args: -sol_type vlap_trig -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
707: test:
708: suffix: 3d_p1_trig_elas
709: requires: ctetgen
710: args: -sol_type elas_trig -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
711: test:
712: suffix: 3d_p2_trig_elas
713: requires: ctetgen
714: args: -sol_type elas_trig -dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
715: test:
716: suffix: 3d_p3_trig_elas
717: requires: ctetgen
718: args: -sol_type elas_trig -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
719: test:
720: suffix: 3d_q1_trig_elas
721: args: -sol_type elas_trig -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
722: test:
723: suffix: 3d_q2_trig_elas
724: args: -sol_type elas_trig -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
725: test:
726: suffix: 3d_q3_trig_elas
727: requires: !__float128
728: args: -sol_type elas_trig -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
730: test:
731: suffix: 2d_p1_axial_elas
732: requires: triangle
733: args: -sol_type elas_axial_disp -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 2 -dmsnes_check .0001 -pc_type lu
734: test:
735: suffix: 2d_p2_axial_elas
736: requires: triangle
737: args: -sol_type elas_axial_disp -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
738: test:
739: suffix: 2d_p3_axial_elas
740: requires: triangle
741: args: -sol_type elas_axial_disp -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
742: test:
743: suffix: 2d_q1_axial_elas
744: args: -sol_type elas_axial_disp -simplex 0 -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 1 -dmsnes_check .0001 -pc_type lu
745: test:
746: suffix: 2d_q2_axial_elas
747: args: -sol_type elas_axial_disp -simplex 0 -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
748: test:
749: suffix: 2d_q3_axial_elas
750: args: -sol_type elas_axial_disp -simplex 0 -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
752: test:
753: suffix: 2d_p1_uniform_elas
754: requires: triangle
755: args: -sol_type elas_uniform_strain -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
756: test:
757: suffix: 2d_p2_uniform_elas
758: requires: triangle
759: args: -sol_type elas_uniform_strain -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
760: test:
761: suffix: 2d_p3_uniform_elas
762: requires: triangle
763: args: -sol_type elas_uniform_strain -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
764: test:
765: suffix: 2d_q1_uniform_elas
766: args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
767: test:
768: suffix: 2d_q2_uniform_elas
769: requires: !single
770: args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
771: test:
772: suffix: 2d_q3_uniform_elas
773: requires: !single
774: args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
776: TEST*/