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)
9: The Lagrangian density (modulo boundary conditions) for this problem is given by
10:
11:   \frac{\mu}{2} (\mathrm{Tr}{C}-3) + J p + \frac{\kappa}{2} (J-1).
12:

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]));

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"));
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"));

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*/