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;
308: PetscScalar *_fields;
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");
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 (PetscInt 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 (PetscInt 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;
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 (PetscInt 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 */
429: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
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 */
492: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
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 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 (PetscInt 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: PetscFunctionBeginUser;
627: for (PetscInt n = 0; n < 4; n++) {
628: 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];
629: 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];
630: }
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: static PetscErrorCode AssembleF_Elasticity(Vec F, DM elas_da, DM properties_da, Vec properties)
635: {
636: DM cda;
637: Vec coords;
638: DMDACoor2d **_coords;
639: MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
640: PetscInt sex, sey, mx, my;
641: PetscScalar Fe[NODES_PER_EL * U_DOFS];
642: PetscScalar el_coords[NODES_PER_EL * NSD];
643: Vec local_properties;
644: GaussPointCoefficients **props;
645: PetscScalar *prop_fx, *prop_fy;
646: Vec local_F;
647: ElasticityDOF **ff;
649: PetscFunctionBeginUser;
650: /* setup for coords */
651: PetscCall(DMGetCoordinateDM(elas_da, &cda));
652: PetscCall(DMGetCoordinatesLocal(elas_da, &coords));
653: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
655: /* setup for coefficients */
656: PetscCall(DMGetLocalVector(properties_da, &local_properties));
657: PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
658: PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
659: PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));
661: /* get access to the vector */
662: PetscCall(DMGetLocalVector(elas_da, &local_F));
663: PetscCall(VecZeroEntries(local_F));
664: PetscCall(DMDAVecGetArray(elas_da, local_F, &ff));
666: PetscCall(DMDAGetElementsCorners(elas_da, &sex, &sey, 0));
667: PetscCall(DMDAGetElementsSizes(elas_da, &mx, &my, 0));
668: for (PetscInt ej = sey; ej < sey + my; ej++) {
669: for (PetscInt ei = sex; ei < sex + mx; ei++) {
670: /* get coords for the element */
671: PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
673: /* get coefficients for the element */
674: prop_fx = props[ej][ei].fx;
675: prop_fy = props[ej][ei].fy;
677: /* initialise element stiffness matrix */
678: PetscCall(PetscMemzero(Fe, sizeof(Fe)));
680: /* form element stiffness matrix */
681: FormMomentumRhsQ1(Fe, el_coords, prop_fx, prop_fy);
683: /* insert element matrix into global matrix */
684: PetscCall(DMDAGetElementEqnums_u(u_eqn, ei, ej));
686: PetscCall(DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, Fe));
687: }
688: }
690: PetscCall(DMDAVecRestoreArray(elas_da, local_F, &ff));
691: PetscCall(DMLocalToGlobalBegin(elas_da, local_F, ADD_VALUES, F));
692: PetscCall(DMLocalToGlobalEnd(elas_da, local_F, ADD_VALUES, F));
693: PetscCall(DMRestoreLocalVector(elas_da, &local_F));
695: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
697: PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
698: PetscCall(DMRestoreLocalVector(properties_da, &local_properties));
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }
702: static PetscErrorCode solve_elasticity_2d(PetscInt mx, PetscInt my)
703: {
704: DM elas_da, da_prop;
705: PetscInt u_dof, dof, stencil_width;
706: Mat A;
707: PetscInt mxl, myl;
708: DM prop_cda, vel_cda;
709: Vec prop_coords, vel_coords;
710: PetscInt si, sj, nx, ny, i, j, p;
711: Vec f, X;
712: PetscInt prop_dof, prop_stencil_width;
713: Vec properties, l_properties;
714: MatNullSpace matnull;
715: PetscReal dx, dy;
716: PetscInt M, N;
717: DMDACoor2d **_prop_coords, **_vel_coords;
718: GaussPointCoefficients **element_props;
719: KSP ksp_E;
720: PetscInt coefficient_structure = 0;
721: PetscInt cpu_x, cpu_y, *lx = NULL, *ly = NULL;
722: PetscBool use_gp_coords = PETSC_FALSE;
723: PetscBool use_nonsymbc = PETSC_FALSE;
724: PetscBool no_view = PETSC_FALSE;
725: PetscBool flg;
727: PetscFunctionBeginUser;
728: /* Generate the da for velocity and pressure */
729: /*
730: We use Q1 elements for the temperature.
731: FEM has a 9-point stencil (BOX) or connectivity pattern
732: Num nodes in each direction is mx+1, my+1
733: */
734: u_dof = U_DOFS; /* Vx, Vy - velocities */
735: dof = u_dof;
736: stencil_width = 1;
737: 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));
739: PetscCall(DMSetMatType(elas_da, MATAIJ));
740: PetscCall(DMSetFromOptions(elas_da));
741: PetscCall(DMSetUp(elas_da));
743: PetscCall(DMDASetFieldName(elas_da, 0, "Ux"));
744: PetscCall(DMDASetFieldName(elas_da, 1, "Uy"));
746: /* unit box [0,1] x [0,1] */
747: PetscCall(DMDASetUniformCoordinates(elas_da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
749: /* Generate element properties, we will assume all material properties are constant over the element */
750: /* local number of elements */
751: PetscCall(DMDAGetElementsSizes(elas_da, &mxl, &myl, NULL));
753: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
754: PetscCall(DMDAGetInfo(elas_da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
755: PetscCall(DMDAGetElementOwnershipRanges2d(elas_da, &lx, &ly));
757: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
758: prop_stencil_width = 0;
759: 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));
760: PetscCall(DMSetFromOptions(da_prop));
761: PetscCall(DMSetUp(da_prop));
763: PetscCall(PetscFree(lx));
764: PetscCall(PetscFree(ly));
766: /* define centroid positions */
767: PetscCall(DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
768: dx = 1.0 / (PetscReal)M;
769: dy = 1.0 / (PetscReal)N;
771: 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));
773: /* define coefficients */
774: PetscCall(PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL));
776: PetscCall(DMCreateGlobalVector(da_prop, &properties));
777: PetscCall(DMCreateLocalVector(da_prop, &l_properties));
778: PetscCall(DMDAVecGetArray(da_prop, l_properties, &element_props));
780: PetscCall(DMGetCoordinateDM(da_prop, &prop_cda));
781: PetscCall(DMGetCoordinatesLocal(da_prop, &prop_coords));
782: PetscCall(DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords));
784: PetscCall(DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0));
786: PetscCall(DMGetCoordinateDM(elas_da, &vel_cda));
787: PetscCall(DMGetCoordinatesLocal(elas_da, &vel_coords));
788: PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));
790: /* interpolate the coordinates */
791: for (j = sj; j < sj + ny; j++) {
792: for (i = si; i < si + nx; i++) {
793: PetscInt ngp;
794: PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
795: PetscScalar el_coords[8];
797: PetscCall(GetElementCoords(_vel_coords, i, j, el_coords));
798: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
800: for (p = 0; p < GAUSS_POINTS; p++) {
801: PetscScalar gp_x, gp_y;
802: PetscScalar xi_p[2], Ni_p[4];
804: xi_p[0] = gp_xi[p][0];
805: xi_p[1] = gp_xi[p][1];
806: ConstructQ12D_Ni(xi_p, Ni_p);
808: gp_x = 0.0;
809: gp_y = 0.0;
810: for (PetscInt n = 0; n < NODES_PER_EL; n++) {
811: gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
812: gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
813: }
814: element_props[j][i].gp_coords[2 * p] = gp_x;
815: element_props[j][i].gp_coords[2 * p + 1] = gp_y;
816: }
817: }
818: }
820: /* define the coefficients */
821: PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, &flg));
823: for (j = sj; j < sj + ny; j++) {
824: for (i = si; i < si + nx; i++) {
825: PetscScalar centroid_x = _prop_coords[j][i].x; /* centroids of cell */
826: PetscScalar centroid_y = _prop_coords[j][i].y;
827: PETSC_UNUSED PetscScalar coord_x, coord_y;
829: if (coefficient_structure == 0) { /* isotropic */
830: PetscScalar opts_E, opts_nu;
832: opts_E = 1.0;
833: opts_nu = 0.33;
834: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-iso_E", &opts_E, &flg));
835: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-iso_nu", &opts_nu, &flg));
837: for (p = 0; p < GAUSS_POINTS; p++) {
838: element_props[j][i].E[p] = opts_E;
839: element_props[j][i].nu[p] = opts_nu;
841: element_props[j][i].fx[p] = 0.0;
842: element_props[j][i].fy[p] = 0.0;
843: }
844: } else if (coefficient_structure == 1) { /* step */
845: PetscScalar opts_E0, opts_nu0, opts_xc;
846: PetscScalar opts_E1, opts_nu1;
848: opts_E0 = opts_E1 = 1.0;
849: opts_nu0 = opts_nu1 = 0.333;
850: opts_xc = 0.5;
851: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_E0", &opts_E0, &flg));
852: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_nu0", &opts_nu0, &flg));
853: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_E1", &opts_E1, &flg));
854: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_nu1", &opts_nu1, &flg));
855: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_xc", &opts_xc, &flg));
857: for (p = 0; p < GAUSS_POINTS; p++) {
858: coord_x = centroid_x;
859: coord_y = centroid_y;
860: if (use_gp_coords) {
861: coord_x = element_props[j][i].gp_coords[2 * p];
862: coord_y = element_props[j][i].gp_coords[2 * p + 1];
863: }
865: element_props[j][i].E[p] = opts_E0;
866: element_props[j][i].nu[p] = opts_nu0;
867: if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
868: element_props[j][i].E[p] = opts_E1;
869: element_props[j][i].nu[p] = opts_nu1;
870: }
872: element_props[j][i].fx[p] = 0.0;
873: element_props[j][i].fy[p] = 0.0;
874: }
875: } else if (coefficient_structure == 2) { /* brick */
876: PetscReal values_E[10];
877: PetscReal values_nu[10];
878: PetscInt nbricks, maxnbricks;
879: PetscInt index, span;
880: PetscInt jj;
882: flg = PETSC_FALSE;
883: maxnbricks = 10;
884: PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-brick_E", values_E, &maxnbricks, &flg));
885: nbricks = maxnbricks;
886: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply a list of E values for each brick");
888: flg = PETSC_FALSE;
889: maxnbricks = 10;
890: PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-brick_nu", values_nu, &maxnbricks, &flg));
891: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply a list of nu values for each brick");
892: PetscCheck(maxnbricks == nbricks, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply equal numbers of values for E and nu");
894: span = 1;
895: PetscCall(PetscOptionsGetInt(NULL, NULL, "-brick_span", &span, &flg));
897: /* cycle through the indices so that no two material properties are repeated in lines of x or y */
898: jj = (j / span) % nbricks;
899: index = (jj + i / span) % nbricks;
900: /*printf("j=%d: index = %d \n", j,index); */
902: for (p = 0; p < GAUSS_POINTS; p++) {
903: element_props[j][i].E[p] = values_E[index];
904: element_props[j][i].nu[p] = values_nu[index];
905: }
906: } else if (coefficient_structure == 3) { /* sponge */
907: PetscScalar opts_E0, opts_nu0;
908: PetscScalar opts_E1, opts_nu1;
909: PetscInt opts_t, opts_w;
910: PetscInt ii, jj, ci, cj;
912: opts_E0 = opts_E1 = 1.0;
913: opts_nu0 = opts_nu1 = 0.333;
914: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_E0", &opts_E0, &flg));
915: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_nu0", &opts_nu0, &flg));
916: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_E1", &opts_E1, &flg));
917: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_nu1", &opts_nu1, &flg));
919: opts_t = opts_w = 1;
920: PetscCall(PetscOptionsGetInt(NULL, NULL, "-sponge_t", &opts_t, &flg));
921: PetscCall(PetscOptionsGetInt(NULL, NULL, "-sponge_w", &opts_w, &flg));
923: ii = (i) / (opts_t + opts_w + opts_t);
924: jj = (j) / (opts_t + opts_w + opts_t);
926: ci = i - ii * (opts_t + opts_w + opts_t);
927: cj = j - jj * (opts_t + opts_w + opts_t);
929: for (p = 0; p < GAUSS_POINTS; p++) {
930: element_props[j][i].E[p] = opts_E0;
931: element_props[j][i].nu[p] = opts_nu0;
932: }
933: if ((ci >= opts_t) && (ci < opts_t + opts_w)) {
934: if ((cj >= opts_t) && (cj < opts_t + opts_w)) {
935: for (p = 0; p < GAUSS_POINTS; p++) {
936: element_props[j][i].E[p] = opts_E1;
937: element_props[j][i].nu[p] = opts_nu1;
938: }
939: }
940: }
941: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
942: }
943: }
944: PetscCall(DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords));
946: PetscCall(DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords));
948: PetscCall(DMDAVecRestoreArray(da_prop, l_properties, &element_props));
949: PetscCall(DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties));
950: PetscCall(DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties));
952: PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL));
953: if (!no_view) {
954: PetscCall(DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for elasticity eqn.", "properties"));
955: PetscCall(DMDACoordViewGnuplot2d(elas_da, "mesh"));
956: }
958: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
959: PetscCall(DMCreateMatrix(elas_da, &A));
960: PetscCall(DMGetCoordinates(elas_da, &vel_coords));
961: PetscCall(MatNullSpaceCreateRigidBody(vel_coords, &matnull));
962: PetscCall(MatSetNearNullSpace(A, matnull));
963: PetscCall(MatNullSpaceDestroy(&matnull));
964: PetscCall(MatCreateVecs(A, &f, &X));
966: /* assemble A11 */
967: PetscCall(MatZeroEntries(A));
969: PetscCall(AssembleA_Elasticity(A, elas_da, da_prop, properties));
970: /* build force vector */
971: PetscCall(AssembleF_Elasticity(f, elas_da, da_prop, properties));
973: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_E));
974: PetscCall(KSPSetOptionsPrefix(ksp_E, "elas_")); /* elasticity */
976: PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_nonsymbc", &use_nonsymbc, &flg));
977: /* solve */
978: if (!use_nonsymbc) {
979: Mat AA;
980: Vec ff, XX;
981: IS is;
982: VecScatter scat;
984: PetscCall(DMDABCApplySymmetricCompression(elas_da, A, f, &is, &AA, &ff));
985: PetscCall(VecDuplicate(ff, &XX));
987: PetscCall(KSPSetOperators(ksp_E, AA, AA));
988: PetscCall(KSPSetFromOptions(ksp_E));
990: PetscCall(KSPSolve(ksp_E, ff, XX));
992: /* push XX back into X */
993: PetscCall(DMDABCApplyCompression(elas_da, NULL, X));
995: PetscCall(VecScatterCreate(XX, NULL, X, is, &scat));
996: PetscCall(VecScatterBegin(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD));
997: PetscCall(VecScatterEnd(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD));
998: PetscCall(VecScatterDestroy(&scat));
1000: PetscCall(MatDestroy(&AA));
1001: PetscCall(VecDestroy(&ff));
1002: PetscCall(VecDestroy(&XX));
1003: PetscCall(ISDestroy(&is));
1004: } else {
1005: PetscCall(DMDABCApplyCompression(elas_da, A, f));
1007: PetscCall(KSPSetOperators(ksp_E, A, A));
1008: PetscCall(KSPSetFromOptions(ksp_E));
1010: PetscCall(KSPSolve(ksp_E, f, X));
1011: }
1013: if (!no_view) PetscCall(DMDAViewGnuplot2d(elas_da, X, "Displacement solution for elasticity eqn.", "X"));
1014: PetscCall(KSPDestroy(&ksp_E));
1016: PetscCall(VecDestroy(&X));
1017: PetscCall(VecDestroy(&f));
1018: PetscCall(MatDestroy(&A));
1020: PetscCall(DMDestroy(&elas_da));
1021: PetscCall(DMDestroy(&da_prop));
1023: PetscCall(VecDestroy(&properties));
1024: PetscCall(VecDestroy(&l_properties));
1025: PetscFunctionReturn(PETSC_SUCCESS);
1026: }
1028: int main(int argc, char **args)
1029: {
1030: PetscInt mx, my;
1032: PetscFunctionBeginUser;
1033: PetscCall(PetscInitialize(&argc, &args, NULL, help));
1034: mx = my = 10;
1035: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1036: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1037: PetscCall(solve_elasticity_2d(mx, my));
1038: PetscCall(PetscFinalize());
1039: return 0;
1040: }
1042: /* -------------------------- helpers for boundary conditions -------------------------------- */
1044: static PetscErrorCode BCApply_EAST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
1045: {
1046: DM cda;
1047: Vec coords;
1048: PetscInt si, sj, nx, ny, i, j;
1049: PetscInt M, N;
1050: DMDACoor2d **_coords;
1051: const PetscInt *g_idx;
1052: PetscInt *bc_global_ids;
1053: PetscScalar *bc_vals;
1054: PetscInt nbcs;
1055: PetscInt n_dofs;
1056: ISLocalToGlobalMapping ltogm;
1058: PetscFunctionBeginUser;
1059: /* enforce bc's */
1060: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1061: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1063: PetscCall(DMGetCoordinateDM(da, &cda));
1064: PetscCall(DMGetCoordinatesLocal(da, &coords));
1065: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1066: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1067: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1069: /* --- */
1071: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1072: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1074: /* init the entries to -1 so VecSetValues will ignore them */
1075: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1077: i = nx - 1;
1078: for (j = 0; j < ny; j++) {
1079: PetscInt local_id;
1080: PETSC_UNUSED PetscScalar coordx, coordy;
1082: local_id = i + j * nx;
1084: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1086: coordx = _coords[j + sj][i + si].x;
1087: coordy = _coords[j + sj][i + si].y;
1089: bc_vals[j] = bc_val;
1090: }
1091: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1092: nbcs = 0;
1093: if ((si + nx) == (M)) nbcs = ny;
1095: if (b) {
1096: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1097: PetscCall(VecAssemblyBegin(b));
1098: PetscCall(VecAssemblyEnd(b));
1099: }
1100: if (A) PetscCall(MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0));
1102: PetscCall(PetscFree(bc_vals));
1103: PetscCall(PetscFree(bc_global_ids));
1105: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1106: PetscFunctionReturn(PETSC_SUCCESS);
1107: }
1109: static PetscErrorCode BCApply_WEST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
1110: {
1111: DM cda;
1112: Vec coords;
1113: PetscInt si, sj, nx, ny, i, j;
1114: PetscInt M, N;
1115: DMDACoor2d **_coords;
1116: const PetscInt *g_idx;
1117: PetscInt *bc_global_ids;
1118: PetscScalar *bc_vals;
1119: PetscInt nbcs;
1120: PetscInt n_dofs;
1121: ISLocalToGlobalMapping ltogm;
1123: PetscFunctionBeginUser;
1124: /* enforce bc's */
1125: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1126: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1128: PetscCall(DMGetCoordinateDM(da, &cda));
1129: PetscCall(DMGetCoordinatesLocal(da, &coords));
1130: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1131: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1132: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1134: /* --- */
1136: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1137: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1139: /* init the entries to -1 so VecSetValues will ignore them */
1140: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1142: i = 0;
1143: for (j = 0; j < ny; j++) {
1144: PetscInt local_id;
1145: PETSC_UNUSED PetscScalar coordx, coordy;
1147: local_id = i + j * nx;
1149: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1151: coordx = _coords[j + sj][i + si].x;
1152: coordy = _coords[j + sj][i + si].y;
1154: bc_vals[j] = bc_val;
1155: }
1156: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1157: nbcs = 0;
1158: if (si == 0) nbcs = ny;
1160: if (b) {
1161: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1162: PetscCall(VecAssemblyBegin(b));
1163: PetscCall(VecAssemblyEnd(b));
1164: }
1165: if (A) PetscCall(MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0));
1167: PetscCall(PetscFree(bc_vals));
1168: PetscCall(PetscFree(bc_global_ids));
1170: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1171: PetscFunctionReturn(PETSC_SUCCESS);
1172: }
1174: static PetscErrorCode DMDABCApplyCompression(DM elas_da, Mat A, Vec f)
1175: {
1176: PetscFunctionBeginUser;
1177: PetscCall(BCApply_EAST(elas_da, 0, -1.0, A, f));
1178: PetscCall(BCApply_EAST(elas_da, 1, 0.0, A, f));
1179: PetscCall(BCApply_WEST(elas_da, 0, 1.0, A, f));
1180: PetscCall(BCApply_WEST(elas_da, 1, 0.0, A, f));
1181: PetscFunctionReturn(PETSC_SUCCESS);
1182: }
1184: static PetscErrorCode Orthogonalize(PetscInt n, Vec *vecs)
1185: {
1186: PetscScalar dot;
1188: PetscFunctionBegin;
1189: for (PetscInt i = 0; i < n; i++) {
1190: PetscCall(VecNormalize(vecs[i], NULL));
1191: for (PetscInt j = i + 1; j < n; j++) {
1192: PetscCall(VecDot(vecs[i], vecs[j], &dot));
1193: PetscCall(VecAXPY(vecs[j], -dot, vecs[i]));
1194: }
1195: }
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da, Mat A, Vec f, IS *dofs, Mat *AA, Vec *ff)
1200: {
1201: PetscInt start, end, m;
1202: PetscInt *unconstrained;
1203: PetscInt cnt;
1204: Vec x;
1205: PetscScalar *_x;
1206: IS is;
1207: VecScatter scat;
1209: PetscFunctionBeginUser;
1210: /* push bc's into f and A */
1211: PetscCall(VecDuplicate(f, &x));
1212: PetscCall(BCApply_EAST(elas_da, 0, -1.0, A, x));
1213: PetscCall(BCApply_EAST(elas_da, 1, 0.0, A, x));
1214: PetscCall(BCApply_WEST(elas_da, 0, 1.0, A, x));
1215: PetscCall(BCApply_WEST(elas_da, 1, 0.0, A, x));
1217: /* define which dofs are not constrained */
1218: PetscCall(VecGetLocalSize(x, &m));
1219: PetscCall(PetscMalloc1(m, &unconstrained));
1220: PetscCall(VecGetOwnershipRange(x, &start, &end));
1221: PetscCall(VecGetArray(x, &_x));
1222: cnt = 0;
1223: for (PetscInt i = 0; i < m; i += 2) {
1224: PetscReal val1, val2;
1226: val1 = PetscRealPart(_x[i]);
1227: val2 = PetscRealPart(_x[i + 1]);
1228: if (PetscAbs(val1) < 0.1 && PetscAbs(val2) < 0.1) {
1229: unconstrained[cnt] = start + i;
1230: cnt++;
1231: unconstrained[cnt] = start + i + 1;
1232: cnt++;
1233: }
1234: }
1235: PetscCall(VecRestoreArray(x, &_x));
1237: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, cnt, unconstrained, PETSC_COPY_VALUES, &is));
1238: PetscCall(PetscFree(unconstrained));
1239: PetscCall(ISSetBlockSize(is, 2));
1241: /* define correction for dirichlet in the rhs */
1242: PetscCall(MatMult(A, x, f));
1243: PetscCall(VecScale(f, -1.0));
1245: /* get new matrix */
1246: PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, AA));
1247: /* get new vector */
1248: PetscCall(MatCreateVecs(*AA, NULL, ff));
1250: PetscCall(VecScatterCreate(f, is, *ff, NULL, &scat));
1251: PetscCall(VecScatterBegin(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD));
1252: PetscCall(VecScatterEnd(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD));
1254: { /* Constrain near-null space */
1255: PetscInt nvecs;
1256: const Vec *vecs;
1257: Vec *uvecs;
1258: PetscBool has_const;
1259: MatNullSpace mnull, unull;
1261: PetscCall(MatGetNearNullSpace(A, &mnull));
1262: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvecs, &vecs));
1263: PetscCall(VecDuplicateVecs(*ff, nvecs, &uvecs));
1264: for (PetscInt i = 0; i < nvecs; i++) {
1265: PetscCall(VecScatterBegin(scat, vecs[i], uvecs[i], INSERT_VALUES, SCATTER_FORWARD));
1266: PetscCall(VecScatterEnd(scat, vecs[i], uvecs[i], INSERT_VALUES, SCATTER_FORWARD));
1267: }
1268: PetscCall(Orthogonalize(nvecs, uvecs));
1269: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)A), PETSC_FALSE, nvecs, uvecs, &unull));
1270: PetscCall(MatSetNearNullSpace(*AA, unull));
1271: PetscCall(MatNullSpaceDestroy(&unull));
1272: PetscCall(VecDestroyVecs(nvecs, &uvecs));
1273: }
1275: PetscCall(VecScatterDestroy(&scat));
1277: *dofs = is;
1278: PetscCall(VecDestroy(&x));
1279: PetscFunctionReturn(PETSC_SUCCESS);
1280: }
1282: /*TEST
1284: build:
1285: requires: !complex !single
1287: test:
1288: args: -mx 20 -my 30 -elas_ksp_monitor -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
1289: output_file: output/ex49_1.out
1291: test:
1292: suffix: 2
1293: nsize: 4
1294: args: -mx 20 -my 30 -elas_ksp_monitor -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
1296: test:
1297: suffix: 3
1298: nsize: 4
1299: args: -mx 20 -my 30 -elas_ksp_monitor -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
1301: test:
1302: suffix: 4
1303: nsize: 4
1304: args: -elas_ksp_monitor -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
1305: requires: ml
1307: test:
1308: suffix: 5
1309: nsize: 3
1310: args: -elas_ksp_monitor -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
1312: test:
1313: suffix: 6
1314: nsize: 4
1315: args: -mx 20 -my 30 -elas_ksp_monitor -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
1317: test:
1318: suffix: 7
1319: nsize: 4
1320: args: -mx 20 -my 30 -elas_ksp_monitor -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
1322: test:
1323: suffix: 8
1324: nsize: 4
1325: args: -mx 20 -my 30 -elas_ksp_monitor -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
1327: test:
1328: suffix: hypre_nullspace
1329: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
1330: args: -elas_ksp_monitor -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
1332: test:
1333: nsize: 4
1334: suffix: bddc
1335: args: -elas_ksp_monitor -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
1337: test:
1338: nsize: 4
1339: suffix: bddc_unsym
1340: args: -elas_ksp_monitor -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
1342: test:
1343: nsize: 4
1344: suffix: bddc_unsym_deluxe
1345: args: -elas_ksp_monitor -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
1347: test:
1348: nsize: 4
1349: suffix: fetidp_unsym_deluxe
1350: args: -elas_ksp_monitor -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
1352: test:
1353: nsize: 4
1354: suffix: bddc_layerjump
1355: args: -mx 40 -my 40 -elas_ksp_monitor -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
1357: test:
1358: nsize: 4
1359: suffix: bddc_subdomainjump
1360: args: -mx 40 -my 40 -elas_ksp_monitor -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
1362: test:
1363: nsize: 9
1364: suffix: bddc_subdomainjump_deluxe
1365: args: -mx 30 -my 30 -elas_ksp_monitor -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
1366: TEST*/