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\

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>

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;
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) {
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;
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:
537:     requires: triangle
538:     args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
539:   test:
541:     requires: triangle
542:     args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
543:   test:
545:     requires: triangle
546:     args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
547:   test:
549:     args: -simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
550:   test:
552:     args: -simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
553:   test:
555:     requires: !single
556:     args: -simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
557:   test:
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:
563:     requires: triangle
564:     args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001
565:   test:
567:     requires: triangle
568:     args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001
569:   test:
571:     args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
572:   test:
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:
577:     args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001
578:   test:
580:     args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 2 -dmsnes_check
581:   test:
583:     args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001
584:   test:
586:     requires: !single
587:     args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 3 -dmsnes_check

589:   test:
591:     requires: ctetgen
592:     args: -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
593:   test:
595:     requires: ctetgen
596:     args: -dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
597:   test:
599:     requires: ctetgen
600:     args: -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
601:   test:
603:     args: -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
604:   test:
606:     args: -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
607:   test:
609:     args: -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
610:   test:
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:
616:     requires: ctetgen
617:     args: -sol_type elas_quad -dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
618:   test:
620:     requires: ctetgen
621:     args: -sol_type elas_quad -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
622:   test:
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:
627:     args: -sol_type elas_quad -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
628:   test:
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*/