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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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;
240: PetscFunctionBeginUser;
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: PetscCall(PetscOptionsEList("-run_type", "The run type", "ex77.c", runTypes, 2, runTypes[options->runType], &run, NULL));
247: options->runType = (RunType)run;
248: PetscCall(PetscOptionsReal("-shear_modulus", "The shear modulus", "ex77.c", options->mu, &options->mu, NULL));
249: PetscCall(PetscOptionsReal("-wall_pressure", "The wall pressure", "ex77.c", options->p_wall, &options->p_wall, NULL));
250: PetscOptionsEnd();
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
255: {
256: PetscFunctionBeginUser;
257: PetscCall(DMCreate(comm, dm));
258: PetscCall(DMSetType(*dm, DMPLEX));
259: PetscCall(DMSetFromOptions(*dm));
260: /* Label the faces (bit of a hack here, until it is properly implemented for simplices) */
261: PetscCall(DMViewFromOptions(*dm, NULL, "-orig_dm_view"));
262: {
263: DM cdm;
264: DMLabel label;
265: IS is;
266: PetscInt d, dim, b, f, Nf;
267: const PetscInt *faces;
268: PetscInt csize;
269: PetscScalar *coords = NULL;
270: PetscSection cs;
271: Vec coordinates;
273: PetscCall(DMGetDimension(*dm, &dim));
274: PetscCall(DMCreateLabel(*dm, "boundary"));
275: PetscCall(DMGetLabel(*dm, "boundary", &label));
276: PetscCall(DMPlexMarkBoundaryFaces(*dm, 1, label));
277: PetscCall(DMGetStratumIS(*dm, "boundary", 1, &is));
278: if (is) {
279: PetscReal faceCoord;
280: PetscInt v;
282: PetscCall(ISGetLocalSize(is, &Nf));
283: PetscCall(ISGetIndices(is, &faces));
285: PetscCall(DMGetCoordinatesLocal(*dm, &coordinates));
286: PetscCall(DMGetCoordinateDM(*dm, &cdm));
287: PetscCall(DMGetLocalSection(cdm, &cs));
289: /* Check for each boundary face if any component of its centroid is either 0.0 or 1.0 */
290: for (f = 0; f < Nf; ++f) {
291: PetscCall(DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords));
292: /* Calculate mean coordinate vector */
293: for (d = 0; d < dim; ++d) {
294: const PetscInt Nv = csize / dim;
295: faceCoord = 0.0;
296: for (v = 0; v < Nv; ++v) faceCoord += PetscRealPart(coords[v * dim + d]);
297: faceCoord /= Nv;
298: for (b = 0; b < 2; ++b) {
299: if (PetscAbs(faceCoord - b * 1.0) < PETSC_SMALL) PetscCall(DMSetLabelValue(*dm, "Faces", faces[f], d * 2 + b + 1));
300: }
301: }
302: PetscCall(DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords));
303: }
304: PetscCall(ISRestoreIndices(is, &faces));
305: }
306: PetscCall(ISDestroy(&is));
307: }
308: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: PetscErrorCode SetupProblem(DM dm, PetscInt dim, AppCtx *user)
313: {
314: PetscDS ds;
315: PetscWeakForm wf;
316: DMLabel label;
317: PetscInt bd;
319: PetscFunctionBeginUser;
320: PetscCall(DMGetDS(dm, &ds));
321: PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u_3d));
322: PetscCall(PetscDSSetResidual(ds, 1, f0_p_3d, NULL));
323: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu_3d));
324: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up_3d, NULL));
325: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu_3d, NULL, NULL));
327: PetscCall(DMGetLabel(dm, "Faces", &label));
328: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "pressure", label, 0, NULL, 0, 0, NULL, NULL, NULL, user, &bd));
329: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
330: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_u_3d, 0, NULL));
331: PetscCall(PetscWeakFormSetIndexBdJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, g1_bd_uu_3d, 0, NULL, 0, NULL));
333: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed", label, 0, NULL, 0, 0, NULL, (void (*)(void))coordinates, NULL, user, NULL));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
338: {
339: PetscErrorCode (*matFuncs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {elasticityMaterial, wallPressure};
340: Vec A;
341: void *ctxs[2];
343: PetscFunctionBegin;
344: ctxs[0] = user;
345: ctxs[1] = user;
346: PetscCall(DMCreateLocalVector(dmAux, &A));
347: PetscCall(DMProjectFunctionLocal(dmAux, 0.0, matFuncs, ctxs, INSERT_ALL_VALUES, A));
348: PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, A));
349: PetscCall(VecDestroy(&A));
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: PetscErrorCode SetupNearNullSpace(DM dm, AppCtx *user)
354: {
355: /* Set up the near null space (a.k.a. rigid body modes) that will be used by the multigrid preconditioner */
356: DM subdm;
357: MatNullSpace nearNullSpace;
358: PetscInt fields = 0;
359: PetscObject deformation;
361: PetscFunctionBeginUser;
362: PetscCall(DMCreateSubDM(dm, 1, &fields, NULL, &subdm));
363: PetscCall(DMPlexCreateRigidBody(subdm, 0, &nearNullSpace));
364: PetscCall(DMGetField(dm, 0, NULL, &deformation));
365: PetscCall(PetscObjectCompose(deformation, "nearnullspace", (PetscObject)nearNullSpace));
366: PetscCall(DMDestroy(&subdm));
367: PetscCall(MatNullSpaceDestroy(&nearNullSpace));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
372: {
373: DM dmAux, coordDM;
374: PetscInt f;
376: PetscFunctionBegin;
377: /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
378: PetscCall(DMGetCoordinateDM(dm, &coordDM));
379: if (!feAux) PetscFunctionReturn(PETSC_SUCCESS);
380: PetscCall(DMClone(dm, &dmAux));
381: PetscCall(DMSetCoordinateDM(dmAux, coordDM));
382: for (f = 0; f < NfAux; ++f) PetscCall(DMSetField(dmAux, f, NULL, (PetscObject)feAux[f]));
383: PetscCall(DMCreateDS(dmAux));
384: PetscCall(SetupMaterial(dm, dmAux, user));
385: PetscCall(DMDestroy(&dmAux));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
390: {
391: DM cdm = dm;
392: PetscFE fe[2], feAux[2];
393: PetscBool simplex;
394: PetscInt dim;
395: MPI_Comm comm;
397: PetscFunctionBeginUser;
398: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
399: PetscCall(DMGetDimension(dm, &dim));
400: PetscCall(DMPlexIsSimplex(dm, &simplex));
401: /* Create finite element */
402: PetscCall(PetscFECreateDefault(comm, dim, dim, simplex, "def_", PETSC_DEFAULT, &fe[0]));
403: PetscCall(PetscObjectSetName((PetscObject)fe[0], "deformation"));
404: PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]));
405: PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
407: PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure"));
409: PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "elastMat_", PETSC_DEFAULT, &feAux[0]));
410: PetscCall(PetscObjectSetName((PetscObject)feAux[0], "elasticityMaterial"));
411: PetscCall(PetscFECopyQuadrature(fe[0], feAux[0]));
412: /* It is not yet possible to define a field on a submesh (e.g. a boundary), so we will use a normal finite element */
413: PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "wall_pres_", PETSC_DEFAULT, &feAux[1]));
414: PetscCall(PetscObjectSetName((PetscObject)feAux[1], "wall_pressure"));
415: PetscCall(PetscFECopyQuadrature(fe[0], feAux[1]));
417: /* Set discretization and boundary conditions for each mesh */
418: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
419: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1]));
420: PetscCall(DMCreateDS(dm));
421: PetscCall(SetupProblem(dm, dim, user));
422: while (cdm) {
423: PetscCall(SetupAuxDM(cdm, 2, feAux, user));
424: PetscCall(DMCopyDisc(dm, cdm));
425: PetscCall(DMGetCoarseDM(cdm, &cdm));
426: }
427: PetscCall(PetscFEDestroy(&fe[0]));
428: PetscCall(PetscFEDestroy(&fe[1]));
429: PetscCall(PetscFEDestroy(&feAux[0]));
430: PetscCall(PetscFEDestroy(&feAux[1]));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: int main(int argc, char **argv)
435: {
436: SNES snes; /* nonlinear solver */
437: DM dm; /* problem definition */
438: Vec u, r; /* solution, residual vectors */
439: Mat A, J; /* Jacobian matrix */
440: AppCtx user; /* user-defined work context */
441: PetscInt its; /* iterations for convergence */
443: PetscFunctionBeginUser;
444: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
445: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
446: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
447: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
448: PetscCall(SNESSetDM(snes, dm));
449: PetscCall(DMSetApplicationContext(dm, &user));
451: PetscCall(SetupDiscretization(dm, &user));
452: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
453: PetscCall(SetupNearNullSpace(dm, &user));
455: PetscCall(DMCreateGlobalVector(dm, &u));
456: PetscCall(PetscObjectSetName((PetscObject)u, "u"));
457: PetscCall(VecDuplicate(u, &r));
459: PetscCall(DMSetMatType(dm, MATAIJ));
460: PetscCall(DMCreateMatrix(dm, &J));
461: A = J;
463: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
464: PetscCall(SNESSetJacobian(snes, A, J, NULL, NULL));
466: PetscCall(SNESSetFromOptions(snes));
468: {
469: PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
470: initialGuess[0] = coordinates;
471: initialGuess[1] = zero_scalar;
472: PetscCall(DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u));
473: }
475: if (user.runType == RUN_FULL) {
476: PetscCall(SNESSolve(snes, NULL, u));
477: PetscCall(SNESGetIterationNumber(snes, &its));
478: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
479: } else {
480: PetscReal res = 0.0;
482: /* Check initial guess */
483: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n"));
484: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));
485: /* Check residual */
486: PetscCall(SNESComputeFunction(snes, u, r));
487: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n"));
488: PetscCall(VecFilter(r, 1.0e-10));
489: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
490: PetscCall(VecNorm(r, NORM_2, &res));
491: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", (double)res));
492: /* Check Jacobian */
493: {
494: Vec b;
496: PetscCall(SNESComputeJacobian(snes, u, A, A));
497: PetscCall(VecDuplicate(u, &b));
498: PetscCall(VecSet(r, 0.0));
499: PetscCall(SNESComputeFunction(snes, r, b));
500: PetscCall(MatMult(A, u, r));
501: PetscCall(VecAXPY(r, 1.0, b));
502: PetscCall(VecDestroy(&b));
503: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n"));
504: PetscCall(VecFilter(r, 1.0e-10));
505: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
506: PetscCall(VecNorm(r, NORM_2, &res));
507: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", (double)res));
508: }
509: }
510: PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));
512: if (A != J) PetscCall(MatDestroy(&A));
513: PetscCall(MatDestroy(&J));
514: PetscCall(VecDestroy(&u));
515: PetscCall(VecDestroy(&r));
516: PetscCall(SNESDestroy(&snes));
517: PetscCall(DMDestroy(&dm));
518: PetscCall(PetscFinalize());
519: return 0;
520: }
522: /*TEST
524: build:
525: requires: !complex
527: testset:
528: requires: ctetgen
529: args: -run_type full -dm_plex_dim 3 \
530: -def_petscspace_degree 2 -pres_petscspace_degree 1 -elastMat_petscspace_degree 0 -wall_pres_petscspace_degree 0 \
531: -snes_rtol 1e-05 -snes_monitor_short -snes_converged_reason \
532: -ksp_type fgmres -ksp_rtol 1e-10 -ksp_monitor_short -ksp_converged_reason \
533: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper \
534: -fieldsplit_deformation_ksp_type preonly -fieldsplit_deformation_pc_type lu \
535: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
537: test:
538: suffix: 0
539: requires: !single
540: args: -dm_refine 2 \
541: -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
543: test:
544: suffix: 1
545: requires: superlu_dist
546: nsize: 2
547: args: -dm_refine 0 -petscpartitioner_type simple \
548: -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
549: timeoutfactor: 2
551: test:
552: suffix: 4
553: requires: superlu_dist
554: nsize: 2
555: args: -dm_refine 0 -petscpartitioner_type simple \
556: -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
557: output_file: output/ex77_1.out
559: test:
560: suffix: 2
561: requires: !single
562: args: -dm_refine 2 \
563: -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0
565: test:
566: suffix: 2_par
567: requires: superlu_dist !single
568: nsize: 4
569: args: -dm_refine 2 -petscpartitioner_type simple \
570: -bc_fixed 3,4,5,6 -bc_pressure 2 -wall_pressure 1.0
571: output_file: output/ex77_2.out
573: TEST*/