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