Actual source code: ex17.c

petsc-main 2021-04-20
Report Typos and Errors
  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:   PetscWeakForm  wf;
383:   DMLabel        label;
384:   PetscInt       id, bd;
385:   PetscInt       dim;

389:   DMGetDS(dm, &prob);
390:   PetscDSGetSpatialDimension(prob, &dim);
391:   switch (user->solType) {
392:   case SOL_VLAP_QUADRATIC:
393:     PetscDSSetResidual(prob, 0, f0_vlap_quadratic_u, f1_vlap_u);
394:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_vlap_uu);
395:     switch (dim) {
396:     case 2: exact = quadratic_2d_u;break;
397:     case 3: exact = quadratic_3d_u;break;
398:     default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
399:     }
400:     break;
401:   case SOL_ELAS_QUADRATIC:
402:     PetscDSSetResidual(prob, 0, f0_elas_quadratic_u, f1_elas_u);
403:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
404:     switch (dim) {
405:     case 2: exact = quadratic_2d_u;break;
406:     case 3: exact = quadratic_3d_u;break;
407:     default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
408:     }
409:     break;
410:   case SOL_VLAP_TRIG:
411:     PetscDSSetResidual(prob, 0, f0_vlap_trig_u, f1_vlap_u);
412:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_vlap_uu);
413:     switch (dim) {
414:     case 2: exact = trig_2d_u;break;
415:     case 3: exact = trig_3d_u;break;
416:     default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
417:     }
418:     break;
419:   case SOL_ELAS_TRIG:
420:     PetscDSSetResidual(prob, 0, f0_elas_trig_u, f1_elas_u);
421:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
422:     switch (dim) {
423:     case 2: exact = trig_2d_u;break;
424:     case 3: exact = trig_3d_u;break;
425:     default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim);
426:     }
427:     break;
428:   case SOL_ELAS_AXIAL_DISP:
429:     PetscDSSetResidual(prob, 0, NULL, f1_elas_u);
430:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
431:     id   = dim == 3 ? 5 : 2;
432:     DMGetLabel(dm, "marker", &label);
433:     DMAddBoundary(dm, DM_BC_NATURAL, "right", label, 1, &id, 0, 0, NULL, (void (*)(void)) NULL, NULL, user, &bd);
434:     PetscDSGetBoundary(prob, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
435:     PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, f0_elas_axial_disp_bd_u, 0, NULL);
436:     exact = axial_disp_u;
437:     break;
438:   case SOL_ELAS_UNIFORM_STRAIN:
439:     PetscDSSetResidual(prob, 0, NULL, f1_elas_u);
440:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);
441:     exact = uniform_strain_u;
442:     break;
443:   default: SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType);
444:   }
445:   PetscDSSetExactSolution(prob, 0, exact, user);
446:   DMGetLabel(dm, "marker", &label);
447:   if (user->solType == SOL_ELAS_AXIAL_DISP) {
448:     PetscInt cmp;

450:     id   = dim == 3 ? 6 : 4;
451:     cmp  = 0;
452:     DMAddBoundary(dm,   DM_BC_ESSENTIAL, "left",   label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);
453:     cmp  = dim == 3 ? 2 : 1;
454:     id   = dim == 3 ? 1 : 1;
455:     DMAddBoundary(dm,   DM_BC_ESSENTIAL, "bottom", label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);
456:     if (dim == 3) {
457:       cmp  = 1;
458:       id   = 3;
459:       DMAddBoundary(dm, DM_BC_ESSENTIAL, "front",  label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);
460:     }
461:   } else {
462:     id = 1;
463:     DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) exact, NULL, user, NULL);
464:   }
465:   return(0);
466: }

468: static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace)
469: {

473:   DMPlexCreateRigidBody(dm, origField, nullspace);
474:   return(0);
475: }

477: PetscErrorCode SetupFE(DM dm, PetscInt Nc, PetscBool simplex, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx)
478: {
479:   AppCtx        *user = (AppCtx *) ctx;
480:   DM             cdm  = dm;
481:   PetscFE        fe;
482:   char           prefix[PETSC_MAX_PATH_LEN];
483:   PetscInt       dim;

487:   /* Create finite element */
488:   DMGetDimension(dm, &dim);
489:   PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);
490:   PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, Nc, simplex, name ? prefix : NULL, -1, &fe);
491:   PetscObjectSetName((PetscObject) fe, name);
492:   /* Set discretization and boundary conditions for each mesh */
493:   DMSetField(dm, 0, NULL, (PetscObject) fe);
494:   DMCreateDS(dm);
495:   (*setup)(dm, user);
496:   while (cdm) {
497:     DMCopyDisc(dm, cdm);
498:     if (user->useNearNullspace) {DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace);}
499:     /* TODO: Check whether the boundary of coarse meshes is marked */
500:     DMGetCoarseDM(cdm, &cdm);
501:   }
502:   PetscFEDestroy(&fe);
503:   return(0);
504: }

