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;
    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, (PetscVoidFn *)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[], PetscCtx ctx) = {elasticityMaterial, wallPressure};
    340:   Vec      A;
    341:   PetscCtx 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, PetscCtx 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*/