Actual source code: ex49.c

```  1: static char help[] = "   Solves the compressible plane strain elasticity equations in 2d on the unit domain using Q1 finite elements. \n\
2:    Material properties E (Youngs modulus) and nu (Poisson ratio) may vary as a function of space. \n\
3:    The model utilises boundary conditions which produce compression in the x direction. \n\
4: Options: \n"
5:                      "\
6:      -mx : number of elements in x-direction \n\
7:      -my : number of elements in y-direction \n\
8:      -c_str : structure of the coefficients to use. \n"
9:                      "\
10:           -c_str 0 => isotropic material with constant coefficients. \n\
11:                          Parameters: \n\
12:                              -iso_E  : Youngs modulus \n\
13:                              -iso_nu : Poisson ratio \n\
14:           -c_str 1 => step function in the material properties in x. \n\
15:                          Parameters: \n\
16:                               -step_E0  : Youngs modulus to the left of the step \n\
17:                               -step_nu0 : Poisson ratio to the left of the step \n\
18:                               -step_E1  : Youngs modulus to the right of the step \n\
19:                               -step_n1  : Poisson ratio to the right of the step \n\
20:                               -step_xc  : x coordinate of the step \n"
21:                      "\
22:           -c_str 2 => checkerboard material with alternating properties. \n\
23:                       Repeats the following pattern throughout the domain. For example with 4 materials specified, we would heve \n\
24:                       -------------------------\n\
25:                       |  D  |  A  |  B  |  C  |\n\
26:                       ------|-----|-----|------\n\
27:                       |  C  |  D  |  A  |  B  |\n\
28:                       ------|-----|-----|------\n\
29:                       |  B  |  C  |  D  |  A  |\n\
30:                       ------|-----|-----|------\n\
31:                       |  A  |  B  |  C  |  D  |\n\
32:                       -------------------------\n\
33:                       \n\
34:                          Parameters: \n\
35:                               -brick_E    : a comma separated list of Young's modulii \n\
36:                               -brick_nu   : a comma separated list of Poisson ratios  \n\
37:                               -brick_span : the number of elements in x and y each brick will span \n\
38:           -c_str 3 => sponge-like material with alternating properties. \n\
39:                       Repeats the following pattern throughout the domain \n"
40:                      "\
41:                       -----------------------------\n\
42:                       |       [background]        |\n\
43:                       |          E0,nu0           |\n\
44:                       |     -----------------     |\n\
45:                       |     |  [inclusion]  |     |\n\
46:                       |     |    E1,nu1     |     |\n\
47:                       |     |               |     |\n\
48:                       |     | <---- w ----> |     |\n\
49:                       |     |               |     |\n\
50:                       |     |               |     |\n\
51:                       |     -----------------     |\n\
52:                       |                           |\n\
53:                       |                           |\n\
54:                       -----------------------------\n\
55:                       <--------  t + w + t ------->\n\
56:                       \n\
57:                          Parameters: \n\
58:                               -sponge_E0  : Youngs modulus of the surrounding material \n\
59:                               -sponge_E1  : Youngs modulus of the inclusion \n\
60:                               -sponge_nu0 : Poisson ratio of the surrounding material \n\
61:                               -sponge_nu1 : Poisson ratio of the inclusion \n\
62:                               -sponge_t   : the number of elements defining the border around each inclusion \n\
63:                               -sponge_w   : the number of elements in x and y each inclusion will span\n\
64:      -use_gp_coords : Evaluate the Youngs modulus, Poisson ratio and the body force at the global coordinates of the quadrature points.\n\
65:      By default, E, nu and the body force are evaluated at the element center and applied as a constant over the entire element.\n\
66:      -use_nonsymbc : Option to use non-symmetric boundary condition imposition. This choice will use less memory.";

68: /* Contributed by Dave May */

70: #include <petscksp.h>
71: #include <petscdm.h>
72: #include <petscdmda.h>

74: static PetscErrorCode DMDABCApplyCompression(DM, Mat, Vec);
75: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da, Mat A, Vec f, IS *dofs, Mat *AA, Vec *ff);

77: #define NSD          2 /* number of spatial dimensions */
78: #define NODES_PER_EL 4 /* nodes per element */
79: #define U_DOFS       2 /* degrees of freedom per displacement node */
80: #define GAUSS_POINTS 4

82: /* cell based evaluation */
83: typedef struct {
84:   PetscScalar E, nu, fx, fy;
85: } Coefficients;

87: /* Gauss point based evaluation 8+4+4+4 = 20 */
88: typedef struct {
89:   PetscScalar gp_coords[2 * GAUSS_POINTS];
90:   PetscScalar E[GAUSS_POINTS];
91:   PetscScalar nu[GAUSS_POINTS];
92:   PetscScalar fx[GAUSS_POINTS];
93:   PetscScalar fy[GAUSS_POINTS];
94: } GaussPointCoefficients;

96: typedef struct {
97:   PetscScalar ux_dof;
98:   PetscScalar uy_dof;
99: } ElasticityDOF;

101: /*

103:  D = E/((1+nu)(1-2nu)) * [ 1-nu   nu        0     ]
104:                          [  nu   1-nu       0     ]
105:                          [  0     0   0.5*(1-2nu) ]

107:  B = [ d_dx   0   ]
108:      [  0    d_dy ]
109:      [ d_dy  d_dx ]

111:  */

113: /* FEM routines */
114: /*
115:  Element: Local basis function ordering
116:  1-----2
117:  |     |
118:  |     |
119:  0-----3
120:  */
121: static void ConstructQ12D_Ni(PetscScalar _xi[], PetscScalar Ni[])
122: {
123:   PetscScalar xi  = _xi[0];
124:   PetscScalar eta = _xi[1];

126:   Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
127:   Ni[1] = 0.25 * (1.0 - xi) * (1.0 + eta);
128:   Ni[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
129:   Ni[3] = 0.25 * (1.0 + xi) * (1.0 - eta);
130: }

132: static void ConstructQ12D_GNi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
133: {
134:   PetscScalar xi  = _xi[0];
135:   PetscScalar eta = _xi[1];

137:   GNi[0][0] = -0.25 * (1.0 - eta);
138:   GNi[0][1] = -0.25 * (1.0 + eta);
139:   GNi[0][2] = 0.25 * (1.0 + eta);
140:   GNi[0][3] = 0.25 * (1.0 - eta);

142:   GNi[1][0] = -0.25 * (1.0 - xi);
143:   GNi[1][1] = 0.25 * (1.0 - xi);
144:   GNi[1][2] = 0.25 * (1.0 + xi);
145:   GNi[1][3] = -0.25 * (1.0 + xi);
146: }

148: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL], PetscScalar GNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
149: {
150:   PetscScalar J00, J01, J10, J11, J;
151:   PetscScalar iJ00, iJ01, iJ10, iJ11;
152:   PetscInt    i;

154:   J00 = J01 = J10 = J11 = 0.0;
155:   for (i = 0; i < NODES_PER_EL; i++) {
156:     PetscScalar cx = coords[2 * i + 0];
157:     PetscScalar cy = coords[2 * i + 1];

159:     J00 = J00 + GNi[0][i] * cx; /* J_xx = dx/dxi */
160:     J01 = J01 + GNi[0][i] * cy; /* J_xy = dy/dxi */
161:     J10 = J10 + GNi[1][i] * cx; /* J_yx = dx/deta */
162:     J11 = J11 + GNi[1][i] * cy; /* J_yy = dy/deta */
163:   }
164:   J = (J00 * J11) - (J01 * J10);

166:   iJ00 = J11 / J;
167:   iJ01 = -J01 / J;
168:   iJ10 = -J10 / J;
169:   iJ11 = J00 / J;

171:   for (i = 0; i < NODES_PER_EL; i++) {
172:     GNx[0][i] = GNi[0][i] * iJ00 + GNi[1][i] * iJ01;
173:     GNx[1][i] = GNi[0][i] * iJ10 + GNi[1][i] * iJ11;
174:   }

176:   if (det_J) *det_J = J;
177: }

179: static void ConstructGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
180: {
181:   *ngp         = 4;
182:   gp_xi[0][0]  = -0.57735026919;
183:   gp_xi[0][1]  = -0.57735026919;
184:   gp_xi[1][0]  = -0.57735026919;
185:   gp_xi[1][1]  = 0.57735026919;
186:   gp_xi[2][0]  = 0.57735026919;
187:   gp_xi[2][1]  = 0.57735026919;
188:   gp_xi[3][0]  = 0.57735026919;
189:   gp_xi[3][1]  = -0.57735026919;
190:   gp_weight[0] = 1.0;
191:   gp_weight[1] = 1.0;
192:   gp_weight[2] = 1.0;
193:   gp_weight[3] = 1.0;
194: }

196: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da, PetscInt **_lx, PetscInt **_ly)
197: {
198:   PetscMPIInt  rank;
199:   PetscInt     proc_I, proc_J;
200:   PetscInt     cpu_x, cpu_y;
201:   PetscInt     local_mx, local_my;
202:   Vec          vlx, vly;
203:   PetscInt    *LX, *LY, i;
204:   PetscScalar *_a;
205:   Vec          V_SEQ;
206:   VecScatter   ctx;

208:   PetscFunctionBeginUser;
209:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

211:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));

213:   proc_J = rank / cpu_x;
214:   proc_I = rank - cpu_x * proc_J;

216:   PetscCall(PetscMalloc1(cpu_x, &LX));
217:   PetscCall(PetscMalloc1(cpu_y, &LY));

219:   PetscCall(DMDAGetElementsSizes(da, &local_mx, &local_my, NULL));
220:   PetscCall(VecCreate(PETSC_COMM_WORLD, &vlx));
221:   PetscCall(VecSetSizes(vlx, PETSC_DECIDE, cpu_x));
222:   PetscCall(VecSetFromOptions(vlx));

224:   PetscCall(VecCreate(PETSC_COMM_WORLD, &vly));
225:   PetscCall(VecSetSizes(vly, PETSC_DECIDE, cpu_y));
226:   PetscCall(VecSetFromOptions(vly));

228:   PetscCall(VecSetValue(vlx, proc_I, (PetscScalar)(local_mx + 1.0e-9), INSERT_VALUES));
229:   PetscCall(VecSetValue(vly, proc_J, (PetscScalar)(local_my + 1.0e-9), INSERT_VALUES));
230:   PetscCall(VecAssemblyBegin(vlx));
231:   PetscCall(VecAssemblyEnd(vlx));
232:   PetscCall(VecAssemblyBegin(vly));
233:   PetscCall(VecAssemblyEnd(vly));

235:   PetscCall(VecScatterCreateToAll(vlx, &ctx, &V_SEQ));
236:   PetscCall(VecScatterBegin(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
237:   PetscCall(VecScatterEnd(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
238:   PetscCall(VecGetArray(V_SEQ, &_a));
239:   for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
240:   PetscCall(VecRestoreArray(V_SEQ, &_a));
241:   PetscCall(VecScatterDestroy(&ctx));
242:   PetscCall(VecDestroy(&V_SEQ));

244:   PetscCall(VecScatterCreateToAll(vly, &ctx, &V_SEQ));
245:   PetscCall(VecScatterBegin(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
246:   PetscCall(VecScatterEnd(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
247:   PetscCall(VecGetArray(V_SEQ, &_a));
248:   for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
249:   PetscCall(VecRestoreArray(V_SEQ, &_a));
250:   PetscCall(VecScatterDestroy(&ctx));
251:   PetscCall(VecDestroy(&V_SEQ));

253:   *_lx = LX;
254:   *_ly = LY;

256:   PetscCall(VecDestroy(&vlx));
257:   PetscCall(VecDestroy(&vly));
258:   PetscFunctionReturn(PETSC_SUCCESS);
259: }

261: static PetscErrorCode DMDACoordViewGnuplot2d(DM da, const char prefix[])
262: {
263:   DM           cda;
264:   Vec          coords;
265:   DMDACoor2d **_coords;
266:   PetscInt     si, sj, nx, ny, i, j;
267:   FILE        *fp;
268:   char         fname[PETSC_MAX_PATH_LEN];
269:   PetscMPIInt  rank;

271:   PetscFunctionBeginUser;
272:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
273:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
274:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
275:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
276:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### Element geometry for processor %1.4d ### \n", rank));

278:   PetscCall(DMGetCoordinateDM(da, &cda));
279:   PetscCall(DMGetCoordinatesLocal(da, &coords));
280:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
281:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
282:   for (j = sj; j < sj + ny - 1; j++) {
283:     for (i = si; i < si + nx - 1; i++) {
284:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
285:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i].x), (double)PetscRealPart(_coords[j + 1][i].y)));
286:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i + 1].x), (double)PetscRealPart(_coords[j + 1][i + 1].y)));
287:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i + 1].x), (double)PetscRealPart(_coords[j][i + 1].y)));
288:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n\n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
289:     }
290:   }
291:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

293:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

297: static PetscErrorCode DMDAViewGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
298: {
299:   DM           cda;
300:   Vec          coords, local_fields;
301:   DMDACoor2d **_coords;
302:   FILE        *fp;
303:   char         fname[PETSC_MAX_PATH_LEN];
304:   const char  *field_name;
305:   PetscMPIInt  rank;
306:   PetscInt     si, sj, nx, ny, i, j;
307:   PetscInt     n_dofs, d;
308:   PetscScalar *_fields;

310:   PetscFunctionBeginUser;
311:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
312:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
313:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
314:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");

316:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
317:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
318:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
319:   for (d = 0; d < n_dofs; d++) {
320:     PetscCall(DMDAGetFieldName(da, d, &field_name));
321:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
322:   }
323:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));

325:   PetscCall(DMGetCoordinateDM(da, &cda));
326:   PetscCall(DMGetCoordinatesLocal(da, &coords));
327:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
328:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));

330:   PetscCall(DMCreateLocalVector(da, &local_fields));
331:   PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
332:   PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
333:   PetscCall(VecGetArray(local_fields, &_fields));

335:   for (j = sj; j < sj + ny; j++) {
336:     for (i = si; i < si + nx; i++) {
337:       PetscScalar coord_x, coord_y;
338:       PetscScalar field_d;

340:       coord_x = _coords[j][i].x;
341:       coord_y = _coords[j][i].y;

343:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
344:       for (d = 0; d < n_dofs; d++) {
345:         field_d = _fields[n_dofs * ((i - si) + (j - sj) * (nx)) + d];
346:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e ", (double)PetscRealPart(field_d)));
347:       }
348:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
349:     }
350:   }
351:   PetscCall(VecRestoreArray(local_fields, &_fields));
352:   PetscCall(VecDestroy(&local_fields));

354:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

356:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
357:   PetscFunctionReturn(PETSC_SUCCESS);
358: }

360: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
361: {
362:   DM                       cda;
363:   Vec                      local_fields;
364:   FILE                    *fp;
365:   char                     fname[PETSC_MAX_PATH_LEN];
366:   const char              *field_name;
367:   PetscMPIInt              rank;
368:   PetscInt                 si, sj, nx, ny, i, j, p;
369:   PetscInt                 n_dofs, d;
370:   GaussPointCoefficients **_coefficients;

372:   PetscFunctionBeginUser;
373:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
374:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
375:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
376:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");

378:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
379:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
380:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
381:   for (d = 0; d < n_dofs; d++) {
382:     PetscCall(DMDAGetFieldName(da, d, &field_name));
383:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
384:   }
385:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));

387:   PetscCall(DMGetCoordinateDM(da, &cda));
388:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));

390:   PetscCall(DMCreateLocalVector(da, &local_fields));
391:   PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
392:   PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
393:   PetscCall(DMDAVecGetArray(da, local_fields, &_coefficients));

395:   for (j = sj; j < sj + ny; j++) {
396:     for (i = si; i < si + nx; i++) {
397:       PetscScalar coord_x, coord_y;

399:       for (p = 0; p < GAUSS_POINTS; p++) {
400:         coord_x = _coefficients[j][i].gp_coords[2 * p];
401:         coord_y = _coefficients[j][i].gp_coords[2 * p + 1];

403:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));

405:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e %1.6e %1.6e\n", (double)PetscRealPart(_coefficients[j][i].E[p]), (double)PetscRealPart(_coefficients[j][i].nu[p]), (double)PetscRealPart(_coefficients[j][i].fx[p]),
406:                                (double)PetscRealPart(_coefficients[j][i].fy[p])));
407:       }
408:     }
409:   }
410:   PetscCall(DMDAVecRestoreArray(da, local_fields, &_coefficients));
411:   PetscCall(VecDestroy(&local_fields));

413:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
414:   PetscFunctionReturn(PETSC_SUCCESS);
415: }

417: static void FormStressOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar E[], PetscScalar nu[])
418: {
419:   PetscInt    ngp;
420:   PetscScalar gp_xi[GAUSS_POINTS][2];
421:   PetscScalar gp_weight[GAUSS_POINTS];
422:   PetscInt    p, i, j, k, l;
423:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
424:   PetscScalar J_p;
425:   PetscScalar B[3][U_DOFS * NODES_PER_EL];
426:   PetscScalar prop_E, prop_nu, factor, constit_D[3][3];

428:   /* define quadrature rule */

431:   /* evaluate integral */
432:   for (p = 0; p < ngp; p++) {
433:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
434:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);

436:     for (i = 0; i < NODES_PER_EL; i++) {
437:       PetscScalar d_dx_i = GNx_p[0][i];
438:       PetscScalar d_dy_i = GNx_p[1][i];

440:       B[0][2 * i]     = d_dx_i;
441:       B[0][2 * i + 1] = 0.0;
442:       B[1][2 * i]     = 0.0;
443:       B[1][2 * i + 1] = d_dy_i;
444:       B[2][2 * i]     = d_dy_i;
445:       B[2][2 * i + 1] = d_dx_i;
446:     }

448:     /* form D for the quadrature point */
449:     prop_E          = E[p];
450:     prop_nu         = nu[p];
451:     factor          = prop_E / ((1.0 + prop_nu) * (1.0 - 2.0 * prop_nu));
452:     constit_D[0][0] = 1.0 - prop_nu;
453:     constit_D[0][1] = prop_nu;
454:     constit_D[0][2] = 0.0;
455:     constit_D[1][0] = prop_nu;
456:     constit_D[1][1] = 1.0 - prop_nu;
457:     constit_D[1][2] = 0.0;
458:     constit_D[2][0] = 0.0;
459:     constit_D[2][1] = 0.0;
460:     constit_D[2][2] = 0.5 * (1.0 - 2.0 * prop_nu);
461:     for (i = 0; i < 3; i++) {
462:       for (j = 0; j < 3; j++) constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
463:     }

465:     /* form Bt tildeD B */
466:     /*
467:      Ke_ij = Bt_ik . D_kl . B_lj
468:      = B_ki . D_kl . B_lj
469:      */
470:     for (i = 0; i < 8; i++) {
471:       for (j = 0; j < 8; j++) {
472:         for (k = 0; k < 3; k++) {
473:           for (l = 0; l < 3; l++) Ke[8 * i + j] = Ke[8 * i + j] + B[k][i] * constit_D[k][l] * B[l][j];
474:         }
475:       }
476:     }

478:   } /* end quadrature */
479: }

481: static void FormMomentumRhsQ1(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
482: {
483:   PetscInt    ngp;
484:   PetscScalar gp_xi[GAUSS_POINTS][2];
485:   PetscScalar gp_weight[GAUSS_POINTS];
486:   PetscInt    p, i;
487:   PetscScalar Ni_p[NODES_PER_EL];
488:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
489:   PetscScalar J_p, fac;

491:   /* define quadrature rule */

494:   /* evaluate integral */
495:   for (p = 0; p < ngp; p++) {
496:     ConstructQ12D_Ni(gp_xi[p], Ni_p);
497:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
498:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
499:     fac = gp_weight[p] * J_p;

501:     for (i = 0; i < NODES_PER_EL; i++) {
502:       Fe[NSD * i] += fac * Ni_p[i] * fx[p];
503:       Fe[NSD * i + 1] += fac * Ni_p[i] * fy[p];
504:     }
505:   }
506: }

508: /*
509:  i,j are the element indices
510:  The unknown is a vector quantity.
511:  The s[].c is used to indicate the degree of freedom.
512:  */
513: static PetscErrorCode DMDAGetElementEqnums_u(MatStencil s_u[], PetscInt i, PetscInt j)
514: {
515:   PetscFunctionBeginUser;
516:   /* displacement */
517:   /* node 0 */
518:   s_u[0].i = i;
519:   s_u[0].j = j;
520:   s_u[0].c = 0; /* Ux0 */
521:   s_u[1].i = i;
522:   s_u[1].j = j;
523:   s_u[1].c = 1; /* Uy0 */

525:   /* node 1 */
526:   s_u[2].i = i;
527:   s_u[2].j = j + 1;
528:   s_u[2].c = 0; /* Ux1 */
529:   s_u[3].i = i;
530:   s_u[3].j = j + 1;
531:   s_u[3].c = 1; /* Uy1 */

533:   /* node 2 */
534:   s_u[4].i = i + 1;
535:   s_u[4].j = j + 1;
536:   s_u[4].c = 0; /* Ux2 */
537:   s_u[5].i = i + 1;
538:   s_u[5].j = j + 1;
539:   s_u[5].c = 1; /* Uy2 */

541:   /* node 3 */
542:   s_u[6].i = i + 1;
543:   s_u[6].j = j;
544:   s_u[6].c = 0; /* Ux3 */
545:   s_u[7].i = i + 1;
546:   s_u[7].j = j;
547:   s_u[7].c = 1; /* Uy3 */
548:   PetscFunctionReturn(PETSC_SUCCESS);
549: }

551: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords, PetscInt ei, PetscInt ej, PetscScalar el_coords[])
552: {
553:   PetscFunctionBeginUser;
554:   /* get coords for the element */
555:   el_coords[NSD * 0 + 0] = _coords[ej][ei].x;
556:   el_coords[NSD * 0 + 1] = _coords[ej][ei].y;
557:   el_coords[NSD * 1 + 0] = _coords[ej + 1][ei].x;
558:   el_coords[NSD * 1 + 1] = _coords[ej + 1][ei].y;
559:   el_coords[NSD * 2 + 0] = _coords[ej + 1][ei + 1].x;
560:   el_coords[NSD * 2 + 1] = _coords[ej + 1][ei + 1].y;
561:   el_coords[NSD * 3 + 0] = _coords[ej][ei + 1].x;
562:   el_coords[NSD * 3 + 1] = _coords[ej][ei + 1].y;
563:   PetscFunctionReturn(PETSC_SUCCESS);
564: }

566: static PetscErrorCode AssembleA_Elasticity(Mat A, DM elas_da, DM properties_da, Vec properties)
567: {
568:   DM                       cda;
569:   Vec                      coords;
570:   DMDACoor2d             **_coords;
571:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
572:   PetscInt                 sex, sey, mx, my;
573:   PetscInt                 ei, ej;
574:   PetscScalar              Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
575:   PetscScalar              el_coords[NODES_PER_EL * NSD];
576:   Vec                      local_properties;
577:   GaussPointCoefficients **props;
578:   PetscScalar             *prop_E, *prop_nu;

580:   PetscFunctionBeginUser;
581:   /* setup for coords */
582:   PetscCall(DMGetCoordinateDM(elas_da, &cda));
583:   PetscCall(DMGetCoordinatesLocal(elas_da, &coords));
584:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

586:   /* setup for coefficients */
587:   PetscCall(DMCreateLocalVector(properties_da, &local_properties));
588:   PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
589:   PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
590:   PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));

592:   PetscCall(DMDAGetElementsCorners(elas_da, &sex, &sey, 0));
593:   PetscCall(DMDAGetElementsSizes(elas_da, &mx, &my, 0));
594:   for (ej = sey; ej < sey + my; ej++) {
595:     for (ei = sex; ei < sex + mx; ei++) {
596:       /* get coords for the element */
597:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));

599:       /* get coefficients for the element */
600:       prop_E  = props[ej][ei].E;
601:       prop_nu = props[ej][ei].nu;

603:       /* initialise element stiffness matrix */
604:       PetscCall(PetscMemzero(Ae, sizeof(Ae)));

606:       /* form element stiffness matrix */
607:       FormStressOperatorQ1(Ae, el_coords, prop_E, prop_nu);

609:       /* insert element matrix into global matrix */
610:       PetscCall(DMDAGetElementEqnums_u(u_eqn, ei, ej));
611:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
612:     }
613:   }
614:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
615:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

617:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

619:   PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
620:   PetscCall(VecDestroy(&local_properties));
621:   PetscFunctionReturn(PETSC_SUCCESS);
622: }

624: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F, MatStencil u_eqn[], PetscScalar Fe_u[])
625: {
626:   PetscInt n;

628:   PetscFunctionBeginUser;
629:   for (n = 0; n < 4; n++) {
630:     fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].ux_dof         = fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].ux_dof + Fe_u[2 * n];
631:     fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].uy_dof = fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].uy_dof + Fe_u[2 * n + 1];
632:   }
633:   PetscFunctionReturn(PETSC_SUCCESS);
634: }

636: static PetscErrorCode AssembleF_Elasticity(Vec F, DM elas_da, DM properties_da, Vec properties)
637: {
638:   DM                       cda;
639:   Vec                      coords;
640:   DMDACoor2d             **_coords;
641:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
642:   PetscInt                 sex, sey, mx, my;
643:   PetscInt                 ei, ej;
644:   PetscScalar              Fe[NODES_PER_EL * U_DOFS];
645:   PetscScalar              el_coords[NODES_PER_EL * NSD];
646:   Vec                      local_properties;
647:   GaussPointCoefficients **props;
648:   PetscScalar             *prop_fx, *prop_fy;
649:   Vec                      local_F;
650:   ElasticityDOF          **ff;

652:   PetscFunctionBeginUser;
653:   /* setup for coords */
654:   PetscCall(DMGetCoordinateDM(elas_da, &cda));
655:   PetscCall(DMGetCoordinatesLocal(elas_da, &coords));
656:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

658:   /* setup for coefficients */
659:   PetscCall(DMGetLocalVector(properties_da, &local_properties));
660:   PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
661:   PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
662:   PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));

665:   PetscCall(DMGetLocalVector(elas_da, &local_F));
666:   PetscCall(VecZeroEntries(local_F));
667:   PetscCall(DMDAVecGetArray(elas_da, local_F, &ff));

669:   PetscCall(DMDAGetElementsCorners(elas_da, &sex, &sey, 0));
670:   PetscCall(DMDAGetElementsSizes(elas_da, &mx, &my, 0));
671:   for (ej = sey; ej < sey + my; ej++) {
672:     for (ei = sex; ei < sex + mx; ei++) {
673:       /* get coords for the element */
674:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));

676:       /* get coefficients for the element */
677:       prop_fx = props[ej][ei].fx;
678:       prop_fy = props[ej][ei].fy;

680:       /* initialise element stiffness matrix */
681:       PetscCall(PetscMemzero(Fe, sizeof(Fe)));

683:       /* form element stiffness matrix */
684:       FormMomentumRhsQ1(Fe, el_coords, prop_fx, prop_fy);

686:       /* insert element matrix into global matrix */
687:       PetscCall(DMDAGetElementEqnums_u(u_eqn, ei, ej));

690:     }
691:   }

693:   PetscCall(DMDAVecRestoreArray(elas_da, local_F, &ff));
696:   PetscCall(DMRestoreLocalVector(elas_da, &local_F));

698:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

700:   PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
701:   PetscCall(DMRestoreLocalVector(properties_da, &local_properties));
702:   PetscFunctionReturn(PETSC_SUCCESS);
703: }

705: static PetscErrorCode solve_elasticity_2d(PetscInt mx, PetscInt my)
706: {
707:   DM                       elas_da, da_prop;
708:   PetscInt                 u_dof, dof, stencil_width;
709:   Mat                      A;
710:   PetscInt                 mxl, myl;
711:   DM                       prop_cda, vel_cda;
712:   Vec                      prop_coords, vel_coords;
713:   PetscInt                 si, sj, nx, ny, i, j, p;
714:   Vec                      f, X;
715:   PetscInt                 prop_dof, prop_stencil_width;
716:   Vec                      properties, l_properties;
717:   MatNullSpace             matnull;
718:   PetscReal                dx, dy;
719:   PetscInt                 M, N;
720:   DMDACoor2d             **_prop_coords, **_vel_coords;
721:   GaussPointCoefficients **element_props;
722:   KSP                      ksp_E;
723:   PetscInt                 coefficient_structure = 0;
724:   PetscInt                 cpu_x, cpu_y, *lx = NULL, *ly = NULL;
725:   PetscBool                use_gp_coords = PETSC_FALSE;
726:   PetscBool                use_nonsymbc  = PETSC_FALSE;
727:   PetscBool                no_view       = PETSC_FALSE;
728:   PetscBool                flg;

730:   PetscFunctionBeginUser;
731:   /* Generate the da for velocity and pressure */
732:   /*
733:    We use Q1 elements for the temperature.
734:    FEM has a 9-point stencil (BOX) or connectivity pattern
735:    Num nodes in each direction is mx+1, my+1
736:    */
737:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
738:   dof           = u_dof;
739:   stencil_width = 1;
740:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, &elas_da));

742:   PetscCall(DMSetMatType(elas_da, MATAIJ));
743:   PetscCall(DMSetFromOptions(elas_da));
744:   PetscCall(DMSetUp(elas_da));

746:   PetscCall(DMDASetFieldName(elas_da, 0, "Ux"));
747:   PetscCall(DMDASetFieldName(elas_da, 1, "Uy"));

749:   /* unit box [0,1] x [0,1] */
750:   PetscCall(DMDASetUniformCoordinates(elas_da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));

752:   /* Generate element properties, we will assume all material properties are constant over the element */
753:   /* local number of elements */
754:   PetscCall(DMDAGetElementsSizes(elas_da, &mxl, &myl, NULL));

756:   /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
757:   PetscCall(DMDAGetInfo(elas_da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
758:   PetscCall(DMDAGetElementOwnershipRanges2d(elas_da, &lx, &ly));

760:   prop_dof           = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
761:   prop_stencil_width = 0;
762:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, cpu_x, cpu_y, prop_dof, prop_stencil_width, lx, ly, &da_prop));
763:   PetscCall(DMSetFromOptions(da_prop));
764:   PetscCall(DMSetUp(da_prop));

766:   PetscCall(PetscFree(lx));
767:   PetscCall(PetscFree(ly));

769:   /* define centroid positions */
770:   PetscCall(DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
771:   dx = 1.0 / ((PetscReal)(M));
772:   dy = 1.0 / ((PetscReal)(N));

774:   PetscCall(DMDASetUniformCoordinates(da_prop, 0.0 + 0.5 * dx, 1.0 - 0.5 * dx, 0.0 + 0.5 * dy, 1.0 - 0.5 * dy, 0.0, 1.0));

776:   /* define coefficients */
777:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL));

779:   PetscCall(DMCreateGlobalVector(da_prop, &properties));
780:   PetscCall(DMCreateLocalVector(da_prop, &l_properties));
781:   PetscCall(DMDAVecGetArray(da_prop, l_properties, &element_props));

783:   PetscCall(DMGetCoordinateDM(da_prop, &prop_cda));
784:   PetscCall(DMGetCoordinatesLocal(da_prop, &prop_coords));
785:   PetscCall(DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords));

787:   PetscCall(DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0));

789:   PetscCall(DMGetCoordinateDM(elas_da, &vel_cda));
790:   PetscCall(DMGetCoordinatesLocal(elas_da, &vel_coords));
791:   PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));

793:   /* interpolate the coordinates */
794:   for (j = sj; j < sj + ny; j++) {
795:     for (i = si; i < si + nx; i++) {
796:       PetscInt    ngp;
797:       PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
798:       PetscScalar el_coords[8];

800:       PetscCall(GetElementCoords(_vel_coords, i, j, el_coords));

803:       for (p = 0; p < GAUSS_POINTS; p++) {
804:         PetscScalar gp_x, gp_y;
805:         PetscInt    n;
806:         PetscScalar xi_p[2], Ni_p[4];

808:         xi_p[0] = gp_xi[p][0];
809:         xi_p[1] = gp_xi[p][1];
810:         ConstructQ12D_Ni(xi_p, Ni_p);

812:         gp_x = 0.0;
813:         gp_y = 0.0;
814:         for (n = 0; n < NODES_PER_EL; n++) {
815:           gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
816:           gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
817:         }
818:         element_props[j][i].gp_coords[2 * p]     = gp_x;
819:         element_props[j][i].gp_coords[2 * p + 1] = gp_y;
820:       }
821:     }
822:   }

824:   /* define the coefficients */
825:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, &flg));

827:   for (j = sj; j < sj + ny; j++) {
828:     for (i = si; i < si + nx; i++) {
829:       PetscScalar              centroid_x = _prop_coords[j][i].x; /* centroids of cell */
830:       PetscScalar              centroid_y = _prop_coords[j][i].y;
831:       PETSC_UNUSED PetscScalar coord_x, coord_y;

833:       if (coefficient_structure == 0) { /* isotropic */
834:         PetscScalar opts_E, opts_nu;

836:         opts_E  = 1.0;
837:         opts_nu = 0.33;
838:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-iso_E", &opts_E, &flg));
839:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-iso_nu", &opts_nu, &flg));

841:         for (p = 0; p < GAUSS_POINTS; p++) {
842:           element_props[j][i].E[p]  = opts_E;
843:           element_props[j][i].nu[p] = opts_nu;

845:           element_props[j][i].fx[p] = 0.0;
846:           element_props[j][i].fy[p] = 0.0;
847:         }
848:       } else if (coefficient_structure == 1) { /* step */
849:         PetscScalar opts_E0, opts_nu0, opts_xc;
850:         PetscScalar opts_E1, opts_nu1;

852:         opts_E0 = opts_E1 = 1.0;
853:         opts_nu0 = opts_nu1 = 0.333;
854:         opts_xc             = 0.5;
855:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_E0", &opts_E0, &flg));
856:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_nu0", &opts_nu0, &flg));
857:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_E1", &opts_E1, &flg));
858:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_nu1", &opts_nu1, &flg));
859:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_xc", &opts_xc, &flg));

861:         for (p = 0; p < GAUSS_POINTS; p++) {
862:           coord_x = centroid_x;
863:           coord_y = centroid_y;
864:           if (use_gp_coords) {
865:             coord_x = element_props[j][i].gp_coords[2 * p];
866:             coord_y = element_props[j][i].gp_coords[2 * p + 1];
867:           }

869:           element_props[j][i].E[p]  = opts_E0;
870:           element_props[j][i].nu[p] = opts_nu0;
871:           if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
872:             element_props[j][i].E[p]  = opts_E1;
873:             element_props[j][i].nu[p] = opts_nu1;
874:           }

876:           element_props[j][i].fx[p] = 0.0;
877:           element_props[j][i].fy[p] = 0.0;
878:         }
879:       } else if (coefficient_structure == 2) { /* brick */
880:         PetscReal values_E[10];
881:         PetscReal values_nu[10];
882:         PetscInt  nbricks, maxnbricks;
883:         PetscInt  index, span;
884:         PetscInt  jj;

886:         flg        = PETSC_FALSE;
887:         maxnbricks = 10;
888:         PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-brick_E", values_E, &maxnbricks, &flg));
889:         nbricks = maxnbricks;
890:         PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply a list of E values for each brick");

892:         flg        = PETSC_FALSE;
893:         maxnbricks = 10;
894:         PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-brick_nu", values_nu, &maxnbricks, &flg));
895:         PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply a list of nu values for each brick");
896:         PetscCheck(maxnbricks == nbricks, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply equal numbers of values for E and nu");

898:         span = 1;
899:         PetscCall(PetscOptionsGetInt(NULL, NULL, "-brick_span", &span, &flg));

901:         /* cycle through the indices so that no two material properties are repeated in lines of x or y */
902:         jj    = (j / span) % nbricks;
903:         index = (jj + i / span) % nbricks;
904:         /*printf("j=%d: index = %d \n", j,index); */

906:         for (p = 0; p < GAUSS_POINTS; p++) {
907:           element_props[j][i].E[p]  = values_E[index];
908:           element_props[j][i].nu[p] = values_nu[index];
909:         }
910:       } else if (coefficient_structure == 3) { /* sponge */
911:         PetscScalar opts_E0, opts_nu0;
912:         PetscScalar opts_E1, opts_nu1;
913:         PetscInt    opts_t, opts_w;
914:         PetscInt    ii, jj, ci, cj;

916:         opts_E0 = opts_E1 = 1.0;
917:         opts_nu0 = opts_nu1 = 0.333;
918:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_E0", &opts_E0, &flg));
919:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_nu0", &opts_nu0, &flg));
920:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_E1", &opts_E1, &flg));
921:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_nu1", &opts_nu1, &flg));

923:         opts_t = opts_w = 1;
924:         PetscCall(PetscOptionsGetInt(NULL, NULL, "-sponge_t", &opts_t, &flg));
925:         PetscCall(PetscOptionsGetInt(NULL, NULL, "-sponge_w", &opts_w, &flg));

927:         ii = (i) / (opts_t + opts_w + opts_t);
928:         jj = (j) / (opts_t + opts_w + opts_t);

930:         ci = i - ii * (opts_t + opts_w + opts_t);
931:         cj = j - jj * (opts_t + opts_w + opts_t);

933:         for (p = 0; p < GAUSS_POINTS; p++) {
934:           element_props[j][i].E[p]  = opts_E0;
935:           element_props[j][i].nu[p] = opts_nu0;
936:         }
937:         if ((ci >= opts_t) && (ci < opts_t + opts_w)) {
938:           if ((cj >= opts_t) && (cj < opts_t + opts_w)) {
939:             for (p = 0; p < GAUSS_POINTS; p++) {
940:               element_props[j][i].E[p]  = opts_E1;
941:               element_props[j][i].nu[p] = opts_nu1;
942:             }
943:           }
944:         }
945:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
946:     }
947:   }
948:   PetscCall(DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords));

950:   PetscCall(DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords));

952:   PetscCall(DMDAVecRestoreArray(da_prop, l_properties, &element_props));

956:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL));
957:   if (!no_view) {
958:     PetscCall(DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for elasticity eqn.", "properties"));
959:     PetscCall(DMDACoordViewGnuplot2d(elas_da, "mesh"));
960:   }

962:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
963:   PetscCall(DMCreateMatrix(elas_da, &A));
964:   PetscCall(DMGetCoordinates(elas_da, &vel_coords));
965:   PetscCall(MatNullSpaceCreateRigidBody(vel_coords, &matnull));
966:   PetscCall(MatSetNearNullSpace(A, matnull));
967:   PetscCall(MatNullSpaceDestroy(&matnull));
968:   PetscCall(MatCreateVecs(A, &f, &X));

970:   /* assemble A11 */
971:   PetscCall(MatZeroEntries(A));
972:   PetscCall(VecZeroEntries(f));

974:   PetscCall(AssembleA_Elasticity(A, elas_da, da_prop, properties));
975:   /* build force vector */
976:   PetscCall(AssembleF_Elasticity(f, elas_da, da_prop, properties));

978:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_E));
979:   PetscCall(KSPSetOptionsPrefix(ksp_E, "elas_")); /* elasticity */

981:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_nonsymbc", &use_nonsymbc, &flg));
982:   /* solve */
983:   if (!use_nonsymbc) {
984:     Mat        AA;
985:     Vec        ff, XX;
986:     IS         is;
987:     VecScatter scat;

989:     PetscCall(DMDABCApplySymmetricCompression(elas_da, A, f, &is, &AA, &ff));
990:     PetscCall(VecDuplicate(ff, &XX));

992:     PetscCall(KSPSetOperators(ksp_E, AA, AA));
993:     PetscCall(KSPSetFromOptions(ksp_E));

995:     PetscCall(KSPSolve(ksp_E, ff, XX));

997:     /* push XX back into X */
998:     PetscCall(DMDABCApplyCompression(elas_da, NULL, X));

1000:     PetscCall(VecScatterCreate(XX, NULL, X, is, &scat));
1001:     PetscCall(VecScatterBegin(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD));
1002:     PetscCall(VecScatterEnd(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD));
1003:     PetscCall(VecScatterDestroy(&scat));

1005:     PetscCall(MatDestroy(&AA));
1006:     PetscCall(VecDestroy(&ff));
1007:     PetscCall(VecDestroy(&XX));
1008:     PetscCall(ISDestroy(&is));
1009:   } else {
1010:     PetscCall(DMDABCApplyCompression(elas_da, A, f));

1012:     PetscCall(KSPSetOperators(ksp_E, A, A));
1013:     PetscCall(KSPSetFromOptions(ksp_E));

1015:     PetscCall(KSPSolve(ksp_E, f, X));
1016:   }

1018:   if (!no_view) PetscCall(DMDAViewGnuplot2d(elas_da, X, "Displacement solution for elasticity eqn.", "X"));
1019:   PetscCall(KSPDestroy(&ksp_E));

1021:   PetscCall(VecDestroy(&X));
1022:   PetscCall(VecDestroy(&f));
1023:   PetscCall(MatDestroy(&A));

1025:   PetscCall(DMDestroy(&elas_da));
1026:   PetscCall(DMDestroy(&da_prop));

1028:   PetscCall(VecDestroy(&properties));
1029:   PetscCall(VecDestroy(&l_properties));
1030:   PetscFunctionReturn(PETSC_SUCCESS);
1031: }

1033: int main(int argc, char **args)
1034: {
1035:   PetscInt mx, my;

1037:   PetscFunctionBeginUser;
1038:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
1039:   mx = my = 10;
1040:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1041:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1042:   PetscCall(solve_elasticity_2d(mx, my));
1043:   PetscCall(PetscFinalize());
1044:   return 0;
1045: }

1047: /* -------------------------- helpers for boundary conditions -------------------------------- */

1049: static PetscErrorCode BCApply_EAST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
1050: {
1051:   DM                     cda;
1052:   Vec                    coords;
1053:   PetscInt               si, sj, nx, ny, i, j;
1054:   PetscInt               M, N;
1055:   DMDACoor2d           **_coords;
1056:   const PetscInt        *g_idx;
1057:   PetscInt              *bc_global_ids;
1058:   PetscScalar           *bc_vals;
1059:   PetscInt               nbcs;
1060:   PetscInt               n_dofs;
1061:   ISLocalToGlobalMapping ltogm;

1063:   PetscFunctionBeginUser;
1064:   /* enforce bc's */
1065:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
1066:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

1068:   PetscCall(DMGetCoordinateDM(da, &cda));
1069:   PetscCall(DMGetCoordinatesLocal(da, &coords));
1070:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1071:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1072:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

1074:   /* --- */

1076:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1077:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

1079:   /* init the entries to -1 so VecSetValues will ignore them */
1080:   for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;

1082:   i = nx - 1;
1083:   for (j = 0; j < ny; j++) {
1084:     PetscInt                 local_id;
1085:     PETSC_UNUSED PetscScalar coordx, coordy;

1087:     local_id = i + j * nx;

1089:     bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];

1091:     coordx = _coords[j + sj][i + si].x;
1092:     coordy = _coords[j + sj][i + si].y;

1094:     bc_vals[j] = bc_val;
1095:   }
1096:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1097:   nbcs = 0;
1098:   if ((si + nx) == (M)) nbcs = ny;

1100:   if (b) {
1101:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1102:     PetscCall(VecAssemblyBegin(b));
1103:     PetscCall(VecAssemblyEnd(b));
1104:   }
1105:   if (A) PetscCall(MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0));

1107:   PetscCall(PetscFree(bc_vals));
1108:   PetscCall(PetscFree(bc_global_ids));

1110:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1111:   PetscFunctionReturn(PETSC_SUCCESS);
1112: }

1114: static PetscErrorCode BCApply_WEST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
1115: {
1116:   DM                     cda;
1117:   Vec                    coords;
1118:   PetscInt               si, sj, nx, ny, i, j;
1119:   PetscInt               M, N;
1120:   DMDACoor2d           **_coords;
1121:   const PetscInt        *g_idx;
1122:   PetscInt              *bc_global_ids;
1123:   PetscScalar           *bc_vals;
1124:   PetscInt               nbcs;
1125:   PetscInt               n_dofs;
1126:   ISLocalToGlobalMapping ltogm;

1128:   PetscFunctionBeginUser;
1129:   /* enforce bc's */
1130:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
1131:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

1133:   PetscCall(DMGetCoordinateDM(da, &cda));
1134:   PetscCall(DMGetCoordinatesLocal(da, &coords));
1135:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1136:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1137:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

1139:   /* --- */

1141:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1142:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

1144:   /* init the entries to -1 so VecSetValues will ignore them */
1145:   for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;

1147:   i = 0;
1148:   for (j = 0; j < ny; j++) {
1149:     PetscInt                 local_id;
1150:     PETSC_UNUSED PetscScalar coordx, coordy;

1152:     local_id = i + j * nx;

1154:     bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];

1156:     coordx = _coords[j + sj][i + si].x;
1157:     coordy = _coords[j + sj][i + si].y;

1159:     bc_vals[j] = bc_val;
1160:   }
1161:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1162:   nbcs = 0;
1163:   if (si == 0) nbcs = ny;

1165:   if (b) {
1166:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1167:     PetscCall(VecAssemblyBegin(b));
1168:     PetscCall(VecAssemblyEnd(b));
1169:   }
1170:   if (A) PetscCall(MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0));

1172:   PetscCall(PetscFree(bc_vals));
1173:   PetscCall(PetscFree(bc_global_ids));

1175:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1176:   PetscFunctionReturn(PETSC_SUCCESS);
1177: }

1179: static PetscErrorCode DMDABCApplyCompression(DM elas_da, Mat A, Vec f)
1180: {
1181:   PetscFunctionBeginUser;
1182:   PetscCall(BCApply_EAST(elas_da, 0, -1.0, A, f));
1183:   PetscCall(BCApply_EAST(elas_da, 1, 0.0, A, f));
1184:   PetscCall(BCApply_WEST(elas_da, 0, 1.0, A, f));
1185:   PetscCall(BCApply_WEST(elas_da, 1, 0.0, A, f));
1186:   PetscFunctionReturn(PETSC_SUCCESS);
1187: }

1189: static PetscErrorCode Orthogonalize(PetscInt n, Vec *vecs)
1190: {
1191:   PetscInt    i, j;
1192:   PetscScalar dot;

1194:   PetscFunctionBegin;
1195:   for (i = 0; i < n; i++) {
1196:     PetscCall(VecNormalize(vecs[i], NULL));
1197:     for (j = i + 1; j < n; j++) {
1198:       PetscCall(VecDot(vecs[i], vecs[j], &dot));
1199:       PetscCall(VecAXPY(vecs[j], -dot, vecs[i]));
1200:     }
1201:   }
1202:   PetscFunctionReturn(PETSC_SUCCESS);
1203: }

1205: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da, Mat A, Vec f, IS *dofs, Mat *AA, Vec *ff)
1206: {
1207:   PetscInt     start, end, m;
1208:   PetscInt    *unconstrained;
1209:   PetscInt     cnt, i;
1210:   Vec          x;
1211:   PetscScalar *_x;
1212:   IS           is;
1213:   VecScatter   scat;

1215:   PetscFunctionBeginUser;
1216:   /* push bc's into f and A */
1217:   PetscCall(VecDuplicate(f, &x));
1218:   PetscCall(BCApply_EAST(elas_da, 0, -1.0, A, x));
1219:   PetscCall(BCApply_EAST(elas_da, 1, 0.0, A, x));
1220:   PetscCall(BCApply_WEST(elas_da, 0, 1.0, A, x));
1221:   PetscCall(BCApply_WEST(elas_da, 1, 0.0, A, x));

1223:   /* define which dofs are not constrained */
1224:   PetscCall(VecGetLocalSize(x, &m));
1225:   PetscCall(PetscMalloc1(m, &unconstrained));
1226:   PetscCall(VecGetOwnershipRange(x, &start, &end));
1227:   PetscCall(VecGetArray(x, &_x));
1228:   cnt = 0;
1229:   for (i = 0; i < m; i += 2) {
1230:     PetscReal val1, val2;

1232:     val1 = PetscRealPart(_x[i]);
1233:     val2 = PetscRealPart(_x[i + 1]);
1234:     if (PetscAbs(val1) < 0.1 && PetscAbs(val2) < 0.1) {
1235:       unconstrained[cnt] = start + i;
1236:       cnt++;
1237:       unconstrained[cnt] = start + i + 1;
1238:       cnt++;
1239:     }
1240:   }
1241:   PetscCall(VecRestoreArray(x, &_x));

1243:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, cnt, unconstrained, PETSC_COPY_VALUES, &is));
1244:   PetscCall(PetscFree(unconstrained));
1245:   PetscCall(ISSetBlockSize(is, 2));

1247:   /* define correction for dirichlet in the rhs */
1248:   PetscCall(MatMult(A, x, f));
1249:   PetscCall(VecScale(f, -1.0));

1251:   /* get new matrix */
1252:   PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, AA));
1253:   /* get new vector */
1254:   PetscCall(MatCreateVecs(*AA, NULL, ff));

1256:   PetscCall(VecScatterCreate(f, is, *ff, NULL, &scat));
1257:   PetscCall(VecScatterBegin(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD));
1258:   PetscCall(VecScatterEnd(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD));

1260:   { /* Constrain near-null space */
1261:     PetscInt     nvecs;
1262:     const Vec   *vecs;
1263:     Vec         *uvecs;
1264:     PetscBool    has_const;
1265:     MatNullSpace mnull, unull;

1267:     PetscCall(MatGetNearNullSpace(A, &mnull));
1268:     PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvecs, &vecs));
1269:     PetscCall(VecDuplicateVecs(*ff, nvecs, &uvecs));
1270:     for (i = 0; i < nvecs; i++) {
1271:       PetscCall(VecScatterBegin(scat, vecs[i], uvecs[i], INSERT_VALUES, SCATTER_FORWARD));
1272:       PetscCall(VecScatterEnd(scat, vecs[i], uvecs[i], INSERT_VALUES, SCATTER_FORWARD));
1273:     }
1274:     PetscCall(Orthogonalize(nvecs, uvecs));
1275:     PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)A), PETSC_FALSE, nvecs, uvecs, &unull));
1276:     PetscCall(MatSetNearNullSpace(*AA, unull));
1277:     PetscCall(MatNullSpaceDestroy(&unull));
1278:     PetscCall(VecDestroyVecs(nvecs, &uvecs));
1279:   }

1281:   PetscCall(VecScatterDestroy(&scat));

1283:   *dofs = is;
1284:   PetscCall(VecDestroy(&x));
1285:   PetscFunctionReturn(PETSC_SUCCESS);
1286: }

1288: /*TEST

1290:    build:
1291:       requires: !complex !single

1293:    test:
1294:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_rtol 5e-3 -elas_ksp_view
1295:       output_file: output/ex49_1.out

1297:    test:
1298:       suffix: 2
1299:       nsize: 4
1300:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type gcr -elas_pc_type asm -elas_sub_pc_type lu -elas_ksp_rtol 5e-3

1302:    test:
1303:       suffix: 3
1304:       nsize: 4
1305:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,10,1000,100 -brick_nu 0.4,0.2,0.3,0.1 -brick_span 3 -elas_pc_type asm -elas_sub_pc_type lu -elas_ksp_rtol 5e-3

1307:    test:
1308:       suffix: 4
1309:       nsize: 4
1310:       args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type unpreconditioned -mx 40 -my 40 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_mg_levels_ksp_type chebyshev -elas_pc_type ml -elas_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -elas_mg_levels_pc_type pbjacobi -elas_mg_levels_ksp_max_it 3 -use_nonsymbc -elas_pc_ml_nullspace user
1311:       requires: ml

1313:    test:
1314:       suffix: 5
1315:       nsize: 3
1316:       args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_pc_type gamg -elas_mg_fine_ksp_type richardson -elas_mg_fine_pc_type jacobi -elas_mg_fine_pc_jacobi_type rowl1 -elas_mg_fine_pc_jacobi_rowl1_scale .25 -elas_mg_levels_ksp_type chebyshev -elas_mg_levels_ksp_max_it 1 -elas_mg_levels_ksp_chebyshev_esteig 0.2,1.1 -elas_mg_levels_pc_type jacobi -elas_pc_gamg_esteig_ksp_type cg

1318:    test:
1319:       suffix: 6
1320:       nsize: 4
1321:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipegcr -elas_pc_type asm -elas_sub_pc_type lu

1323:    test:
1324:       suffix: 7
1325:       nsize: 4
1326:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipegcr -elas_pc_type asm -elas_sub_pc_type ksp -elas_sub_ksp_ksp_type cg -elas_sub_ksp_ksp_max_it 15

1328:    test:
1329:       suffix: 8
1330:       nsize: 4
1331:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipefgmres -elas_pc_type asm -elas_sub_pc_type ksp -elas_sub_ksp_ksp_type cg -elas_sub_ksp_ksp_max_it 15

1333:    test:
1334:       suffix: hypre_nullspace
1335:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
1336:       args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_pc_type hypre -elas_pc_hypre_boomeramg_nodal_coarsen 6 -elas_pc_hypre_boomeramg_vec_interp_variant 3 -elas_pc_hypre_boomeramg_interp_type ext+i -elas_ksp_view

1338:    test:
1339:       nsize: 4
1340:       suffix: bddc
1341:       args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic

1343:    test:
1344:       nsize: 4
1345:       suffix: bddc_unsym
1346:       args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic -use_nonsymbc -elas_pc_bddc_symmetric 0

1348:    test:
1349:       nsize: 4
1350:       suffix: bddc_unsym_deluxe
1351:       args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic -use_nonsymbc -elas_pc_bddc_symmetric 0 -elas_pc_bddc_use_deluxe_scaling -elas_sub_schurs_symmetric 0

1353:    test:
1354:       nsize: 4
1355:       suffix: fetidp_unsym_deluxe
1356:       args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type fetidp -elas_fetidp_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_fetidp_bddc_pc_bddc_monolithic -use_nonsymbc -elas_fetidp_bddc_pc_bddc_use_deluxe_scaling -elas_fetidp_bddc_sub_schurs_symmetric 0 -elas_fetidp_bddc_pc_bddc_deluxe_singlemat

1358:    test:
1359:       nsize: 4
1360:       suffix: bddc_layerjump
1361:       args: -mx 40 -my 40 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_ksp_norm_type natural

1363:    test:
1364:       nsize: 4
1365:       suffix: bddc_subdomainjump
1366:       args: -mx 40 -my 40 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,1000 -brick_nu 0.4,0.2 -brick_span 20 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_pc_is_use_stiffness_scaling -elas_ksp_norm_type natural

1368:    test:
1369:       nsize: 9
1370:       suffix: bddc_subdomainjump_deluxe
1371:       args: -mx 30 -my 30 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,1000 -brick_nu 0.4,0.2 -brick_span 10 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_pc_bddc_use_deluxe_scaling -elas_ksp_norm_type natural -elas_pc_bddc_schur_layers 1
1372: TEST*/
```