506: int main(int argc, char **argv)
507: {
508:   DM             dm;   /* Problem specification */
509:   SNES           snes; /* Nonlinear solver */
510:   Vec            u;    /* Solutions */
511:   AppCtx         user; /* User-defined work context */

514:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
515:   ProcessOptions(PETSC_COMM_WORLD, &user);
516:   /* Primal system */
517:   SNESCreate(PETSC_COMM_WORLD, &snes);
518:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
519:   SNESSetDM(snes, dm);
520:   SetupFE(dm, user.dim, user.simplex, "displacement", SetupPrimalProblem, &user);
521:   DMCreateGlobalVector(dm, &u);
522:   VecSet(u, 0.0);
523:   PetscObjectSetName((PetscObject) u, "displacement");
524:   DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
525:   SNESSetFromOptions(snes);
526:   DMSNESCheckFromOptions(snes, u);
527:   SNESSolve(snes, NULL, u);
528:   SNESGetSolution(snes, &u);
529:   VecViewFromOptions(u, NULL, "-displacement_view");
530:   /* Cleanup */
531:   VecDestroy(&u);
532:   SNESDestroy(&snes);
533:   DMDestroy(&dm);
534:   PetscFinalize();
535:   return ierr;
536: }

538: /*TEST

540:   test:
541:     suffix: 2d_p1_quad_vlap
542:     requires: triangle
543:     args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
544:   test:
545:     suffix: 2d_p2_quad_vlap
546:     requires: triangle
547:     args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
548:   test:
549:     suffix: 2d_p3_quad_vlap
550:     requires: triangle
551:     args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
552:   test:
553:     suffix: 2d_q1_quad_vlap
554:     args: -simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
555:   test:
556:     suffix: 2d_q2_quad_vlap
557:     args: -simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001
558:   test:
559:     suffix: 2d_q3_quad_vlap
560:     requires: !single
561:     args: -simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001
562:   test:
563:     suffix: 2d_p1_quad_elas
564:     requires: triangle
565:     args: -sol_type elas_quad -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate
566:   test:
567:     suffix: 2d_p2_quad_elas
568:     requires: triangle
569:     args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001
570:   test:
571:     suffix: 2d_p3_quad_elas
572:     requires: triangle
573:     args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001
574:   test:
575:     suffix: 2d_q1_quad_elas
576:     args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
577:   test:
578:     suffix: 2d_q1_quad_elas_shear
579:     args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
580:   test:
581:     suffix: 2d_q2_quad_elas
582:     args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001
583:   test:
584:     suffix: 2d_q2_quad_elas_shear
585:     args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 2 -dmsnes_check
586:   test:
587:     suffix: 2d_q3_quad_elas
588:     args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001
589:   test:
590:     suffix: 2d_q3_quad_elas_shear
591:     requires: !single
592:     args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 3 -dmsnes_check

594:   test:
595:     suffix: 3d_p1_quad_vlap
596:     requires: ctetgen
597:     args: -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
598:   test:
599:     suffix: 3d_p2_quad_vlap
600:     requires: ctetgen
601:     args: -dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
602:   test:
603:     suffix: 3d_p3_quad_vlap
604:     requires: ctetgen
605:     args: -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
606:   test:
607:     suffix: 3d_q1_quad_vlap
608:     args: -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
609:   test:
610:     suffix: 3d_q2_quad_vlap
611:     args: -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
612:   test:
613:     suffix: 3d_q3_quad_vlap
614:     args: -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
615:   test:
616:     suffix: 3d_p1_quad_elas
617:     requires: ctetgen
618:     args: -sol_type elas_quad -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
619:   test:
620:     suffix: 3d_p2_quad_elas
621:     requires: ctetgen
622:     args: -sol_type elas_quad -dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
623:   test:
624:     suffix: 3d_p3_quad_elas
625:     requires: ctetgen
626:     args: -sol_type elas_quad -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001
627:   test:
628:     suffix: 3d_q1_quad_elas
629:     args: -sol_type elas_quad -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
630:   test:
631:     suffix: 3d_q2_quad_elas
632:     args: -sol_type elas_quad -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001
633:   test:
634:     suffix: 3d_q3_quad_elas
635:     requires: !single
636:     args: -sol_type elas_quad -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001

638:   test:
639:     suffix: 2d_p1_trig_vlap
640:     requires: triangle
641:     args: -sol_type vlap_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
642:   test:
643:     suffix: 2d_p2_trig_vlap
644:     requires: triangle
645:     args: -sol_type vlap_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
646:   test:
647:     suffix: 2d_p3_trig_vlap
648:     requires: triangle
649:     args: -sol_type vlap_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
650:   test:
651:     suffix: 2d_q1_trig_vlap
652:     args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
653:   test:
654:     suffix: 2d_q2_trig_vlap
655:     args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
656:   test:
657:     suffix: 2d_q3_trig_vlap
658:     args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
659:   test:
660:     suffix: 2d_p1_trig_elas
661:     requires: triangle
662:     args: -sol_type elas_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
663:   test:
664:     suffix: 2d_p2_trig_elas
665:     requires: triangle
666:     args: -sol_type elas_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
667:   test:
668:     suffix: 2d_p3_trig_elas
669:     requires: triangle
670:     args: -sol_type elas_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
671:   test:
672:     suffix: 2d_q1_trig_elas
673:     args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
674:   test:
675:     suffix: 2d_q1_trig_elas_shear
676:     args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
677:   test:
678:     suffix: 2d_q2_trig_elas
679:     args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
680:   test:
681:     suffix: 2d_q2_trig_elas_shear
682:     args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
683:   test:
684:     suffix: 2d_q3_trig_elas
685:     args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate
686:   test:
687:     suffix: 2d_q3_trig_elas_shear
688:     args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate

690:   test:
691:     suffix: 3d_p1_trig_vlap
692:     requires: ctetgen
693:     args: -sol_type vlap_trig -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
694:   test:
695:     suffix: 3d_p2_trig_vlap
696:     requires: ctetgen
697:     args: -sol_type vlap_trig -dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
698:   test:
699:     suffix: 3d_p3_trig_vlap
700:     requires: ctetgen
701:     args: -sol_type vlap_trig -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
702:   test:
703:     suffix: 3d_q1_trig_vlap
704:     args: -sol_type vlap_trig -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
705:   test:
706:     suffix: 3d_q2_trig_vlap
707:     args: -sol_type vlap_trig -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
708:   test:
709:     suffix: 3d_q3_trig_vlap
710:     requires: !__float128
711:     args: -sol_type vlap_trig -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
712:   test:
713:     suffix: 3d_p1_trig_elas
714:     requires: ctetgen
715:     args: -sol_type elas_trig -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
716:   test:
717:     suffix: 3d_p2_trig_elas
718:     requires: ctetgen
719:     args: -sol_type elas_trig -dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
720:   test:
721:     suffix: 3d_p3_trig_elas
722:     requires: ctetgen
723:     args: -sol_type elas_trig -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
724:   test:
725:     suffix: 3d_q1_trig_elas
726:     args: -sol_type elas_trig -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate
727:   test:
728:     suffix: 3d_q2_trig_elas
729:     args: -sol_type elas_trig -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate
730:   test:
731:     suffix: 3d_q3_trig_elas
732:     requires: !__float128
733:     args: -sol_type elas_trig -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate

735:   test:
736:     suffix: 2d_p1_axial_elas
737:     requires: triangle
738:     args: -sol_type elas_axial_disp -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 2 -dmsnes_check .0001 -pc_type lu
739:   test:
740:     suffix: 2d_p2_axial_elas
741:     requires: triangle
742:     args: -sol_type elas_axial_disp -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
743:   test:
744:     suffix: 2d_p3_axial_elas
745:     requires: triangle
746:     args: -sol_type elas_axial_disp -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
747:   test:
748:     suffix: 2d_q1_axial_elas
749:     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
750:   test:
751:     suffix: 2d_q2_axial_elas
752:     args: -sol_type elas_axial_disp -simplex 0 -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu
753:   test:
754:     suffix: 2d_q3_axial_elas
755:     args: -sol_type elas_axial_disp -simplex 0 -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu

757:   test:
758:     suffix: 2d_p1_uniform_elas
759:     requires: triangle
760:     args: -sol_type elas_uniform_strain -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
761:   test:
762:     suffix: 2d_p2_uniform_elas
763:     requires: triangle
764:     args: -sol_type elas_uniform_strain -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
765:   test:
766:     suffix: 2d_p3_uniform_elas
767:     requires: triangle
768:     args: -sol_type elas_uniform_strain -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
769:   test:
770:     suffix: 2d_q1_uniform_elas
771:     args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
772:   test:
773:     suffix: 2d_q2_uniform_elas
774:     requires: !single
775:     args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu
776:   test:
777:     suffix: 2d_q3_uniform_elas
778:     requires: !single
779:     args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu

781: TEST*/