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