Actual source code: ex77.c
petsc-3.7.5 2017-01-01
1: static char help[] = "Nonlinear elasticity problem in 3d with simplicial finite elements.\n\
2: We solve a nonlinear elasticity problem, modelled as an incompressible Neo-Hookean solid, \n\
3: with pressure loading in a rectangular domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";
5: /*
6: Nonlinear elasticity problem, which we discretize using the finite
7: element method on an unstructured mesh. This uses both Dirichlet boundary conditions (fixed faces)
8: and nonlinear Neumann boundary conditions (pressure loading).
9: The Lagrangian density (modulo boundary conditions) for this problem is given by
10: \begin{equation}
11: \frac{\mu}{2} (\mathrm{Tr}{C}-3) + J p + \frac{\kappa}{2} (J-1).
12: \end{equation}
14: Discretization:
16: We use PetscFE to generate a tabulation of the finite element basis functions
17: at quadrature points. We can currently generate an arbitrary order Lagrange
18: element.
20: Field Data:
22: DMPLEX data is organized by point, and the closure operation just stacks up the
23: data from each sieve point in the closure. Thus, for a P_2-P_1 Stokes element, we
24: have
26: cl{e} = {f e_0 e_1 e_2 v_0 v_1 v_2}
27: x = [u_{e_0} v_{e_0} u_{e_1} v_{e_1} u_{e_2} v_{e_2} u_{v_0} v_{v_0} p_{v_0} u_{v_1} v_{v_1} p_{v_1} u_{v_2} v_{v_2} p_{v_2}]
29: The problem here is that we would like to loop over each field separately for
30: integration. Therefore, the closure visitor in DMPlexVecGetClosure() reorders
31: the data so that each field is contiguous
33: x' = [u_{e_0} v_{e_0} u_{e_1} v_{e_1} u_{e_2} v_{e_2} u_{v_0} v_{v_0} u_{v_1} v_{v_1} u_{v_2} v_{v_2} p_{v_0} p_{v_1} p_{v_2}]
35: Likewise, DMPlexVecSetClosure() takes data partitioned by field, and correctly
36: puts it into the Sieve ordering.
38: */
40: #include <petscdmplex.h>
41: #include <petscsnes.h>
42: #include <petscds.h>
44: typedef enum {RUN_FULL, RUN_TEST} RunType;
46: typedef struct {
47: PetscInt debug; /* The debugging level */
48: RunType runType; /* Whether to run tests, or solve the full problem */
49: PetscLogEvent createMeshEvent;
50: PetscBool showInitial, showSolution;
51: /* Domain and mesh definition */
52: PetscInt dim; /* The topological mesh dimension */
53: PetscBool interpolate; /* Generate intermediate mesh elements */
54: PetscBool simplex; /* Use simplices or tensor product cells */
55: PetscReal refinementLimit; /* The largest allowable cell volume */
56: PetscBool testPartition; /* Use a fixed partitioning for testing */
57: PetscReal mu; /* The shear modulus */
58: } AppCtx;
60: /* Kronecker-delta */
61: static const PetscInt delta2D[2*2] = {1,0,0,1};
62: static const PetscInt delta3D[3*3] = {1,0,0,0,1,0,0,0,1};
64: /* Levi-Civita symbol */
65: static const PetscInt epsilon2D[2*2] = {0,1,-1,0};
66: static const PetscInt epsilon3D[3*3*3] = {0,0,0,0,0,1,0,-1,0,0,0,-1,0,0,0,1,0,0,0,1,0,-1,0,0,0,0,0};
68: PETSC_STATIC_INLINE void Det2D(PetscReal *detJ, const PetscReal J[])
69: {
70: *detJ = J[0]*J[3] - J[1]*J[2];
71: }
73: PETSC_STATIC_INLINE void Det3D(PetscReal *detJ, const PetscReal J[])
74: {
75: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
76: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
77: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
78: }
80: PETSC_STATIC_INLINE void Cof2D(PetscReal C[], const PetscReal A[])
81: {
82: C[0] = A[3];
83: C[1] = -A[2];
84: C[2] = -A[1];
85: C[3] = A[0];
86: }
88: PETSC_STATIC_INLINE void Cof3D(PetscReal C[], const PetscReal A[])
89: {
90: C[0*3+0] = A[1*3+1]*A[2*3+2] - A[1*3+2]*A[2*3+1];
91: C[0*3+1] = A[1*3+2]*A[2*3+0] - A[1*3+0]*A[2*3+2];
92: C[0*3+2] = A[1*3+0]*A[2*3+1] - A[1*3+1]*A[2*3+0];
93: C[1*3+0] = A[0*3+2]*A[2*3+1] - A[0*3+1]*A[2*3+2];
94: C[1*3+1] = A[0*3+0]*A[2*3+2] - A[0*3+2]*A[2*3+0];
95: C[1*3+2] = A[0*3+1]*A[2*3+0] - A[0*3+0]*A[2*3+1];
96: C[2*3+0] = A[0*3+1]*A[1*3+2] - A[0*3+2]*A[1*3+1];
97: C[2*3+1] = A[0*3+2]*A[1*3+0] - A[0*3+0]*A[1*3+2];
98: C[2*3+2] = A[0*3+0]*A[1*3+1] - A[0*3+1]*A[1*3+0];
99: }
101: PetscErrorCode zero_scalar(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
102: {
103: u[0] = 0.0;
104: return 0;
105: }
107: PetscErrorCode zero_vector(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
108: {
109: const PetscInt Ncomp = dim;
111: PetscInt comp;
112: for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0.0;
113: return 0;
114: }
116: PetscErrorCode coordinates(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
117: {
118: const PetscInt Ncomp = dim;
120: PetscInt comp;
121: for (comp = 0; comp < Ncomp; ++comp) u[comp] = x[comp];
122: return 0;
123: }
125: PetscErrorCode elasticityMaterial(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
126: {
127: AppCtx *user = (AppCtx *) ctx;
128: u[0] = user->mu;
129: return 0;
130: }
132: void f0_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[], PetscScalar f0[])
136: {
137: const PetscInt Ncomp = dim;
139: PetscInt comp;
140: for (comp = 0; comp < Ncomp; ++comp) f0[comp] = 0.0;
141: }
143: void f1_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
144: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
145: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
146: PetscReal t, const PetscReal x[], PetscScalar f1[])
147: {
148: const PetscInt Ncomp = dim;
149: const PetscReal mu = a[0], p = u[Ncomp], kappa = 3.0;
150: PetscReal cofu_x[Ncomp*dim], detu_x;
152: Cof3D(cofu_x, u_x);
153: Det3D(&detu_x, u_x);
155: /* f1 is the first Piola-Kirchhoff tensor */
156: PetscInt comp, d;
157: for (comp = 0; comp < Ncomp; ++comp) {
158: for (d = 0; d < dim; ++d) {
159: f1[comp*dim+d] = mu * u_x[comp*dim+d] + cofu_x[comp*dim+d] * (p + kappa * (detu_x - 1.0));
160: }
161: }
162: }
164: void g3_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
165: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
166: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
167: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
168: {
169: const PetscInt Ncomp = dim;
170: const PetscReal mu = a[0], p = u[Ncomp], kappa = 3.0;
171: PetscReal cofu_x[Ncomp*dim], detu_x;
173: Cof3D(cofu_x, u_x);
174: Det3D(&detu_x, u_x);
176: /* g3 is the first elasticity tensor, i.e. A_i^I_j^J = S^{IJ} g_{ij} + C^{KILJ} F^k_K F^l_L g_{kj} g_{li} */
177: PetscInt compI, compJ, compK, d1, d2, d3;
178: for (compI = 0; compI < Ncomp; ++compI) {
179: for (compJ = 0; compJ < Ncomp; ++compJ) {
180: for (d1 = 0; d1 < dim; ++d1) {
181: for (d2 = 0; d2 < dim; ++d2) {
182: for (compK = 0; compK < Ncomp; ++compK) {
183: for (d3 = 0; d3 < dim; ++d3) {
184: g3[((compI*Ncomp+compJ)*dim+d1)*dim+d2] += (p + kappa * (detu_x - 1.0)) * epsilon3D[(compI*Ncomp+compJ)*Ncomp+compK] * epsilon3D[(d1*dim+d2)*dim+d3] * u_x[compK*dim+d3];
185: }
186: }
187: g3[((compI*Ncomp+compJ)*dim+d1)*dim+d2] += kappa * cofu_x[compI*dim+d1] * cofu_x[compJ*dim+d2] + delta3D[compI*Ncomp+compJ] * delta3D[d1*dim+d2] * mu;
188: }
189: }
190: }
191: }
192: }
194: void f0_bd_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
195: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
196: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
197: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
198: {
199: const PetscInt Ncomp = dim;
200: const PetscReal p = 0.4; /* Should this be specified as an auxiliary field? */
201: PetscReal cofu_x[Ncomp*dim];
203: Cof3D(cofu_x, u_x);
205: PetscInt comp, d;
206: for (comp = 0; comp < Ncomp; ++comp) {
207: for (d = 0, f0[comp] = 0.0; d < dim; ++d) f0[comp] += cofu_x[comp*dim+d] * n[d];
208: f0[comp] *= p;
209: }
210: }
212: void f1_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
213: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
214: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
215: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[])
216: {
217: const PetscInt Ncomp = dim;
219: PetscInt comp, d;
220: for (comp = 0; comp < Ncomp; ++comp) {
221: for (d = 0; d < dim; ++d) {
222: f1[comp*dim+d] = 0.0;
223: }
224: }
225: }
227: void g1_bd_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
228: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
229: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
230: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscScalar g1[])
231: {
232: const PetscInt Ncomp = dim;
233: const PetscReal p = 0.4; /* Should this be specified as an auxiliary field? */
235: PetscInt compI, compJ, compK, d1, d2, d3;
236: for (compI = 0; compI < Ncomp; ++compI) {
237: for (compJ = 0; compJ < Ncomp; ++compJ) {
238: for (d2 = 0; d2 < dim; ++d2) {
239: for (compK = 0; compK < Ncomp; ++compK) {
240: for (d3 = 0; d3 < dim; ++d3) {
241: for (d1 = 0; d1 < dim; ++d1) {
242: g1[(compI*Ncomp+compJ)*dim+d2] += epsilon3D[(compI*Ncomp+compJ)*Ncomp+compK] * epsilon3D[(d2*dim+d3)*dim+d1] * u_x[compK*dim+d3] * n[d1];
243: }
244: }
245: }
246: g1[(compI*Ncomp+compJ)*dim+d2] *= p;
247: }
248: }
249: }
250: }
252: void f0_p_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
253: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
254: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
255: PetscReal t, const PetscReal x[], PetscScalar f0[])
256: {
257: PetscReal detu_x;
258: Det3D(&detu_x, u_x);
259: f0[0] = detu_x - 1.0;
260: }
262: void f1_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
263: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
264: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
265: PetscReal t, const PetscReal x[], PetscScalar f1[])
266: {
267: PetscInt d;
268: for (d = 0; d < dim; ++d) f1[d] = 0.0;
269: }
271: void f0_bd_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
272: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
273: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
274: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
275: {
276: f0[0] = 0.0;
277: }
278: void f1_bd_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
279: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
280: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
281: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[])
282: {
283: PetscInt d;
284: for (d = 0; d < dim; ++d) f1[d] = 0.0;
285: }
287: void g1_pu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
288: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
289: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
290: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g1[])
291: {
292: Cof3D(g1, u_x);
293: }
295: void g2_up_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
296: const PetscInt uOff[], const PetscInt uOff_x[], const PetscReal u[], const PetscReal u_t[], const PetscReal u_x[],
297: const PetscInt aOff[], const PetscInt aOff_x[], const PetscReal a[], const PetscReal a_t[], const PetscReal a_x[],
298: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscReal g2[])
299: {
300: Cof3D(g2, u_x);
301: }
305: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
306: {
307: const char *runTypes[2] = {"full", "test"};
308: PetscInt run;
312: options->debug = 0;
313: options->runType = RUN_FULL;
314: options->dim = 3;
315: options->interpolate = PETSC_FALSE;
316: options->simplex = PETSC_TRUE;
317: options->refinementLimit = 0.0;
318: options->mu = 1.0;
319: options->testPartition = PETSC_FALSE;
320: options->showInitial = PETSC_FALSE;
321: options->showSolution = PETSC_TRUE;
323: PetscOptionsBegin(comm, "", "Nonlinear elasticity problem options", "DMPLEX");
324: PetscOptionsInt("-debug", "The debugging level", "ex77.c", options->debug, &options->debug, NULL);
325: run = options->runType;
326: PetscOptionsEList("-run_type", "The run type", "ex77.c", runTypes, 2, runTypes[options->runType], &run, NULL);
328: options->runType = (RunType) run;
330: PetscOptionsInt("-dim", "The topological mesh dimension", "ex77.c", options->dim, &options->dim, NULL);
331: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex77.c", options->interpolate, &options->interpolate, NULL);
332: PetscOptionsBool("-simplex", "Use simplices or tensor product cells", "ex77.c", options->simplex, &options->simplex, NULL);
333: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex77.c", options->refinementLimit, &options->refinementLimit, NULL);
334: PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex77.c", options->testPartition, &options->testPartition, NULL);
335: PetscOptionsReal("-shear_modulus", "The shear modulus", "ex77.c", options->mu, &options->mu, NULL);
337: PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex77.c", options->showInitial, &options->showInitial, NULL);
338: PetscOptionsBool("-show_solution", "Output the solution for verification", "ex77.c", options->showSolution, &options->showSolution, NULL);
339: PetscOptionsEnd();
341: PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
342: return(0);
343: }
347: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
348: {
349: Vec lv;
350: PetscInt p;
351: PetscMPIInt rank, numProcs;
355: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
356: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
357: DMGetLocalVector(dm, &lv);
358: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
359: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
360: PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
361: for (p = 0; p < numProcs; ++p) {
362: if (p == rank) {VecView(lv, PETSC_VIEWER_STDOUT_SELF);}
363: PetscBarrier((PetscObject) dm);
364: }
365: DMRestoreLocalVector(dm, &lv);
366: return(0);
367: }
371: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
372: {
373: PetscInt dim = user->dim;
374: PetscBool interpolate = user->interpolate;
375: PetscReal refinementLimit = user->refinementLimit;
376: const PetscInt cells[3] = {3, 3, 3};
380: PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
381: if (user->simplex) {DMPlexCreateBoxMesh(comm, dim, interpolate, dm);}
382: else {DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);}
383: /* Label the faces (bit of a hack here, until it is properly implemented for simplices) */
384: {
385: DM cdm;
386: DMLabel label;
387: IS is;
388: PetscInt d, dim = user->dim, b, f, Nf;
389: const PetscInt *faces;
390: PetscInt csize;
391: PetscReal *coords = NULL;
392: PetscSection cs;
393: Vec coordinates ;
395: DMCreateLabel(*dm, "boundary");
396: DMGetLabel(*dm, "boundary", &label);
397: DMPlexMarkBoundaryFaces(*dm, label);
398: DMGetStratumIS(*dm, "boundary", 1, &is);
399: DMCreateLabel(*dm, "Faces");
400: if (is) {
401: ISGetLocalSize(is, &Nf);
402: ISGetIndices(is, &faces);
404: DMGetCoordinatesLocal(*dm, &coordinates);
405: DMGetCoordinateDM(*dm, &cdm);
406: DMGetDefaultSection(cdm, &cs);
408: /* Check for each boundary face if any component of its centroid is either 0.0 or 1.0 */
409: PetscReal faceCoord;
410: PetscInt v;
411: for (f = 0; f < Nf; ++f) {
412: DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
413: /* Calculate mean coordinate vector */
414: const PetscInt Nv = csize/dim;
415: for (d = 0; d < dim; ++d) {
416: faceCoord = 0.0;
417: for (v = 0; v < Nv; ++v) faceCoord += coords[v*dim+d];
418: faceCoord /= Nv;
419: for (b = 0; b < 2; ++b) {
420: if (PetscAbs(faceCoord - b*1.0) < PETSC_SMALL) {
421: DMSetLabelValue(*dm, "Faces", faces[f], d*2+b+1);
422: }
423: }
424: }
425: DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
426: }
427: ISRestoreIndices(is, &faces);
428: }
429: ISDestroy(&is);
430: DMGetLabel(*dm, "Faces", &label);
431: DMPlexLabelComplete(*dm, label);
432: }
433: {
434: DM refinedMesh = NULL;
435: DM distributedMesh = NULL;
437: /* Refine mesh using a volume constraint */
438: DMPlexSetRefinementLimit(*dm, refinementLimit);
439: if (user->simplex) {DMRefine(*dm, comm, &refinedMesh);}
440: if (refinedMesh) {
441: DMDestroy(dm);
442: *dm = refinedMesh;
443: }
444: /* Setup test partitioning */
445: if (user->testPartition) {
446: PetscInt triSizes_n2[2] = {4, 4};
447: PetscInt triPoints_n2[8] = {3, 5, 6, 7, 0, 1, 2, 4};
448: PetscInt triSizes_n3[3] = {2, 3, 3};
449: PetscInt triPoints_n3[8] = {3, 5, 1, 6, 7, 0, 2, 4};
450: PetscInt triSizes_n5[5] = {1, 2, 2, 1, 2};
451: PetscInt triPoints_n5[8] = {3, 5, 6, 4, 7, 0, 1, 2};
452: PetscInt triSizes_ref_n2[2] = {8, 8};
453: PetscInt triPoints_ref_n2[16] = {1, 5, 6, 7, 10, 11, 14, 15, 0, 2, 3, 4, 8, 9, 12, 13};
454: PetscInt triSizes_ref_n3[3] = {5, 6, 5};
455: PetscInt triPoints_ref_n3[16] = {1, 7, 10, 14, 15, 2, 6, 8, 11, 12, 13, 0, 3, 4, 5, 9};
456: PetscInt triSizes_ref_n5[5] = {3, 4, 3, 3, 3};
457: PetscInt triPoints_ref_n5[16] = {1, 7, 10, 2, 11, 13, 14, 5, 6, 15, 0, 8, 9, 3, 4, 12};
458: const PetscInt *sizes = NULL;
459: const PetscInt *points = NULL;
460: PetscPartitioner part;
461: PetscInt cEnd;
462: PetscMPIInt rank, numProcs;
464: MPI_Comm_rank(comm, &rank);
465: MPI_Comm_size(comm, &numProcs);
466: DMPlexGetHeightStratum(*dm, 0, NULL, &cEnd);
467: if (!rank) {
468: if (dim == 2 && user->simplex && numProcs == 2 && cEnd == 8) {
469: sizes = triSizes_n2; points = triPoints_n2;
470: } else if (dim == 2 && user->simplex && numProcs == 3 && cEnd == 8) {
471: sizes = triSizes_n3; points = triPoints_n3;
472: } else if (dim == 2 && user->simplex && numProcs == 5 && cEnd == 8) {
473: sizes = triSizes_n5; points = triPoints_n5;
474: } else if (dim == 2 && user->simplex && numProcs == 2 && cEnd == 16) {
475: sizes = triSizes_ref_n2; points = triPoints_ref_n2;
476: } else if (dim == 2 && user->simplex && numProcs == 3 && cEnd == 16) {
477: sizes = triSizes_ref_n3; points = triPoints_ref_n3;
478: } else if (dim == 2 && user->simplex && numProcs == 5 && cEnd == 16) {
479: sizes = triSizes_ref_n5; points = triPoints_ref_n5;
480: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "No stored partition matching run parameters");
481: }
482: DMPlexGetPartitioner(*dm, &part);
483: PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);
484: PetscPartitionerShellSetPartition(part, numProcs, sizes, points);
485: }
486: /* Distribute mesh over processes */
487: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
488: if (distributedMesh) {
489: DMDestroy(dm);
490: *dm = distributedMesh;
491: }
492: }
493: DMSetFromOptions(*dm);
494: DMViewFromOptions(*dm, NULL, "-dm_view");
495:
496: PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
497: return(0);
498: }
502: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
503: {
504: PetscDS prob;
508: DMGetDS(dm, &prob);
509: PetscDSSetResidual(prob, 0, f0_u, f1_u_3d);
510: PetscDSSetResidual(prob, 1, f0_p_3d, f1_p);
511: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu_3d);
512: PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_up_3d, NULL);
513: PetscDSSetJacobian(prob, 1, 0, NULL, g1_pu_3d, NULL, NULL);
515: PetscDSSetBdResidual(prob, 0, f0_bd_u_3d, f1_bd_u);
516: PetscDSSetBdResidual(prob, 1, f0_bd_p, f1_bd_p);
517: PetscDSSetBdJacobian(prob, 0, 0, NULL, g1_bd_uu_3d, NULL, NULL);
518: return(0);
519: }
523: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
524: {
525: PetscErrorCode (*matFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {elasticityMaterial};
526: Vec nu;
527: void *ctxs[] = {user};
531: DMCreateLocalVector(dmAux, &nu);
532: DMProjectFunctionLocal(dmAux, 0.0, matFuncs, ctxs, INSERT_ALL_VALUES, nu);
533: PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
534: VecDestroy(&nu);
535: return(0);
536: }
540: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
541: {
542: DM cdm = dm, coordDM, dmAux;
543: const PetscInt dim = user->dim;
544: const PetscBool simplex = user->simplex;
545: PetscFE fe[2], feBd[2], feAux;
546: PetscQuadrature q;
547: PetscDS prob, probAux;
548: PetscInt order;
549: PetscErrorCode ierr;
552: /* Create finite element */
553: PetscFECreateDefault(dm, dim, dim, simplex, "def_", -1, &fe[0]);
554: PetscObjectSetName((PetscObject) fe[0], "deformation");
555: PetscFEGetQuadrature(fe[0], &q);
556: PetscQuadratureGetOrder(q, &order);
557: PetscFECreateDefault(dm, dim, 1, simplex, "pres_", order, &fe[1]);
558: PetscObjectSetName((PetscObject) fe[1], "pressure");
559: PetscFECreateDefault(dm, dim-1, dim, simplex, "bd_def_", order, &feBd[0]);
560: PetscObjectSetName((PetscObject) feBd[0], "deformation");
561: PetscFECreateDefault(dm, dim-1, 1, simplex, "bd_pres_", order, &feBd[1]);
562: PetscObjectSetName((PetscObject) feBd[1], "pressure");
564: PetscFECreateDefault(dm, dim, 1, simplex, "elastMat_", order, &feAux);
565: PetscObjectSetName((PetscObject) feAux, "elasticityMaterial");
566:
568: /* Set discretization and boundary conditions for each mesh */
569: while (cdm) {
570: DMGetDS(cdm, &prob);
571: PetscDSSetDiscretization(prob, 0, (PetscObject) fe[0]);
572: PetscDSSetDiscretization(prob, 1, (PetscObject) fe[1]);
573: PetscDSSetBdDiscretization(prob, 0, (PetscObject) feBd[0]);
574: PetscDSSetBdDiscretization(prob, 1, (PetscObject) feBd[1]);
575: SetupProblem(cdm, user);
577: DMClone(cdm, &dmAux);
578: DMGetCoordinateDM(cdm, &coordDM);
579: DMSetCoordinateDM(dmAux, coordDM);
580: DMGetDS(dmAux, &probAux);
581: PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
582: PetscObjectCompose((PetscObject) cdm, "dmAux", (PetscObject) dmAux);
583: SetupMaterial(cdm, dmAux, user);
584: DMDestroy(&dmAux);
586: const PetscInt Ncomp = dim;
587: const PetscInt components[] = {0,1,2};
588: const PetscInt Nfid = 1;
589: const PetscInt fid[] = {1}; /* The fixed faces */
590: const PetscInt Npid = 1;
591: const PetscInt pid[] = {2}; /* The faces with pressure loading */
592: DMAddBoundary(cdm, PETSC_TRUE, "fixed", "Faces", 0, Ncomp, components, (void (*)()) coordinates, Nfid, fid, user);
593: DMAddBoundary(cdm, PETSC_FALSE, "pressure", "Faces", 0, Ncomp, components, NULL, Npid, pid, user);
594: DMGetCoarseDM(cdm, &cdm);
595: }
597: {
598: /* Set up the near null space (a.k.a. rigid body modes) that will be used by the multigrid preconditioner */
599: DM subdm;
600: MatNullSpace nearNullSpace;
601: PetscInt fields = 0;
602: PetscObject deformation;
603: DMCreateSubDM(dm, 1, &fields, NULL, &subdm);
604: DMPlexCreateRigidBody(subdm, &nearNullSpace);
605: DMGetField(dm, 0, &deformation);
606: PetscObjectCompose(deformation, "nearnullspace", (PetscObject) nearNullSpace);
607: DMDestroy(&subdm);
608: MatNullSpaceDestroy(&nearNullSpace);
609: }
611: PetscFEDestroy(&fe[0]);
612: PetscFEDestroy(&fe[1]);
613: PetscFEDestroy(&feBd[0]);
614: PetscFEDestroy(&feBd[1]);
615: PetscFEDestroy(&feAux);
616: return(0);
617: }
622: int main(int argc, char **argv)
623: {
624: SNES snes; /* nonlinear solver */
625: DM dm; /* problem definition */
626: Vec u,r; /* solution, residual vectors */
627: Mat A,J; /* Jacobian matrix */
628: MatNullSpace nullSpace; /* May be necessary for pressure */
629: AppCtx user; /* user-defined work context */
630: PetscInt its; /* iterations for convergence */
633: PetscInitialize(&argc, &argv, NULL, help);
634: ProcessOptions(PETSC_COMM_WORLD, &user);
635: SNESCreate(PETSC_COMM_WORLD, &snes);
636: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
637: SNESSetDM(snes, dm);
638: DMSetApplicationContext(dm, &user);
640: SetupDiscretization(dm, &user);
641: DMPlexCreateClosureIndex(dm, NULL);
643: DMCreateGlobalVector(dm, &u);
644: VecDuplicate(u, &r);
646: DMSetMatType(dm,MATAIJ);
647: DMCreateMatrix(dm, &J);
648: A = J;
650: DMPlexSetSNESLocalFEM(dm,&user,&user,&user);
651: SNESSetJacobian(snes, A, J, NULL, NULL);
653: SNESSetFromOptions(snes);
655: PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void* ctx) = {coordinates, zero_scalar};
656: DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
657: if (user.showInitial) {DMVecViewLocal(dm, u, PETSC_VIEWER_STDOUT_SELF);}
659: if (user.runType == RUN_FULL) {
660: if (user.debug) {
661: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
662: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
663: }
664: SNESSolve(snes, NULL, u);
665: SNESGetIterationNumber(snes, &its);
666: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
667: if (user.showSolution) {
668: PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
669: VecChop(u, 3.0e-9);
670: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
671: }
672: } else {
673: PetscReal res = 0.0;
675: /* Check initial guess */
676: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
677: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
678: /* Check residual */
679: SNESComputeFunction(snes, u, r);
680: PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
681: VecChop(r, 1.0e-10);
682: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
683: VecNorm(r, NORM_2, &res);
684: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
685: /* Check Jacobian */
686: {
687: Vec b;
688: PetscBool isNull;
690: SNESComputeJacobian(snes, u, A, A);
691: VecDuplicate(u, &b);
692: VecSet(r, 0.0);
693: SNESComputeFunction(snes, r, b);
694: MatMult(A, u, r);
695: VecAXPY(r, 1.0, b);
696: VecDestroy(&b);
697: PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
698: VecChop(r, 1.0e-10);
699: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
700: VecNorm(r, NORM_2, &res);
701: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
702: }
703: }
704: VecViewFromOptions(u, NULL, "-sol_vec_view");
706: if (A != J) {MatDestroy(&A);}
707: MatDestroy(&J);
708: VecDestroy(&u);
709: VecDestroy(&r);
710: SNESDestroy(&snes);
711: DMDestroy(&dm);
712: PetscFinalize();
713: return 0;
714: }