Actual source code: ex77.c
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: RunType runType; /* Whether to run tests, or solve the full problem */
48: PetscReal mu; /* The shear modulus */
49: PetscReal p_wall; /* The wall pressure */
50: } AppCtx;
52: #if 0
53: PETSC_STATIC_INLINE void Det2D(PetscReal *detJ, const PetscReal J[])
54: {
55: *detJ = J[0]*J[3] - J[1]*J[2];
56: }
57: #endif
59: PETSC_STATIC_INLINE void Det3D(PetscReal *detJ, const PetscScalar J[])
60: {
61: *detJ = PetscRealPart(J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
62: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
63: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
64: }
66: #if 0
67: PETSC_STATIC_INLINE void Cof2D(PetscReal C[], const PetscReal A[])
68: {
69: C[0] = A[3];
70: C[1] = -A[2];
71: C[2] = -A[1];
72: C[3] = A[0];
73: }
74: #endif
76: PETSC_STATIC_INLINE void Cof3D(PetscReal C[], const PetscScalar A[])
77: {
78: C[0*3+0] = PetscRealPart(A[1*3+1]*A[2*3+2] - A[1*3+2]*A[2*3+1]);
79: C[0*3+1] = PetscRealPart(A[1*3+2]*A[2*3+0] - A[1*3+0]*A[2*3+2]);
80: C[0*3+2] = PetscRealPart(A[1*3+0]*A[2*3+1] - A[1*3+1]*A[2*3+0]);
81: C[1*3+0] = PetscRealPart(A[0*3+2]*A[2*3+1] - A[0*3+1]*A[2*3+2]);
82: C[1*3+1] = PetscRealPart(A[0*3+0]*A[2*3+2] - A[0*3+2]*A[2*3+0]);
83: C[1*3+2] = PetscRealPart(A[0*3+1]*A[2*3+0] - A[0*3+0]*A[2*3+1]);
84: C[2*3+0] = PetscRealPart(A[0*3+1]*A[1*3+2] - A[0*3+2]*A[1*3+1]);
85: C[2*3+1] = PetscRealPart(A[0*3+2]*A[1*3+0] - A[0*3+0]*A[1*3+2]);
86: C[2*3+2] = PetscRealPart(A[0*3+0]*A[1*3+1] - A[0*3+1]*A[1*3+0]);
87: }
89: PetscErrorCode zero_scalar(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
90: {
91: u[0] = 0.0;
92: return 0;
93: }
95: PetscErrorCode zero_vector(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
96: {
97: const PetscInt Ncomp = dim;
99: PetscInt comp;
100: for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0.0;
101: return 0;
102: }
104: PetscErrorCode coordinates(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] = x[comp];
110: return 0;
111: }
113: PetscErrorCode elasticityMaterial(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
114: {
115: AppCtx *user = (AppCtx *) ctx;
116: u[0] = user->mu;
117: return 0;
118: }
120: PetscErrorCode wallPressure(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
121: {
122: AppCtx *user = (AppCtx *) ctx;
123: u[0] = user->p_wall;
124: return 0;
125: }
127: void f1_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
128: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
129: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
130: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
131: {
132: const PetscInt Ncomp = dim;
133: const PetscReal mu = PetscRealPart(a[0]), kappa = 3.0;
134: PetscReal cofu_x[9/*Ncomp*dim*/], detu_x, p = PetscRealPart(u[Ncomp]);
135: PetscInt comp, d;
137: Cof3D(cofu_x, u_x);
138: Det3D(&detu_x, u_x);
139: p += kappa * (detu_x - 1.0);
141: /* f1 is the first Piola-Kirchhoff tensor */
142: for (comp = 0; comp < Ncomp; ++comp) {
143: for (d = 0; d < dim; ++d) {
144: f1[comp*dim+d] = mu * u_x[comp*dim+d] + p * cofu_x[comp*dim+d];
145: }
146: }
147: }
149: void g3_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
150: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
151: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
152: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
153: {
154: const PetscInt Ncomp = dim;
155: const PetscReal mu = PetscRealPart(a[0]), kappa = 3.0;
156: PetscReal cofu_x[9/*Ncomp*dim*/], detu_x, pp, pm, p = PetscRealPart(u[Ncomp]);
157: PetscInt compI, compJ, d1, d2;
159: Cof3D(cofu_x, u_x);
160: Det3D(&detu_x, u_x);
161: p += kappa * (detu_x - 1.0);
162: pp = p/detu_x + kappa;
163: pm = p/detu_x;
165: /* 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} */
166: for (compI = 0; compI < Ncomp; ++compI) {
167: for (compJ = 0; compJ < Ncomp; ++compJ) {
168: const PetscReal G = (compI == compJ) ? 1.0 : 0.0;
170: for (d1 = 0; d1 < dim; ++d1) {
171: for (d2 = 0; d2 < dim; ++d2) {
172: const PetscReal g = (d1 == d2) ? 1.0 : 0.0;
174: 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];
175: }
176: }
177: }
178: }
179: }
181: void f0_bd_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
182: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
183: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
184: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
185: {
186: const PetscInt Ncomp = dim;
187: const PetscScalar p = a[aOff[1]];
188: PetscReal cofu_x[9/*Ncomp*dim*/];
189: PetscInt comp, d;
191: Cof3D(cofu_x, u_x);
192: for (comp = 0; comp < Ncomp; ++comp) {
193: for (d = 0, f0[comp] = 0.0; d < dim; ++d) f0[comp] += cofu_x[comp*dim+d] * n[d];
194: f0[comp] *= p;
195: }
196: }
198: void g1_bd_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
199: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
200: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
201: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
202: {
203: const PetscInt Ncomp = dim;
204: PetscScalar p = a[aOff[1]];
205: PetscReal cofu_x[9/*Ncomp*dim*/], m[3/*Ncomp*/], detu_x;
206: PetscInt comp, compI, compJ, d;
208: Cof3D(cofu_x, u_x);
209: Det3D(&detu_x, u_x);
210: p /= detu_x;
212: 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];
213: for (compI = 0; compI < Ncomp; ++compI) {
214: for (compJ = 0; compJ < Ncomp; ++compJ) {
215: for (d = 0; d < dim; ++d) {
216: g1[(compI*Ncomp+compJ)*dim+d] = p * (m[compI] * cofu_x[compJ*dim+d] - cofu_x[compI*dim+d] * m[compJ]);
217: }
218: }
219: }
220: }
222: void f0_p_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
223: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
224: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
225: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
226: {
227: PetscReal detu_x;
228: Det3D(&detu_x, u_x);
229: f0[0] = detu_x - 1.0;
230: }
232: void g1_pu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
233: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
234: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
235: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
236: {
237: PetscReal cofu_x[9/*Ncomp*dim*/];
238: PetscInt compI, d;
240: Cof3D(cofu_x, u_x);
241: for (compI = 0; compI < dim; ++compI)
242: for (d = 0; d < dim; ++d) g1[compI*dim+d] = cofu_x[compI*dim+d];
243: }
245: void g2_up_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
246: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
247: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
248: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
249: {
250: PetscReal cofu_x[9/*Ncomp*dim*/];
251: PetscInt compI, d;
253: Cof3D(cofu_x, u_x);
254: for (compI = 0; compI < dim; ++compI)
255: for (d = 0; d < dim; ++d) g2[compI*dim+d] = cofu_x[compI*dim+d];
256: }
258: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
259: {
260: const char *runTypes[2] = {"full", "test"};
261: PetscInt run;
265: options->runType = RUN_FULL;
266: options->mu = 1.0;
267: options->p_wall = 0.4;
269: PetscOptionsBegin(comm, "", "Nonlinear elasticity problem options", "DMPLEX");
270: run = options->runType;
271: PetscOptionsEList("-run_type", "The run type", "ex77.c", runTypes, 2, runTypes[options->runType], &run, NULL);
272: options->runType = (RunType) run;
273: PetscOptionsReal("-shear_modulus", "The shear modulus", "ex77.c", options->mu, &options->mu, NULL);
274: PetscOptionsReal("-wall_pressure", "The wall pressure", "ex77.c", options->p_wall, &options->p_wall, NULL);
275: PetscOptionsEnd();
276: return(0);
277: }
279: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
280: {
284: /* TODO The P1 coordinate space gives wrong results when compared to the affine version. Track this down */
285: if (0) {
286: DMPlexCreateBoxMesh(comm, 3, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);
287: } else {
288: DMCreate(comm, dm);
289: DMSetType(*dm, DMPLEX);
290: }
291: DMSetFromOptions(*dm);
292: /* Label the faces (bit of a hack here, until it is properly implemented for simplices) */
293: DMViewFromOptions(*dm, NULL, "-orig_dm_view");
294: {
295: DM cdm;
296: DMLabel label;
297: IS is;
298: PetscInt d, dim, b, f, Nf;
299: const PetscInt *faces;
300: PetscInt csize;
301: PetscScalar *coords = NULL;
302: PetscSection cs;
303: Vec coordinates ;
305: DMGetDimension(*dm, &dim);
306: DMCreateLabel(*dm, "boundary");
307: DMGetLabel(*dm, "boundary", &label);
308: DMPlexMarkBoundaryFaces(*dm, 1, label);
309: DMGetStratumIS(*dm, "boundary", 1, &is);
310: if (is) {
311: PetscReal faceCoord;
312: PetscInt v;
314: ISGetLocalSize(is, &Nf);
315: ISGetIndices(is, &faces);
317: DMGetCoordinatesLocal(*dm, &coordinates);
318: DMGetCoordinateDM(*dm, &cdm);
319: DMGetLocalSection(cdm, &cs);
321: /* Check for each boundary face if any component of its centroid is either 0.0 or 1.0 */
322: for (f = 0; f < Nf; ++f) {
323: DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
324: /* Calculate mean coordinate vector */
325: for (d = 0; d < dim; ++d) {
326: const PetscInt Nv = csize/dim;
327: faceCoord = 0.0;
328: for (v = 0; v < Nv; ++v) faceCoord += PetscRealPart(coords[v*dim+d]);
329: faceCoord /= Nv;
330: for (b = 0; b < 2; ++b) {
331: if (PetscAbs(faceCoord - b*1.0) < PETSC_SMALL) {
332: DMSetLabelValue(*dm, "Faces", faces[f], d*2+b+1);
333: }
334: }
335: }
336: DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
337: }
338: ISRestoreIndices(is, &faces);
339: }
340: ISDestroy(&is);
341: }
342: DMViewFromOptions(*dm, NULL, "-dm_view");
343: return(0);
344: }
346: PetscErrorCode SetupProblem(DM dm, PetscInt dim, AppCtx *user)
347: {
348: PetscDS ds;
349: PetscWeakForm wf;
350: DMLabel label;
351: PetscInt bd;
355: DMGetDS(dm, &ds);
356: PetscDSSetResidual(ds, 0, NULL, f1_u_3d);
357: PetscDSSetResidual(ds, 1, f0_p_3d, NULL);
358: PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu_3d);
359: PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up_3d, NULL);
360: PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu_3d, NULL, NULL);
362: DMGetLabel(dm, "Faces", &label);
363: DMAddBoundary(dm, DM_BC_NATURAL, "pressure", label, 0, NULL, 0, 0, NULL, NULL, NULL, user, &bd);
364: PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
365: PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_u_3d, 0, NULL);
366: PetscWeakFormSetIndexBdJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, g1_bd_uu_3d, 0, NULL, 0, NULL);
368: DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed", label, 0, NULL, 0, 0, NULL, (void (*)(void)) coordinates, NULL, user, NULL);
369: return(0);
370: }
372: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
373: {
374: PetscErrorCode (*matFuncs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {elasticityMaterial, wallPressure};
375: Vec A;
376: void *ctxs[2];
380: ctxs[0] = user; ctxs[1] = user;
381: DMCreateLocalVector(dmAux, &A);
382: DMProjectFunctionLocal(dmAux, 0.0, matFuncs, ctxs, INSERT_ALL_VALUES, A);
383: DMSetAuxiliaryVec(dm, NULL, 0, A);
384: VecDestroy(&A);
385: return(0);
386: }
388: PetscErrorCode SetupNearNullSpace(DM dm, AppCtx *user)
389: {
390: /* Set up the near null space (a.k.a. rigid body modes) that will be used by the multigrid preconditioner */
391: DM subdm;
392: MatNullSpace nearNullSpace;
393: PetscInt fields = 0;
394: PetscObject deformation;
398: DMCreateSubDM(dm, 1, &fields, NULL, &subdm);
399: DMPlexCreateRigidBody(subdm, 0, &nearNullSpace);
400: DMGetField(dm, 0, NULL, &deformation);
401: PetscObjectCompose(deformation, "nearnullspace", (PetscObject) nearNullSpace);
402: DMDestroy(&subdm);
403: MatNullSpaceDestroy(&nearNullSpace);
404: return(0);
405: }
407: static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
408: {
409: DM dmAux, coordDM;
410: PetscInt f;
414: /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
415: DMGetCoordinateDM(dm, &coordDM);
416: if (!feAux) return(0);
417: DMClone(dm, &dmAux);
418: DMSetCoordinateDM(dmAux, coordDM);
419: for (f = 0; f < NfAux; ++f) {DMSetField(dmAux, f, NULL, (PetscObject) feAux[f]);}
420: DMCreateDS(dmAux);
421: SetupMaterial(dm, dmAux, user);
422: DMDestroy(&dmAux);
423: return(0);
424: }
426: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
427: {
428: DM cdm = dm;
429: PetscFE fe[2], feAux[2];
430: PetscBool simplex;
431: PetscInt dim;
432: MPI_Comm comm;
433: PetscErrorCode ierr;
436: PetscObjectGetComm((PetscObject) dm, &comm);
437: DMGetDimension(dm, &dim);
438: DMPlexIsSimplex(dm, &simplex);
439: /* Create finite element */
440: PetscFECreateDefault(comm, dim, dim, simplex, "def_", PETSC_DEFAULT, &fe[0]);
441: PetscObjectSetName((PetscObject) fe[0], "deformation");
442: PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);
443: PetscFECopyQuadrature(fe[0], fe[1]);
445: PetscObjectSetName((PetscObject) fe[1], "pressure");
447: PetscFECreateDefault(comm, dim, 1, simplex, "elastMat_", PETSC_DEFAULT, &feAux[0]);
448: PetscObjectSetName((PetscObject) feAux[0], "elasticityMaterial");
449: PetscFECopyQuadrature(fe[0], feAux[0]);
450: /* It is not yet possible to define a field on a submesh (e.g. a boundary), so we will use a normal finite element */
451: PetscFECreateDefault(comm, dim, 1, simplex, "wall_pres_", PETSC_DEFAULT, &feAux[1]);
452: PetscObjectSetName((PetscObject) feAux[1], "wall_pressure");
453: PetscFECopyQuadrature(fe[0], feAux[1]);
455: /* Set discretization and boundary conditions for each mesh */
456: DMSetField(dm, 0, NULL, (PetscObject) fe[0]);
457: DMSetField(dm, 1, NULL, (PetscObject) fe[1]);
458: DMCreateDS(dm);
459: SetupProblem(dm, dim, user);
460: while (cdm) {
461: SetupAuxDM(cdm, 2, feAux, user);
462: DMCopyDisc(dm, cdm);
463: DMGetCoarseDM(cdm, &cdm);
464: }
465: PetscFEDestroy(&fe[0]);
466: PetscFEDestroy(&fe[1]);
467: PetscFEDestroy(&feAux[0]);
468: PetscFEDestroy(&feAux[1]);
469: return(0);
470: }
472: int main(int argc, char **argv)
473: {
474: SNES snes; /* nonlinear solver */
475: DM dm; /* problem definition */
476: Vec u,r; /* solution, residual vectors */
477: Mat A,J; /* Jacobian matrix */
478: AppCtx user; /* user-defined work context */
479: PetscInt its; /* iterations for convergence */
482: PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
483: ProcessOptions(PETSC_COMM_WORLD, &user);
484: SNESCreate(PETSC_COMM_WORLD, &snes);
485: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
486: SNESSetDM(snes, dm);
487: DMSetApplicationContext(dm, &user);
489: SetupDiscretization(dm, &user);
490: DMPlexCreateClosureIndex(dm, NULL);
491: SetupNearNullSpace(dm, &user);
493: DMCreateGlobalVector(dm, &u);
494: VecDuplicate(u, &r);
496: DMSetMatType(dm,MATAIJ);
497: DMCreateMatrix(dm, &J);
498: A = J;
500: DMPlexSetSNESLocalFEM(dm,&user,&user,&user);
501: SNESSetJacobian(snes, A, J, NULL, NULL);
503: SNESSetFromOptions(snes);
505: {
506: PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void* ctx);
507: initialGuess[0] = coordinates; initialGuess[1] = zero_scalar;
508: DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
509: }
511: if (user.runType == RUN_FULL) {
512: SNESSolve(snes, NULL, u);
513: SNESGetIterationNumber(snes, &its);
514: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
515: } else {
516: PetscReal res = 0.0;
518: /* Check initial guess */
519: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
520: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
521: /* Check residual */
522: SNESComputeFunction(snes, u, r);
523: PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
524: VecChop(r, 1.0e-10);
525: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
526: VecNorm(r, NORM_2, &res);
527: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", (double)res);
528: /* Check Jacobian */
529: {
530: Vec b;
532: SNESComputeJacobian(snes, u, A, A);
533: VecDuplicate(u, &b);
534: VecSet(r, 0.0);
535: SNESComputeFunction(snes, r, b);
536: MatMult(A, u, r);
537: VecAXPY(r, 1.0, b);
538: VecDestroy(&b);
539: PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
540: VecChop(r, 1.0e-10);
541: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
542: VecNorm(r, NORM_2, &res);
543: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", (double)res);
544: }
545: }
546: VecViewFromOptions(u, NULL, "-sol_vec_view");
548: if (A != J) {MatDestroy(&A);}
549: MatDestroy(&J);
550: VecDestroy(&u);
551: VecDestroy(&r);
552: SNESDestroy(&snes);
553: DMDestroy(&dm);
554: PetscFinalize();
555: return ierr;
556: }
558: /*TEST
560: build:
561: requires: !complex
563: testset:
564: requires: ctetgen
565: args: -run_type full -dm_plex_dim 3 \
566: -def_petscspace_degree 2 -pres_petscspace_degree 1 -elastMat_petscspace_degree 0 -wall_pres_petscspace_degree 0 \
567: -snes_rtol 1e-05 -snes_monitor_short -snes_converged_reason \
568: -ksp_type fgmres -ksp_rtol 1e-10 -ksp_monitor_short -ksp_converged_reason \
569: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper \
570: -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu \
571: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
573: test:
574: suffix: 0
575: requires: !single
576: args: -dm_refine 2 \
577: -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
579: test:
580: suffix: 1
581: requires: superlu_dist
582: nsize: 2
583: args: -dm_distribute -dm_refine 0 -petscpartitioner_type simple \
584: -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
585: timeoutfactor: 2
587: test:
588: suffix: 4
589: requires: superlu_dist
590: nsize: 2
591: args: -dm_distribute -dm_refine 0 -petscpartitioner_type simple \
592: -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
593: output_file: output/ex77_1.out
595: test:
596: suffix: 2
597: requires: !single
598: args: -dm_refine 2 \
599: -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0
601: test:
602: suffix: 2_par
603: requires: superlu_dist !single
604: nsize: 4
605: args: -dm_distribute -dm_refine 2 -petscpartitioner_type simple \
606: -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0
607: output_file: output/ex77_2.out
609: TEST*/