Actual source code: ex77.c
petsc-3.9.4 2018-09-11
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: PetscReal p_wall; /* The wall pressure */
59: } AppCtx;
61: #if 0
62: PETSC_STATIC_INLINE void Det2D(PetscReal *detJ, const PetscReal J[])
63: {
64: *detJ = J[0]*J[3] - J[1]*J[2];
65: }
66: #endif
68: PETSC_STATIC_INLINE void Det3D(PetscReal *detJ, const PetscScalar J[])
69: {
70: *detJ = PetscRealPart(J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
71: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
72: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
73: }
75: #if 0
76: PETSC_STATIC_INLINE void Cof2D(PetscReal C[], const PetscReal A[])
77: {
78: C[0] = A[3];
79: C[1] = -A[2];
80: C[2] = -A[1];
81: C[3] = A[0];
82: }
83: #endif
85: PETSC_STATIC_INLINE void Cof3D(PetscReal C[], const PetscScalar A[])
86: {
87: C[0*3+0] = PetscRealPart(A[1*3+1]*A[2*3+2] - A[1*3+2]*A[2*3+1]);
88: C[0*3+1] = PetscRealPart(A[1*3+2]*A[2*3+0] - A[1*3+0]*A[2*3+2]);
89: C[0*3+2] = PetscRealPart(A[1*3+0]*A[2*3+1] - A[1*3+1]*A[2*3+0]);
90: C[1*3+0] = PetscRealPart(A[0*3+2]*A[2*3+1] - A[0*3+1]*A[2*3+2]);
91: C[1*3+1] = PetscRealPart(A[0*3+0]*A[2*3+2] - A[0*3+2]*A[2*3+0]);
92: C[1*3+2] = PetscRealPart(A[0*3+1]*A[2*3+0] - A[0*3+0]*A[2*3+1]);
93: C[2*3+0] = PetscRealPart(A[0*3+1]*A[1*3+2] - A[0*3+2]*A[1*3+1]);
94: C[2*3+1] = PetscRealPart(A[0*3+2]*A[1*3+0] - A[0*3+0]*A[1*3+2]);
95: C[2*3+2] = PetscRealPart(A[0*3+0]*A[1*3+1] - A[0*3+1]*A[1*3+0]);
96: }
98: PetscErrorCode zero_scalar(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
99: {
100: u[0] = 0.0;
101: return 0;
102: }
104: PetscErrorCode zero_vector(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
105: {
106: const PetscInt Ncomp = dim;
108: PetscInt comp;
109: for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0.0;
110: return 0;
111: }
113: PetscErrorCode coordinates(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
114: {
115: const PetscInt Ncomp = dim;
117: PetscInt comp;
118: for (comp = 0; comp < Ncomp; ++comp) u[comp] = x[comp];
119: return 0;
120: }
122: PetscErrorCode elasticityMaterial(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
123: {
124: AppCtx *user = (AppCtx *) ctx;
125: u[0] = user->mu;
126: return 0;
127: }
129: PetscErrorCode wallPressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
130: {
131: AppCtx *user = (AppCtx *) ctx;
132: u[0] = user->p_wall;
133: return 0;
134: }
136: void f1_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
137: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
138: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
139: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
140: {
141: const PetscInt Ncomp = dim;
142: const PetscReal mu = PetscRealPart(a[0]), kappa = 3.0;
143: PetscReal cofu_x[9/*Ncomp*dim*/], detu_x, p = PetscRealPart(u[Ncomp]);
144: PetscInt comp, d;
146: Cof3D(cofu_x, u_x);
147: Det3D(&detu_x, u_x);
148: p += kappa * (detu_x - 1.0);
150: /* f1 is the first Piola-Kirchhoff tensor */
151: for (comp = 0; comp < Ncomp; ++comp) {
152: for (d = 0; d < dim; ++d) {
153: f1[comp*dim+d] = mu * u_x[comp*dim+d] + p * cofu_x[comp*dim+d];
154: }
155: }
156: }
158: void g3_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
159: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
160: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
161: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
162: {
163: const PetscInt Ncomp = dim;
164: const PetscReal mu = PetscRealPart(a[0]), kappa = 3.0;
165: PetscReal cofu_x[9/*Ncomp*dim*/], detu_x, pp, pm, p = PetscRealPart(u[Ncomp]);
166: PetscInt compI, compJ, d1, d2;
168: Cof3D(cofu_x, u_x);
169: Det3D(&detu_x, u_x);
170: p += kappa * (detu_x - 1.0);
171: pp = p/detu_x + kappa;
172: pm = p/detu_x;
174: /* 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} */
175: for (compI = 0; compI < Ncomp; ++compI) {
176: for (compJ = 0; compJ < Ncomp; ++compJ) {
177: const PetscReal G = (compI == compJ) ? 1.0 : 0.0;
179: for (d1 = 0; d1 < dim; ++d1) {
180: for (d2 = 0; d2 < dim; ++d2) {
181: const PetscReal g = (d1 == d2) ? 1.0 : 0.0;
183: g3[((compI*Ncomp+compJ)*dim+d1)*dim+d2] = g * G * mu + pp * cofu_x[compI*dim+d1] * cofu_x[compJ*dim+d2] - pm * cofu_x[compI*dim+d2] * cofu_x[compJ*dim+d1];
184: }
185: }
186: }
187: }
188: }
190: void f0_bd_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
191: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
192: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
193: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
194: {
195: const PetscInt Ncomp = dim;
196: const PetscScalar p = a[aOff[1]];
197: PetscReal cofu_x[9/*Ncomp*dim*/];
198: PetscInt comp, d;
200: Cof3D(cofu_x, u_x);
201: for (comp = 0; comp < Ncomp; ++comp) {
202: for (d = 0, f0[comp] = 0.0; d < dim; ++d) f0[comp] += cofu_x[comp*dim+d] * n[d];
203: f0[comp] *= p;
204: }
205: }
207: void g1_bd_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
208: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
209: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
210: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
211: {
212: const PetscInt Ncomp = dim;
213: PetscScalar p = a[aOff[1]];
214: PetscReal cofu_x[9/*Ncomp*dim*/], m[3/*Ncomp*/], detu_x;
215: PetscInt comp, compI, compJ, d;
217: Cof3D(cofu_x, u_x);
218: Det3D(&detu_x, u_x);
219: p /= detu_x;
221: for (comp = 0; comp < Ncomp; ++comp) for (d = 0, m[comp] = 0.0; d < dim; ++d) m[comp] += cofu_x[comp*dim+d] * n[d];
222: for (compI = 0; compI < Ncomp; ++compI) {
223: for (compJ = 0; compJ < Ncomp; ++compJ) {
224: for (d = 0; d < dim; ++d) {
225: g1[(compI*Ncomp+compJ)*dim+d] = p * (m[compI] * cofu_x[compJ*dim+d] - cofu_x[compI*dim+d] * m[compJ]);
226: }
227: }
228: }
229: }
231: void f0_p_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
232: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
233: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
234: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
235: {
236: PetscReal detu_x;
237: Det3D(&detu_x, u_x);
238: f0[0] = detu_x - 1.0;
239: }
241: void g1_pu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
242: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
243: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
244: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
245: {
246: PetscReal cofu_x[9/*Ncomp*dim*/];
247: PetscInt compI, d;
249: Cof3D(cofu_x, u_x);
250: for (compI = 0; compI < dim; ++compI)
251: for (d = 0; d < dim; ++d) g1[compI*dim+d] = cofu_x[compI*dim+d];
252: }
254: void g2_up_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
255: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
256: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
257: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
258: {
259: PetscReal cofu_x[9/*Ncomp*dim*/];
260: PetscInt compI, d;
262: Cof3D(cofu_x, u_x);
263: for (compI = 0; compI < dim; ++compI)
264: for (d = 0; d < dim; ++d) g2[compI*dim+d] = cofu_x[compI*dim+d];
265: }
267: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
268: {
269: const char *runTypes[2] = {"full", "test"};
270: PetscInt run;
274: options->debug = 0;
275: options->runType = RUN_FULL;
276: options->dim = 3;
277: options->interpolate = PETSC_FALSE;
278: options->simplex = PETSC_TRUE;
279: options->refinementLimit = 0.0;
280: options->mu = 1.0;
281: options->p_wall = 0.4;
282: options->testPartition = PETSC_FALSE;
283: options->showInitial = PETSC_FALSE;
284: options->showSolution = PETSC_TRUE;
286: PetscOptionsBegin(comm, "", "Nonlinear elasticity problem options", "DMPLEX");
287: PetscOptionsInt("-debug", "The debugging level", "ex77.c", options->debug, &options->debug, NULL);
288: run = options->runType;
289: PetscOptionsEList("-run_type", "The run type", "ex77.c", runTypes, 2, runTypes[options->runType], &run, NULL);
291: options->runType = (RunType) run;
293: PetscOptionsInt("-dim", "The topological mesh dimension", "ex77.c", options->dim, &options->dim, NULL);
294: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex77.c", options->interpolate, &options->interpolate, NULL);
295: PetscOptionsBool("-simplex", "Use simplices or tensor product cells", "ex77.c", options->simplex, &options->simplex, NULL);
296: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex77.c", options->refinementLimit, &options->refinementLimit, NULL);
297: PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex77.c", options->testPartition, &options->testPartition, NULL);
298: PetscOptionsReal("-shear_modulus", "The shear modulus", "ex77.c", options->mu, &options->mu, NULL);
299: PetscOptionsReal("-wall_pressure", "The wall pressure", "ex77.c", options->p_wall, &options->p_wall, NULL);
301: PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex77.c", options->showInitial, &options->showInitial, NULL);
302: PetscOptionsBool("-show_solution", "Output the solution for verification", "ex77.c", options->showSolution, &options->showSolution, NULL);
303: PetscOptionsEnd();
305: PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
306: return(0);
307: }
309: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
310: {
311: Vec lv;
312: PetscInt p;
313: PetscMPIInt rank, size;
317: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
318: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
319: DMGetLocalVector(dm, &lv);
320: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
321: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
322: PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
323: for (p = 0; p < size; ++p) {
324: if (p == rank) {VecView(lv, PETSC_VIEWER_STDOUT_SELF);}
325: PetscBarrier((PetscObject) dm);
326: }
327: DMRestoreLocalVector(dm, &lv);
328: return(0);
329: }
331: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
332: {
333: PetscInt dim = user->dim;
334: PetscBool interpolate = user->interpolate;
335: PetscReal refinementLimit = user->refinementLimit;
336: const PetscInt cells[3] = {3, 3, 3};
340: PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
341: DMPlexCreateBoxMesh(comm, dim, user->simplex, user->simplex ? NULL : cells, NULL, NULL, NULL, interpolate, dm);
342: /* Label the faces (bit of a hack here, until it is properly implemented for simplices) */
343: {
344: DM cdm;
345: DMLabel label;
346: IS is;
347: PetscInt d, dim = user->dim, b, f, Nf;
348: const PetscInt *faces;
349: PetscInt csize;
350: PetscScalar *coords = NULL;
351: PetscSection cs;
352: Vec coordinates ;
354: DMCreateLabel(*dm, "boundary");
355: DMGetLabel(*dm, "boundary", &label);
356: DMPlexMarkBoundaryFaces(*dm, 1, label);
357: DMGetStratumIS(*dm, "boundary", 1, &is);
358: if (is) {
359: PetscReal faceCoord;
360: PetscInt v;
362: ISGetLocalSize(is, &Nf);
363: ISGetIndices(is, &faces);
365: DMGetCoordinatesLocal(*dm, &coordinates);
366: DMGetCoordinateDM(*dm, &cdm);
367: DMGetDefaultSection(cdm, &cs);
369: /* Check for each boundary face if any component of its centroid is either 0.0 or 1.0 */
370: for (f = 0; f < Nf; ++f) {
371: DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
372: /* Calculate mean coordinate vector */
373: for (d = 0; d < dim; ++d) {
374: const PetscInt Nv = csize/dim;
375: faceCoord = 0.0;
376: for (v = 0; v < Nv; ++v) faceCoord += PetscRealPart(coords[v*dim+d]);
377: faceCoord /= Nv;
378: for (b = 0; b < 2; ++b) {
379: if (PetscAbs(faceCoord - b*1.0) < PETSC_SMALL) {
380: DMSetLabelValue(*dm, "Faces", faces[f], d*2+b+1);
381: }
382: }
383: }
384: DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
385: }
386: ISRestoreIndices(is, &faces);
387: }
388: ISDestroy(&is);
389: }
390: {
391: DM refinedMesh = NULL;
392: DM distributedMesh = NULL;
393: PetscPartitioner part;
395: DMPlexGetPartitioner(*dm, &part);
397: /* Refine mesh using a volume constraint */
398: DMPlexSetRefinementLimit(*dm, refinementLimit);
399: if (user->simplex) {DMRefine(*dm, comm, &refinedMesh);}
400: if (refinedMesh) {
401: DMDestroy(dm);
402: *dm = refinedMesh;
403: }
404: /* Setup test partitioning */
405: if (user->testPartition) {
406: PetscInt triSizes_n2[2] = {4, 4};
407: PetscInt triPoints_n2[8] = {3, 5, 6, 7, 0, 1, 2, 4};
408: PetscInt triSizes_n3[3] = {2, 3, 3};
409: PetscInt triPoints_n3[8] = {3, 5, 1, 6, 7, 0, 2, 4};
410: PetscInt triSizes_n5[5] = {1, 2, 2, 1, 2};
411: PetscInt triPoints_n5[8] = {3, 5, 6, 4, 7, 0, 1, 2};
412: PetscInt triSizes_ref_n2[2] = {8, 8};
413: PetscInt triPoints_ref_n2[16] = {1, 5, 6, 7, 10, 11, 14, 15, 0, 2, 3, 4, 8, 9, 12, 13};
414: PetscInt triSizes_ref_n3[3] = {5, 6, 5};
415: PetscInt triPoints_ref_n3[16] = {1, 7, 10, 14, 15, 2, 6, 8, 11, 12, 13, 0, 3, 4, 5, 9};
416: PetscInt triSizes_ref_n5[5] = {3, 4, 3, 3, 3};
417: PetscInt triPoints_ref_n5[16] = {1, 7, 10, 2, 11, 13, 14, 5, 6, 15, 0, 8, 9, 3, 4, 12};
418: PetscInt tetSizes_n2[2] = {3, 3};
419: PetscInt tetPoints_n2[6] = {1, 2, 3, 0, 4, 5};
420: const PetscInt *sizes = NULL;
421: const PetscInt *points = NULL;
422: PetscInt cEnd;
423: PetscMPIInt rank, size;
425: MPI_Comm_rank(comm, &rank);
426: MPI_Comm_size(comm, &size);
427: DMPlexGetHeightStratum(*dm, 0, NULL, &cEnd);
428: if (!rank) {
429: if (dim == 2 && user->simplex && size == 2 && cEnd == 8) {
430: sizes = triSizes_n2; points = triPoints_n2;
431: } else if (dim == 2 && user->simplex && size == 3 && cEnd == 8) {
432: sizes = triSizes_n3; points = triPoints_n3;
433: } else if (dim == 2 && user->simplex && size == 5 && cEnd == 8) {
434: sizes = triSizes_n5; points = triPoints_n5;
435: } else if (dim == 2 && user->simplex && size == 2 && cEnd == 16) {
436: sizes = triSizes_ref_n2; points = triPoints_ref_n2;
437: } else if (dim == 2 && user->simplex && size == 3 && cEnd == 16) {
438: sizes = triSizes_ref_n3; points = triPoints_ref_n3;
439: } else if (dim == 2 && user->simplex && size == 5 && cEnd == 16) {
440: sizes = triSizes_ref_n5; points = triPoints_ref_n5;
441: } else if (dim == 3 && user->simplex && size == 2 && cEnd == 6) {
442: sizes = tetSizes_n2; points = tetPoints_n2;
443: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "No stored partition matching run parameters");
444: }
445: PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);
446: PetscPartitionerShellSetPartition(part, size, sizes, points);
447: } else {
448: PetscPartitionerSetFromOptions(part);
449: }
450: /* Distribute mesh over processes */
451: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
452: if (distributedMesh) {
453: DMDestroy(dm);
454: *dm = distributedMesh;
455: }
456: }
457: DMSetFromOptions(*dm);
458: DMViewFromOptions(*dm, NULL, "-dm_view");
459:
460: PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
461: return(0);
462: }
464: PetscErrorCode SetupProblem(PetscDS prob, PetscInt dim, AppCtx *user)
465: {
469: PetscDSSetResidual(prob, 0, NULL, f1_u_3d);
470: PetscDSSetResidual(prob, 1, f0_p_3d, NULL);
471: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu_3d);
472: PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_up_3d, NULL);
473: PetscDSSetJacobian(prob, 1, 0, NULL, g1_pu_3d, NULL, NULL);
475: PetscDSSetBdResidual(prob, 0, f0_bd_u_3d, NULL);
476: PetscDSSetBdJacobian(prob, 0, 0, NULL, g1_bd_uu_3d, NULL, NULL);
478: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "fixed", "Faces", 0, 0, NULL, (void (*)(void)) coordinates, 0, NULL, user);
479: PetscDSAddBoundary(prob, DM_BC_NATURAL, "pressure", "Faces", 0, 0, NULL, NULL, 0, NULL, user);
480: return(0);
481: }
483: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
484: {
485: PetscErrorCode (*matFuncs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {elasticityMaterial, wallPressure};
486: Vec A;
487: void *ctxs[2];
491: ctxs[0] = user; ctxs[1] = user;
492: DMCreateLocalVector(dmAux, &A);
493: DMProjectFunctionLocal(dmAux, 0.0, matFuncs, ctxs, INSERT_ALL_VALUES, A);
494: PetscObjectCompose((PetscObject) dm, "A", (PetscObject) A);
495: VecDestroy(&A);
496: return(0);
497: }
499: PetscErrorCode SetupNearNullSpace(DM dm, AppCtx *user)
500: {
501: /* Set up the near null space (a.k.a. rigid body modes) that will be used by the multigrid preconditioner */
502: DM subdm;
503: MatNullSpace nearNullSpace;
504: PetscInt fields = 0;
505: PetscObject deformation;
509: DMCreateSubDM(dm, 1, &fields, NULL, &subdm);
510: DMPlexCreateRigidBody(subdm, &nearNullSpace);
511: DMGetField(dm, 0, &deformation);
512: PetscObjectCompose(deformation, "nearnullspace", (PetscObject) nearNullSpace);
513: DMDestroy(&subdm);
514: MatNullSpaceDestroy(&nearNullSpace);
515: return(0);
516: }
518: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
519: {
520: DM cdm = dm, dmAux;
521: const PetscInt dim = user->dim;
522: const PetscBool simplex = user->simplex;
523: PetscFE fe[2], feAux[2];
524: PetscQuadrature q, fq;
525: PetscDS prob, probAux;
526: PetscErrorCode ierr;
529: /* Create finite element */
530: PetscFECreateDefault(dm, dim, dim, simplex, "def_", PETSC_DEFAULT, &fe[0]);
531: PetscObjectSetName((PetscObject) fe[0], "deformation");
532: PetscFEGetQuadrature(fe[0], &q);
533: PetscFEGetFaceQuadrature(fe[0], &fq);
534: PetscFECreateDefault(dm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);
535: PetscFESetQuadrature(fe[1], q);
536: PetscFESetFaceQuadrature(fe[1], fq);
538: PetscObjectSetName((PetscObject) fe[1], "pressure");
540: PetscFECreateDefault(dm, dim, 1, simplex, "elastMat_", PETSC_DEFAULT, &feAux[0]);
541: PetscObjectSetName((PetscObject) feAux[0], "elasticityMaterial");
542: PetscFESetQuadrature(feAux[0], q);
543: PetscFESetFaceQuadrature(feAux[0], fq);
544: /* It is not yet possible to define a field on a submesh (e.g. a boundary), so we will use a normal finite element */
545: PetscFECreateDefault(dm, dim, 1, simplex, "wall_pres_", PETSC_DEFAULT, &feAux[1]);
546: PetscObjectSetName((PetscObject) feAux[1], "wall_pressure");
547: PetscFESetQuadrature(feAux[1], q);
548: PetscFESetFaceQuadrature(feAux[1], fq);
550: /* Set discretization and boundary conditions for each mesh */
551: DMGetDS(dm, &prob);
552: PetscDSSetDiscretization(prob, 0, (PetscObject) fe[0]);
553: PetscDSSetDiscretization(prob, 1, (PetscObject) fe[1]);
554: SetupProblem(prob, dim, user);
555: PetscDSSetFromOptions(prob);
556: PetscDSCreate(PetscObjectComm((PetscObject)dm),&probAux);
557: PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux[0]);
558: PetscDSSetDiscretization(probAux, 1, (PetscObject) feAux[1]);
559: PetscDSSetFromOptions(probAux);
560: while (cdm) {
561: DMSetDS(cdm, prob);
562: DMClone(cdm, &dmAux);
563: DMSetDS(dmAux, probAux);
564: PetscObjectCompose((PetscObject) cdm, "dmAux", (PetscObject) dmAux);
565: SetupMaterial(cdm, dmAux, user);
566: DMDestroy(&dmAux);
568: DMGetCoarseDM(cdm, &cdm);
569: }
570: PetscDSDestroy(&probAux);
571: PetscFEDestroy(&fe[0]);
572: PetscFEDestroy(&fe[1]);
573: PetscFEDestroy(&feAux[0]);
574: PetscFEDestroy(&feAux[1]);
575: return(0);
576: }
579: int main(int argc, char **argv)
580: {
581: SNES snes; /* nonlinear solver */
582: DM dm; /* problem definition */
583: Vec u,r; /* solution, residual vectors */
584: Mat A,J; /* Jacobian matrix */
585: AppCtx user; /* user-defined work context */
586: PetscInt its; /* iterations for convergence */
589: PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
590: ProcessOptions(PETSC_COMM_WORLD, &user);
591: SNESCreate(PETSC_COMM_WORLD, &snes);
592: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
593: SNESSetDM(snes, dm);
594: DMSetApplicationContext(dm, &user);
596: SetupDiscretization(dm, &user);
597: DMPlexCreateClosureIndex(dm, NULL);
598: SetupNearNullSpace(dm, &user);
600: DMCreateGlobalVector(dm, &u);
601: VecDuplicate(u, &r);
603: DMSetMatType(dm,MATAIJ);
604: DMCreateMatrix(dm, &J);
605: A = J;
607: DMPlexSetSNESLocalFEM(dm,&user,&user,&user);
608: SNESSetJacobian(snes, A, J, NULL, NULL);
610: SNESSetFromOptions(snes);
612: {
613: PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void* ctx);
614: initialGuess[0] = coordinates; initialGuess[1] = zero_scalar;
615: DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
616: }
617: if (user.showInitial) {DMVecViewLocal(dm, u, PETSC_VIEWER_STDOUT_SELF);}
619: if (user.runType == RUN_FULL) {
620: if (user.debug) {
621: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
622: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
623: }
624: SNESSolve(snes, NULL, u);
625: SNESGetIterationNumber(snes, &its);
626: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
627: if (user.showSolution) {
628: PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
629: VecChop(u, 3.0e-9);
630: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
631: }
632: } else {
633: PetscReal res = 0.0;
635: /* Check initial guess */
636: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
637: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
638: /* Check residual */
639: SNESComputeFunction(snes, u, r);
640: PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
641: VecChop(r, 1.0e-10);
642: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
643: VecNorm(r, NORM_2, &res);
644: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
645: /* Check Jacobian */
646: {
647: Vec b;
649: SNESComputeJacobian(snes, u, A, A);
650: VecDuplicate(u, &b);
651: VecSet(r, 0.0);
652: SNESComputeFunction(snes, r, b);
653: MatMult(A, u, r);
654: VecAXPY(r, 1.0, b);
655: VecDestroy(&b);
656: PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
657: VecChop(r, 1.0e-10);
658: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
659: VecNorm(r, NORM_2, &res);
660: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
661: }
662: }
663: VecViewFromOptions(u, NULL, "-sol_vec_view");
665: if (A != J) {MatDestroy(&A);}
666: MatDestroy(&J);
667: VecDestroy(&u);
668: VecDestroy(&r);
669: SNESDestroy(&snes);
670: DMDestroy(&dm);
671: PetscFinalize();
672: return ierr;
673: }
675: /*TEST
677: build:
678: requires: !complex
680: test:
681: suffix: 0
682: requires: ctetgen !single
683: args: -run_type full -dim 3 -dm_refine 2 -interpolate 1 -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4 -def_petscspace_order 2 -pres_petscspace_order 1 -elastMat_petscspace_order 0 -wall_pres_petscspace_order 0 -snes_rtol 1e-05 -ksp_type fgmres -ksp_rtol 1e-10 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -ksp_converged_reason -show_solution 0
685: test:
686: suffix: 1
687: requires: ctetgen superlu_dist
688: nsize: 2
689: args: -run_type full -dim 3 -dm_refine 0 -interpolate 1 -test_partition -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4 -def_petscspace_order 2 -pres_petscspace_order 1 -elastMat_petscspace_order 0 -wall_pres_petscspace_order 0 -snes_rtol 1e-05 -ksp_type fgmres -ksp_rtol 1e-10 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -ksp_converged_reason -show_solution 0
690: timeoutfactor: 2
692: test:
693: suffix: 2
694: requires: ctetgen !single
695: args: -run_type full -dim 3 -dm_refine 2 -interpolate 1 -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0 -def_petscspace_order 2 -pres_petscspace_order 1 -elastMat_petscspace_order 0 -wall_pres_petscspace_order 0 -snes_rtol 1e-05 -ksp_type fgmres -ksp_rtol 1e-10 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -ksp_converged_reason -show_solution 0
697: test:
698: requires: ctetgen !single
699: args: -run_type full -dim 3 -dm_refine 2 -interpolate 1 -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4 -def_petscspace_order 2 -pres_petscspace_order 1 -elastMat_petscspace_order 0 -wall_pres_petscspace_order 0 -snes_rtol 1e-05 -ksp_type fgmres -ksp_rtol 1e-10 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -ksp_converged_reason -show_solution 0
700: output_file: output/ex77_0.out
702: test:
703: requires: ctetgen superlu_dist
704: suffix: 4
705: nsize: 2
706: args: -run_type full -dim 3 -dm_refine 0 -interpolate 1 -test_partition -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4 -def_petscspace_order 2 -pres_petscspace_order 1 -elastMat_petscspace_order 0 -wall_pres_petscspace_order 0 -snes_rtol 1e-05 -ksp_type fgmres -ksp_rtol 1e-10 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -ksp_converged_reason -show_solution 0
707: output_file: output/ex77_1.out
709: test:
710: requires: ctetgen !single
711: suffix: 3
712: args: -run_type full -dim 3 -dm_refine 2 -interpolate 1 -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0 -def_petscspace_order 2 -pres_petscspace_order 1 -elastMat_petscspace_order 0 -wall_pres_petscspace_order 0 -snes_rtol 1e-05 -ksp_type fgmres -ksp_rtol 1e-10 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -ksp_converged_reason -show_solution 0
713: output_file: output/ex77_2.out
715: test:
716: requires: ctetgen !single
717: suffix: 2_par
718: nsize: 4
719: args: -run_type full -dim 3 -dm_refine 2 -interpolate 1 -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0 -def_petscspace_order 2 -pres_petscspace_order 1 -elastMat_petscspace_order 0 -wall_pres_petscspace_order 0 -snes_rtol 1e-05 -ksp_type fgmres -ksp_rtol 1e-10 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -ksp_converged_reason -show_solution 0
720: output_file: output/ex77_2.out
722: TEST*/