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";
4
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}
13
14 Discretization:
15
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.
19
20 Field Data:
21
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
25
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}]
28
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
32
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}]
34
35 Likewise, DMPlexVecSetClosure() takes data partitioned by field, and correctly
36 puts it into the Sieve ordering.
37
38 */
39
40 #include <petscdmplex.h>
41 #include <petscsnes.h>
42 #include <petscds.h>
43
44 typedef enum {
45 RUN_FULL,
46 RUN_TEST
47 } RunType;
48
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;
54
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
61
Det3D(PetscReal * detJ,const PetscScalar J[])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 }
66
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
76
Cof3D(PetscReal C[],const PetscScalar A[])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 }
89
zero_scalar(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,PetscCtx ctx)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 }
95
zero_vector(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,PetscCtx ctx)96 PetscErrorCode zero_vector(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
97 {
98 const PetscInt Ncomp = dim;
99
100 PetscInt comp;
101 for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0.0;
102 return PETSC_SUCCESS;
103 }
104
coordinates(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,PetscCtx ctx)105 PetscErrorCode coordinates(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
106 {
107 const PetscInt Ncomp = dim;
108
109 PetscInt comp;
110 for (comp = 0; comp < Ncomp; ++comp) u[comp] = x[comp];
111 return PETSC_SUCCESS;
112 }
113
elasticityMaterial(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,PetscCtx ctx)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 }
120
wallPressure(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,PetscCtx ctx)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 }
127
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[])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;
134
135 Cof3D(cofu_x, u_x);
136 Det3D(&detu_x, u_x);
137 p += kappa * (detu_x - 1.0);
138
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 }
144
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[])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;
151
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;
157
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;
162
163 for (d1 = 0; d1 < dim; ++d1) {
164 for (d2 = 0; d2 < dim; ++d2) {
165 const PetscReal g = (d1 == d2) ? 1.0 : 0.0;
166
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 }
173
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[])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;
180
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 }
187
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[])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;
194
195 Cof3D(cofu_x, u_x);
196 Det3D(&detu_x, u_x);
197 p /= detu_x;
198
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 }
207
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[])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 }
214
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[])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;
219
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 }
224
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[])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;
229
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 }
234
ProcessOptions(MPI_Comm comm,AppCtx * options)235 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
236 {
237 const char *runTypes[2] = {"full", "test"};
238 PetscInt run;
239
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 }
253
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)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;
272
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;
281
282 PetscCall(ISGetLocalSize(is, &Nf));
283 PetscCall(ISGetIndices(is, &faces));
284
285 PetscCall(DMGetCoordinatesLocal(*dm, &coordinates));
286 PetscCall(DMGetCoordinateDM(*dm, &cdm));
287 PetscCall(DMGetLocalSection(cdm, &cs));
288
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 }
311
SetupProblem(DM dm,PetscInt dim,AppCtx * user)312 PetscErrorCode SetupProblem(DM dm, PetscInt dim, AppCtx *user)
313 {
314 PetscDS ds;
315 PetscWeakForm wf;
316 DMLabel label;
317 PetscInt bd;
318
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));
326
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));
332
333 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed", label, 0, NULL, 0, 0, NULL, (PetscVoidFn *)coordinates, NULL, user, NULL));
334 PetscFunctionReturn(PETSC_SUCCESS);
335 }
336
SetupMaterial(DM dm,DM dmAux,AppCtx * user)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];
342
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 }
352
SetupNearNullSpace(DM dm,AppCtx * user)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;
360
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 }
370
SetupAuxDM(DM dm,PetscInt NfAux,PetscFE feAux[],AppCtx * user)371 static PetscErrorCode SetupAuxDM(DM dm, PetscInt NfAux, PetscFE feAux[], AppCtx *user)
372 {
373 DM dmAux, coordDM;
374 PetscInt f;
375
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 }
388
SetupDiscretization(DM dm,AppCtx * user)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;
396
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]));
406
407 PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure"));
408
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]));
416
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 }
433
main(int argc,char ** argv)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 */
442
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));
450
451 PetscCall(SetupDiscretization(dm, &user));
452 PetscCall(DMPlexCreateClosureIndex(dm, NULL));
453 PetscCall(SetupNearNullSpace(dm, &user));
454
455 PetscCall(DMCreateGlobalVector(dm, &u));
456 PetscCall(PetscObjectSetName((PetscObject)u, "u"));
457 PetscCall(VecDuplicate(u, &r));
458
459 PetscCall(DMSetMatType(dm, MATAIJ));
460 PetscCall(DMCreateMatrix(dm, &J));
461 A = J;
462
463 PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
464 PetscCall(SNESSetJacobian(snes, A, J, NULL, NULL));
465
466 PetscCall(SNESSetFromOptions(snes));
467
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 }
474
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;
481
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;
495
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"));
511
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 }
521
522 /*TEST
523
524 build:
525 requires: !complex
526
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
536
537 test:
538 suffix: 0
539 requires: !single
540 args: -dm_refine 2 \
541 -bc_fixed 1 -bc_pressure 2 -wall_pressure 0.4
542
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
550
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
558
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
564
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
572
573 TEST*/
574