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.";
67
68 /* Contributed by Dave May */
69
70 #include <petscksp.h>
71 #include <petscdm.h>
72 #include <petscdmda.h>
73
74 static PetscErrorCode DMDABCApplyCompression(DM, Mat, Vec);
75 static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da, Mat A, Vec f, IS *dofs, Mat *AA, Vec *ff);
76
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
81
82 /* cell based evaluation */
83 typedef struct {
84 PetscScalar E, nu, fx, fy;
85 } Coefficients;
86
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;
95
96 typedef struct {
97 PetscScalar ux_dof;
98 PetscScalar uy_dof;
99 } ElasticityDOF;
100
101 /*
102
103 D = E/((1+nu)(1-2nu)) * [ 1-nu nu 0 ]
104 [ nu 1-nu 0 ]
105 [ 0 0 0.5*(1-2nu) ]
106
107 B = [ d_dx 0 ]
108 [ 0 d_dy ]
109 [ d_dy d_dx ]
110
111 */
112
113 /* FEM routines */
114 /*
115 Element: Local basis function ordering
116 1-----2
117 | |
118 | |
119 0-----3
120 */
ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])121 static void ConstructQ12D_Ni(PetscScalar _xi[], PetscScalar Ni[])
122 {
123 PetscScalar xi = _xi[0];
124 PetscScalar eta = _xi[1];
125
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 }
131
ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])132 static void ConstructQ12D_GNi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
133 {
134 PetscScalar xi = _xi[0];
135 PetscScalar eta = _xi[1];
136
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);
141
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 }
147
ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar * det_J)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;
153
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];
158
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);
165
166 iJ00 = J11 / J;
167 iJ01 = -J01 / J;
168 iJ10 = -J10 / J;
169 iJ11 = J00 / J;
170
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 }
175
176 if (det_J) *det_J = J;
177 }
178
ConstructGaussQuadrature(PetscInt * ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])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 }
195
DMDAGetElementOwnershipRanges2d(DM da,PetscInt ** _lx,PetscInt ** _ly)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;
207
208 PetscFunctionBeginUser;
209 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
210
211 PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
212
213 proc_J = rank / cpu_x;
214 proc_I = rank - cpu_x * proc_J;
215
216 PetscCall(PetscMalloc1(cpu_x, &LX));
217 PetscCall(PetscMalloc1(cpu_y, &LY));
218
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));
223
224 PetscCall(VecCreate(PETSC_COMM_WORLD, &vly));
225 PetscCall(VecSetSizes(vly, PETSC_DECIDE, cpu_y));
226 PetscCall(VecSetFromOptions(vly));
227
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));
234
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));
243
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));
252
253 *_lx = LX;
254 *_ly = LY;
255
256 PetscCall(VecDestroy(&vlx));
257 PetscCall(VecDestroy(&vly));
258 PetscFunctionReturn(PETSC_SUCCESS);
259 }
260
DMDACoordViewGnuplot2d(DM da,const char prefix[])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;
270
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));
277
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));
292
293 PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
294 PetscFunctionReturn(PETSC_SUCCESS);
295 }
296
DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])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;
309
310 PetscFunctionBeginUser;
311 PetscCallMPI(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");
315
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"));
324
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));
329
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));
334
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;
339
340 coord_x = _coords[j][i].x;
341 coord_y = _coords[j][i].y;
342
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));
353
354 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
355
356 PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
357 PetscFunctionReturn(PETSC_SUCCESS);
358 }
359
DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])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;
371
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");
377
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"));
386
387 PetscCall(DMGetCoordinateDM(da, &cda));
388 PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
389
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));
394
395 for (j = sj; j < sj + ny; j++) {
396 for (i = si; i < si + nx; i++) {
397 PetscScalar coord_x, coord_y;
398
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];
402
403 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
404
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));
412
413 PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
414 PetscFunctionReturn(PETSC_SUCCESS);
415 }
416
FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar E[],PetscScalar nu[])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];
427
428 /* define quadrature rule */
429 ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
430
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);
435
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];
439
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 }
447
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 }
464
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 }
477
478 } /* end quadrature */
479 }
480
FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])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;
490
491 /* define quadrature rule */
492 ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
493
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;
500
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 }
507
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 */
DMDAGetElementEqnums_u(MatStencil s_u[],PetscInt i,PetscInt j)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 */
524
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 */
532
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 */
540
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 }
550
GetElementCoords(DMDACoor2d ** _coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])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 }
565
AssembleA_Elasticity(Mat A,DM elas_da,DM properties_da,Vec properties)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;
579
580 PetscFunctionBeginUser;
581 /* setup for coords */
582 PetscCall(DMGetCoordinateDM(elas_da, &cda));
583 PetscCall(DMGetCoordinatesLocal(elas_da, &coords));
584 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
585
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));
591
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));
598
599 /* get coefficients for the element */
600 prop_E = props[ej][ei].E;
601 prop_nu = props[ej][ei].nu;
602
603 /* initialise element stiffness matrix */
604 PetscCall(PetscMemzero(Ae, sizeof(Ae)));
605
606 /* form element stiffness matrix */
607 FormStressOperatorQ1(Ae, el_coords, prop_E, prop_nu);
608
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));
616
617 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
618
619 PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
620 PetscCall(VecDestroy(&local_properties));
621 PetscFunctionReturn(PETSC_SUCCESS);
622 }
623
DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF ** fields_F,MatStencil u_eqn[],PetscScalar Fe_u[])624 static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F, MatStencil u_eqn[], PetscScalar Fe_u[])
625 {
626 PetscInt n;
627
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 }
635
AssembleF_Elasticity(Vec F,DM elas_da,DM properties_da,Vec properties)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;
651
652 PetscFunctionBeginUser;
653 /* setup for coords */
654 PetscCall(DMGetCoordinateDM(elas_da, &cda));
655 PetscCall(DMGetCoordinatesLocal(elas_da, &coords));
656 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
657
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));
663
664 /* get access to the vector */
665 PetscCall(DMGetLocalVector(elas_da, &local_F));
666 PetscCall(VecZeroEntries(local_F));
667 PetscCall(DMDAVecGetArray(elas_da, local_F, &ff));
668
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));
675
676 /* get coefficients for the element */
677 prop_fx = props[ej][ei].fx;
678 prop_fy = props[ej][ei].fy;
679
680 /* initialise element stiffness matrix */
681 PetscCall(PetscMemzero(Fe, sizeof(Fe)));
682
683 /* form element stiffness matrix */
684 FormMomentumRhsQ1(Fe, el_coords, prop_fx, prop_fy);
685
686 /* insert element matrix into global matrix */
687 PetscCall(DMDAGetElementEqnums_u(u_eqn, ei, ej));
688
689 PetscCall(DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, Fe));
690 }
691 }
692
693 PetscCall(DMDAVecRestoreArray(elas_da, local_F, &ff));
694 PetscCall(DMLocalToGlobalBegin(elas_da, local_F, ADD_VALUES, F));
695 PetscCall(DMLocalToGlobalEnd(elas_da, local_F, ADD_VALUES, F));
696 PetscCall(DMRestoreLocalVector(elas_da, &local_F));
697
698 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
699
700 PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
701 PetscCall(DMRestoreLocalVector(properties_da, &local_properties));
702 PetscFunctionReturn(PETSC_SUCCESS);
703 }
704
solve_elasticity_2d(PetscInt mx,PetscInt my)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;
729
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));
741
742 PetscCall(DMSetMatType(elas_da, MATAIJ));
743 PetscCall(DMSetFromOptions(elas_da));
744 PetscCall(DMSetUp(elas_da));
745
746 PetscCall(DMDASetFieldName(elas_da, 0, "Ux"));
747 PetscCall(DMDASetFieldName(elas_da, 1, "Uy"));
748
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));
751
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));
755
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));
759
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));
765
766 PetscCall(PetscFree(lx));
767 PetscCall(PetscFree(ly));
768
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;
773
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));
775
776 /* define coefficients */
777 PetscCall(PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL));
778
779 PetscCall(DMCreateGlobalVector(da_prop, &properties));
780 PetscCall(DMCreateLocalVector(da_prop, &l_properties));
781 PetscCall(DMDAVecGetArray(da_prop, l_properties, &element_props));
782
783 PetscCall(DMGetCoordinateDM(da_prop, &prop_cda));
784 PetscCall(DMGetCoordinatesLocal(da_prop, &prop_coords));
785 PetscCall(DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords));
786
787 PetscCall(DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0));
788
789 PetscCall(DMGetCoordinateDM(elas_da, &vel_cda));
790 PetscCall(DMGetCoordinatesLocal(elas_da, &vel_coords));
791 PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));
792
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];
799
800 PetscCall(GetElementCoords(_vel_coords, i, j, el_coords));
801 ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
802
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];
807
808 xi_p[0] = gp_xi[p][0];
809 xi_p[1] = gp_xi[p][1];
810 ConstructQ12D_Ni(xi_p, Ni_p);
811
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 }
823
824 /* define the coefficients */
825 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, &flg));
826
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;
832
833 if (coefficient_structure == 0) { /* isotropic */
834 PetscScalar opts_E, opts_nu;
835
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));
840
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;
844
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;
851
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));
860
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 }
868
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 }
875
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;
885
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");
891
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");
897
898 span = 1;
899 PetscCall(PetscOptionsGetInt(NULL, NULL, "-brick_span", &span, &flg));
900
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); */
905
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;
915
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));
922
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));
926
927 ii = (i) / (opts_t + opts_w + opts_t);
928 jj = (j) / (opts_t + opts_w + opts_t);
929
930 ci = i - ii * (opts_t + opts_w + opts_t);
931 cj = j - jj * (opts_t + opts_w + opts_t);
932
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));
949
950 PetscCall(DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords));
951
952 PetscCall(DMDAVecRestoreArray(da_prop, l_properties, &element_props));
953 PetscCall(DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties));
954 PetscCall(DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties));
955
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 }
961
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));
969
970 /* assemble A11 */
971 PetscCall(MatZeroEntries(A));
972 PetscCall(VecZeroEntries(f));
973
974 PetscCall(AssembleA_Elasticity(A, elas_da, da_prop, properties));
975 /* build force vector */
976 PetscCall(AssembleF_Elasticity(f, elas_da, da_prop, properties));
977
978 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_E));
979 PetscCall(KSPSetOptionsPrefix(ksp_E, "elas_")); /* elasticity */
980
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;
988
989 PetscCall(DMDABCApplySymmetricCompression(elas_da, A, f, &is, &AA, &ff));
990 PetscCall(VecDuplicate(ff, &XX));
991
992 PetscCall(KSPSetOperators(ksp_E, AA, AA));
993 PetscCall(KSPSetFromOptions(ksp_E));
994
995 PetscCall(KSPSolve(ksp_E, ff, XX));
996
997 /* push XX back into X */
998 PetscCall(DMDABCApplyCompression(elas_da, NULL, X));
999
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));
1004
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));
1011
1012 PetscCall(KSPSetOperators(ksp_E, A, A));
1013 PetscCall(KSPSetFromOptions(ksp_E));
1014
1015 PetscCall(KSPSolve(ksp_E, f, X));
1016 }
1017
1018 if (!no_view) PetscCall(DMDAViewGnuplot2d(elas_da, X, "Displacement solution for elasticity eqn.", "X"));
1019 PetscCall(KSPDestroy(&ksp_E));
1020
1021 PetscCall(VecDestroy(&X));
1022 PetscCall(VecDestroy(&f));
1023 PetscCall(MatDestroy(&A));
1024
1025 PetscCall(DMDestroy(&elas_da));
1026 PetscCall(DMDestroy(&da_prop));
1027
1028 PetscCall(VecDestroy(&properties));
1029 PetscCall(VecDestroy(&l_properties));
1030 PetscFunctionReturn(PETSC_SUCCESS);
1031 }
1032
main(int argc,char ** args)1033 int main(int argc, char **args)
1034 {
1035 PetscInt mx, my;
1036
1037 PetscFunctionBeginUser;
1038 PetscCall(PetscInitialize(&argc, &args, NULL, 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 }
1046
1047 /* -------------------------- helpers for boundary conditions -------------------------------- */
1048
BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)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;
1062
1063 PetscFunctionBeginUser;
1064 /* enforce bc's */
1065 PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1066 PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1067
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));
1073
1074 /* --- */
1075
1076 PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1077 PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1078
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;
1081
1082 i = nx - 1;
1083 for (j = 0; j < ny; j++) {
1084 PetscInt local_id;
1085 PETSC_UNUSED PetscScalar coordx, coordy;
1086
1087 local_id = i + j * nx;
1088
1089 bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1090
1091 coordx = _coords[j + sj][i + si].x;
1092 coordy = _coords[j + sj][i + si].y;
1093
1094 bc_vals[j] = bc_val;
1095 }
1096 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1097 nbcs = 0;
1098 if ((si + nx) == (M)) nbcs = ny;
1099
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));
1106
1107 PetscCall(PetscFree(bc_vals));
1108 PetscCall(PetscFree(bc_global_ids));
1109
1110 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1111 PetscFunctionReturn(PETSC_SUCCESS);
1112 }
1113
BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)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;
1127
1128 PetscFunctionBeginUser;
1129 /* enforce bc's */
1130 PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1131 PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1132
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));
1138
1139 /* --- */
1140
1141 PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1142 PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1143
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;
1146
1147 i = 0;
1148 for (j = 0; j < ny; j++) {
1149 PetscInt local_id;
1150 PETSC_UNUSED PetscScalar coordx, coordy;
1151
1152 local_id = i + j * nx;
1153
1154 bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1155
1156 coordx = _coords[j + sj][i + si].x;
1157 coordy = _coords[j + sj][i + si].y;
1158
1159 bc_vals[j] = bc_val;
1160 }
1161 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1162 nbcs = 0;
1163 if (si == 0) nbcs = ny;
1164
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));
1171
1172 PetscCall(PetscFree(bc_vals));
1173 PetscCall(PetscFree(bc_global_ids));
1174
1175 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1176 PetscFunctionReturn(PETSC_SUCCESS);
1177 }
1178
DMDABCApplyCompression(DM elas_da,Mat A,Vec f)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 }
1188
Orthogonalize(PetscInt n,Vec * vecs)1189 static PetscErrorCode Orthogonalize(PetscInt n, Vec *vecs)
1190 {
1191 PetscInt i, j;
1192 PetscScalar dot;
1193
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 }
1204
DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS * dofs,Mat * AA,Vec * ff)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;
1214
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));
1222
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;
1231
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));
1242
1243 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, cnt, unconstrained, PETSC_COPY_VALUES, &is));
1244 PetscCall(PetscFree(unconstrained));
1245 PetscCall(ISSetBlockSize(is, 2));
1246
1247 /* define correction for dirichlet in the rhs */
1248 PetscCall(MatMult(A, x, f));
1249 PetscCall(VecScale(f, -1.0));
1250
1251 /* get new matrix */
1252 PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, AA));
1253 /* get new vector */
1254 PetscCall(MatCreateVecs(*AA, NULL, ff));
1255
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));
1259
1260 { /* Constrain near-null space */
1261 PetscInt nvecs;
1262 const Vec *vecs;
1263 Vec *uvecs;
1264 PetscBool has_const;
1265 MatNullSpace mnull, unull;
1266
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 }
1280
1281 PetscCall(VecScatterDestroy(&scat));
1282
1283 *dofs = is;
1284 PetscCall(VecDestroy(&x));
1285 PetscFunctionReturn(PETSC_SUCCESS);
1286 }
1287
1288 /*TEST
1289
1290 build:
1291 requires: !complex !single
1292
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
1296
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
1301
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
1306
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
1312
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
1317
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
1322
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
1327
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
1332
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
1337
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
1342
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
1347
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
1352
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
1357
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
1362
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
1367
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*/
1373