Actual source code: ex49.c
1: static char help[] = " Solves the compressible plane strain elasticity equations in 2d on the unit domain using Q1 finite elements. \n\
2: Material properties E (Youngs modulus) and nu (Poisson ratio) may vary as a function of space. \n\
3: The model utilises boundary conditions which produce compression in the x direction. \n\
4: Options: \n"
5: "\
6: -mx : number of elements in x-direction \n\
7: -my : number of elements in y-direction \n\
8: -c_str : structure of the coefficients to use. \n"
9: "\
10: -c_str 0 => isotropic material with constant coefficients. \n\
11: Parameters: \n\
12: -iso_E : Youngs modulus \n\
13: -iso_nu : Poisson ratio \n\
14: -c_str 1 => step function in the material properties in x. \n\
15: Parameters: \n\
16: -step_E0 : Youngs modulus to the left of the step \n\
17: -step_nu0 : Poisson ratio to the left of the step \n\
18: -step_E1 : Youngs modulus to the right of the step \n\
19: -step_n1 : Poisson ratio to the right of the step \n\
20: -step_xc : x coordinate of the step \n"
21: "\
22: -c_str 2 => checkerboard material with alternating properties. \n\
23: Repeats the following pattern throughout the domain. For example with 4 materials specified, we would heve \n\
24: -------------------------\n\
25: | D | A | B | C |\n\
26: ------|-----|-----|------\n\
27: | C | D | A | B |\n\
28: ------|-----|-----|------\n\
29: | B | C | D | A |\n\
30: ------|-----|-----|------\n\
31: | A | B | C | D |\n\
32: -------------------------\n\
33: \n\
34: Parameters: \n\
35: -brick_E : a comma separated list of Young's modulii \n\
36: -brick_nu : a comma separated list of Poisson ratios \n\
37: -brick_span : the number of elements in x and y each brick will span \n\
38: -c_str 3 => sponge-like material with alternating properties. \n\
39: Repeats the following pattern throughout the domain \n"
40: "\
41: -----------------------------\n\
42: | [background] |\n\
43: | E0,nu0 |\n\
44: | ----------------- |\n\
45: | | [inclusion] | |\n\
46: | | E1,nu1 | |\n\
47: | | | |\n\
48: | | <---- w ----> | |\n\
49: | | | |\n\
50: | | | |\n\
51: | ----------------- |\n\
52: | |\n\
53: | |\n\
54: -----------------------------\n\
55: <-------- t + w + t ------->\n\
56: \n\
57: Parameters: \n\
58: -sponge_E0 : Youngs modulus of the surrounding material \n\
59: -sponge_E1 : Youngs modulus of the inclusion \n\
60: -sponge_nu0 : Poisson ratio of the surrounding material \n\
61: -sponge_nu1 : Poisson ratio of the inclusion \n\
62: -sponge_t : the number of elements defining the border around each inclusion \n\
63: -sponge_w : the number of elements in x and y each inclusion will span\n\
64: -use_gp_coords : Evaluate the Youngs modulus, Poisson ratio and the body force at the global coordinates of the quadrature points.\n\
65: By default, E, nu and the body force are evaluated at the element center and applied as a constant over the entire element.\n\
66: -use_nonsymbc : Option to use non-symmetric boundary condition imposition. This choice will use less memory.";
68: /* Contributed by Dave May */
70: #include <petscksp.h>
71: #include <petscdm.h>
72: #include <petscdmda.h>
74: static PetscErrorCode DMDABCApplyCompression(DM, Mat, Vec);
75: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da, Mat A, Vec f, IS *dofs, Mat *AA, Vec *ff);
77: #define NSD 2 /* number of spatial dimensions */
78: #define NODES_PER_EL 4 /* nodes per element */
79: #define U_DOFS 2 /* degrees of freedom per displacement node */
80: #define GAUSS_POINTS 4
82: /* cell based evaluation */
83: typedef struct {
84: PetscScalar E, nu, fx, fy;
85: } Coefficients;
87: /* Gauss point based evaluation 8+4+4+4 = 20 */
88: typedef struct {
89: PetscScalar gp_coords[2 * GAUSS_POINTS];
90: PetscScalar E[GAUSS_POINTS];
91: PetscScalar nu[GAUSS_POINTS];
92: PetscScalar fx[GAUSS_POINTS];
93: PetscScalar fy[GAUSS_POINTS];
94: } GaussPointCoefficients;
96: typedef struct {
97: PetscScalar ux_dof;
98: PetscScalar uy_dof;
99: } ElasticityDOF;
101: /*
103: D = E/((1+nu)(1-2nu)) * [ 1-nu nu 0 ]
104: [ nu 1-nu 0 ]
105: [ 0 0 0.5*(1-2nu) ]
107: B = [ d_dx 0 ]
108: [ 0 d_dy ]
109: [ d_dy d_dx ]
111: */
113: /* FEM routines */
114: /*
115: Element: Local basis function ordering
116: 1-----2
117: | |
118: | |
119: 0-----3
120: */
121: static void ConstructQ12D_Ni(PetscScalar _xi[], PetscScalar Ni[])
122: {
123: PetscScalar xi = _xi[0];
124: PetscScalar eta = _xi[1];
126: Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
127: Ni[1] = 0.25 * (1.0 - xi) * (1.0 + eta);
128: Ni[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
129: Ni[3] = 0.25 * (1.0 + xi) * (1.0 - eta);
130: }
132: static void ConstructQ12D_GNi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
133: {
134: PetscScalar xi = _xi[0];
135: PetscScalar eta = _xi[1];
137: GNi[0][0] = -0.25 * (1.0 - eta);
138: GNi[0][1] = -0.25 * (1.0 + eta);
139: GNi[0][2] = 0.25 * (1.0 + eta);
140: GNi[0][3] = 0.25 * (1.0 - eta);
142: GNi[1][0] = -0.25 * (1.0 - xi);
143: GNi[1][1] = 0.25 * (1.0 - xi);
144: GNi[1][2] = 0.25 * (1.0 + xi);
145: GNi[1][3] = -0.25 * (1.0 + xi);
146: }
148: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL], PetscScalar GNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
149: {
150: PetscScalar J00, J01, J10, J11, J;
151: PetscScalar iJ00, iJ01, iJ10, iJ11;
152: PetscInt i;
154: J00 = J01 = J10 = J11 = 0.0;
155: for (i = 0; i < NODES_PER_EL; i++) {
156: PetscScalar cx = coords[2 * i + 0];
157: PetscScalar cy = coords[2 * i + 1];
159: J00 = J00 + GNi[0][i] * cx; /* J_xx = dx/dxi */
160: J01 = J01 + GNi[0][i] * cy; /* J_xy = dy/dxi */
161: J10 = J10 + GNi[1][i] * cx; /* J_yx = dx/deta */
162: J11 = J11 + GNi[1][i] * cy; /* J_yy = dy/deta */
163: }
164: J = (J00 * J11) - (J01 * J10);
166: iJ00 = J11 / J;
167: iJ01 = -J01 / J;
168: iJ10 = -J10 / J;
169: iJ11 = J00 / J;
171: for (i = 0; i < NODES_PER_EL; i++) {
172: GNx[0][i] = GNi[0][i] * iJ00 + GNi[1][i] * iJ01;
173: GNx[1][i] = GNi[0][i] * iJ10 + GNi[1][i] * iJ11;
174: }
176: if (det_J) *det_J = J;
177: }
179: static void ConstructGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
180: {
181: *ngp = 4;
182: gp_xi[0][0] = -0.57735026919;
183: gp_xi[0][1] = -0.57735026919;
184: gp_xi[1][0] = -0.57735026919;
185: gp_xi[1][1] = 0.57735026919;
186: gp_xi[2][0] = 0.57735026919;
187: gp_xi[2][1] = 0.57735026919;
188: gp_xi[3][0] = 0.57735026919;
189: gp_xi[3][1] = -0.57735026919;
190: gp_weight[0] = 1.0;
191: gp_weight[1] = 1.0;
192: gp_weight[2] = 1.0;
193: gp_weight[3] = 1.0;
194: }
196: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da, PetscInt **_lx, PetscInt **_ly)
197: {
198: PetscMPIInt rank;
199: PetscInt proc_I, proc_J;
200: PetscInt cpu_x, cpu_y;
201: PetscInt local_mx, local_my;
202: Vec vlx, vly;
203: PetscInt *LX, *LY, i;
204: PetscScalar *_a;
205: Vec V_SEQ;
206: VecScatter ctx;
208: PetscFunctionBeginUser;
209: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
211: PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
213: proc_J = rank / cpu_x;
214: proc_I = rank - cpu_x * proc_J;
216: PetscCall(PetscMalloc1(cpu_x, &LX));
217: PetscCall(PetscMalloc1(cpu_y, &LY));
219: PetscCall(DMDAGetElementsSizes(da, &local_mx, &local_my, NULL));
220: PetscCall(VecCreate(PETSC_COMM_WORLD, &vlx));
221: PetscCall(VecSetSizes(vlx, PETSC_DECIDE, cpu_x));
222: PetscCall(VecSetFromOptions(vlx));
224: PetscCall(VecCreate(PETSC_COMM_WORLD, &vly));
225: PetscCall(VecSetSizes(vly, PETSC_DECIDE, cpu_y));
226: PetscCall(VecSetFromOptions(vly));
228: PetscCall(VecSetValue(vlx, proc_I, (PetscScalar)(local_mx + 1.0e-9), INSERT_VALUES));
229: PetscCall(VecSetValue(vly, proc_J, (PetscScalar)(local_my + 1.0e-9), INSERT_VALUES));
230: PetscCall(VecAssemblyBegin(vlx));
231: PetscCall(VecAssemblyEnd(vlx));
232: PetscCall(VecAssemblyBegin(vly));
233: PetscCall(VecAssemblyEnd(vly));
235: PetscCall(VecScatterCreateToAll(vlx, &ctx, &V_SEQ));
236: PetscCall(VecScatterBegin(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
237: PetscCall(VecScatterEnd(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
238: PetscCall(VecGetArray(V_SEQ, &_a));
239: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
240: PetscCall(VecRestoreArray(V_SEQ, &_a));
241: PetscCall(VecScatterDestroy(&ctx));
242: PetscCall(VecDestroy(&V_SEQ));
244: PetscCall(VecScatterCreateToAll(vly, &ctx, &V_SEQ));
245: PetscCall(VecScatterBegin(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
246: PetscCall(VecScatterEnd(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
247: PetscCall(VecGetArray(V_SEQ, &_a));
248: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
249: PetscCall(VecRestoreArray(V_SEQ, &_a));
250: PetscCall(VecScatterDestroy(&ctx));
251: PetscCall(VecDestroy(&V_SEQ));
253: *_lx = LX;
254: *_ly = LY;
256: PetscCall(VecDestroy(&vlx));
257: PetscCall(VecDestroy(&vly));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: static PetscErrorCode DMDACoordViewGnuplot2d(DM da, const char prefix[])
262: {
263: DM cda;
264: Vec coords;
265: DMDACoor2d **_coords;
266: PetscInt si, sj, nx, ny, i, j;
267: FILE *fp;
268: char fname[PETSC_MAX_PATH_LEN];
269: PetscMPIInt rank;
271: PetscFunctionBeginUser;
272: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
273: PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
274: PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
275: PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
276: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### Element geometry for processor %1.4d ### \n", rank));
278: PetscCall(DMGetCoordinateDM(da, &cda));
279: PetscCall(DMGetCoordinatesLocal(da, &coords));
280: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
281: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
282: for (j = sj; j < sj + ny - 1; j++) {
283: for (i = si; i < si + nx - 1; i++) {
284: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
285: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i].x), (double)PetscRealPart(_coords[j + 1][i].y)));
286: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i + 1].x), (double)PetscRealPart(_coords[j + 1][i + 1].y)));
287: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i + 1].x), (double)PetscRealPart(_coords[j][i + 1].y)));
288: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n\n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
289: }
290: }
291: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
293: PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: static PetscErrorCode DMDAViewGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
298: {
299: DM cda;
300: Vec coords, local_fields;
301: DMDACoor2d **_coords;
302: FILE *fp;
303: char fname[PETSC_MAX_PATH_LEN];
304: const char *field_name;
305: PetscMPIInt rank;
306: PetscInt si, sj, nx, ny, i, j;
307: PetscInt n_dofs, d;
308: PetscScalar *_fields;
310: PetscFunctionBeginUser;
311: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
312: PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
313: PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
314: PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
316: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
317: PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
318: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
319: for (d = 0; d < n_dofs; d++) {
320: PetscCall(DMDAGetFieldName(da, d, &field_name));
321: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
322: }
323: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));
325: PetscCall(DMGetCoordinateDM(da, &cda));
326: PetscCall(DMGetCoordinatesLocal(da, &coords));
327: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
328: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
330: PetscCall(DMCreateLocalVector(da, &local_fields));
331: PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
332: PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
333: PetscCall(VecGetArray(local_fields, &_fields));
335: for (j = sj; j < sj + ny; j++) {
336: for (i = si; i < si + nx; i++) {
337: PetscScalar coord_x, coord_y;
338: PetscScalar field_d;
340: coord_x = _coords[j][i].x;
341: coord_y = _coords[j][i].y;
343: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
344: for (d = 0; d < n_dofs; d++) {
345: field_d = _fields[n_dofs * ((i - si) + (j - sj) * (nx)) + d];
346: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e ", (double)PetscRealPart(field_d)));
347: }
348: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
349: }
350: }
351: PetscCall(VecRestoreArray(local_fields, &_fields));
352: PetscCall(VecDestroy(&local_fields));
354: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
356: PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
361: {
362: DM cda;
363: Vec local_fields;
364: FILE *fp;
365: char fname[PETSC_MAX_PATH_LEN];
366: const char *field_name;
367: PetscMPIInt rank;
368: PetscInt si, sj, nx, ny, i, j, p;
369: PetscInt n_dofs, d;
370: GaussPointCoefficients **_coefficients;
372: PetscFunctionBeginUser;
373: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
374: PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
375: PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
376: PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
378: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
379: PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
380: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
381: for (d = 0; d < n_dofs; d++) {
382: PetscCall(DMDAGetFieldName(da, d, &field_name));
383: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
384: }
385: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));
387: PetscCall(DMGetCoordinateDM(da, &cda));
388: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
390: PetscCall(DMCreateLocalVector(da, &local_fields));
391: PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
392: PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
393: PetscCall(DMDAVecGetArray(da, local_fields, &_coefficients));
395: for (j = sj; j < sj + ny; j++) {
396: for (i = si; i < si + nx; i++) {
397: PetscScalar coord_x, coord_y;
399: for (p = 0; p < GAUSS_POINTS; p++) {
400: coord_x = _coefficients[j][i].gp_coords[2 * p];
401: coord_y = _coefficients[j][i].gp_coords[2 * p + 1];
403: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
405: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e %1.6e %1.6e\n", (double)PetscRealPart(_coefficients[j][i].E[p]), (double)PetscRealPart(_coefficients[j][i].nu[p]), (double)PetscRealPart(_coefficients[j][i].fx[p]),
406: (double)PetscRealPart(_coefficients[j][i].fy[p])));
407: }
408: }
409: }
410: PetscCall(DMDAVecRestoreArray(da, local_fields, &_coefficients));
411: PetscCall(VecDestroy(&local_fields));
413: PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: static void FormStressOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar E[], PetscScalar nu[])
418: {
419: PetscInt ngp;
420: PetscScalar gp_xi[GAUSS_POINTS][2];
421: PetscScalar gp_weight[GAUSS_POINTS];
422: PetscInt p, i, j, k, l;
423: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
424: PetscScalar J_p;
425: PetscScalar B[3][U_DOFS * NODES_PER_EL];
426: PetscScalar prop_E, prop_nu, factor, constit_D[3][3];
428: /* define quadrature rule */
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 ei, ej;
574: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
575: PetscScalar el_coords[NODES_PER_EL * NSD];
576: Vec local_properties;
577: GaussPointCoefficients **props;
578: PetscScalar *prop_E, *prop_nu;
580: PetscFunctionBeginUser;
581: /* setup for coords */
582: PetscCall(DMGetCoordinateDM(elas_da, &cda));
583: PetscCall(DMGetCoordinatesLocal(elas_da, &coords));
584: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
586: /* setup for coefficients */
587: PetscCall(DMCreateLocalVector(properties_da, &local_properties));
588: PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
589: PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
590: PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));
592: PetscCall(DMDAGetElementsCorners(elas_da, &sex, &sey, 0));
593: PetscCall(DMDAGetElementsSizes(elas_da, &mx, &my, 0));
594: for (ej = sey; ej < sey + my; ej++) {
595: for (ei = sex; ei < sex + mx; ei++) {
596: /* get coords for the element */
597: PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
599: /* get coefficients for the element */
600: prop_E = props[ej][ei].E;
601: prop_nu = props[ej][ei].nu;
603: /* initialise element stiffness matrix */
604: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
606: /* form element stiffness matrix */
607: FormStressOperatorQ1(Ae, el_coords, prop_E, prop_nu);
609: /* insert element matrix into global matrix */
610: PetscCall(DMDAGetElementEqnums_u(u_eqn, ei, ej));
611: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
612: }
613: }
614: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
615: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
617: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
619: PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
620: PetscCall(VecDestroy(&local_properties));
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F, MatStencil u_eqn[], PetscScalar Fe_u[])
625: {
626: PetscInt n;
628: PetscFunctionBeginUser;
629: for (n = 0; n < 4; n++) {
630: fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].ux_dof = fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].ux_dof + Fe_u[2 * n];
631: fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].uy_dof = fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].uy_dof + Fe_u[2 * n + 1];
632: }
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: static PetscErrorCode AssembleF_Elasticity(Vec F, DM elas_da, DM properties_da, Vec properties)
637: {
638: DM cda;
639: Vec coords;
640: DMDACoor2d **_coords;
641: MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
642: PetscInt sex, sey, mx, my;
643: PetscInt ei, ej;
644: PetscScalar Fe[NODES_PER_EL * U_DOFS];
645: PetscScalar el_coords[NODES_PER_EL * NSD];
646: Vec local_properties;
647: GaussPointCoefficients **props;
648: PetscScalar *prop_fx, *prop_fy;
649: Vec local_F;
650: ElasticityDOF **ff;
652: PetscFunctionBeginUser;
653: /* setup for coords */
654: PetscCall(DMGetCoordinateDM(elas_da, &cda));
655: PetscCall(DMGetCoordinatesLocal(elas_da, &coords));
656: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
658: /* setup for coefficients */
659: PetscCall(DMGetLocalVector(properties_da, &local_properties));
660: PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
661: PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
662: PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));
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));
669: PetscCall(DMDAGetElementsCorners(elas_da, &sex, &sey, 0));
670: PetscCall(DMDAGetElementsSizes(elas_da, &mx, &my, 0));
671: for (ej = sey; ej < sey + my; ej++) {
672: for (ei = sex; ei < sex + mx; ei++) {
673: /* get coords for the element */
674: PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
676: /* get coefficients for the element */
677: prop_fx = props[ej][ei].fx;
678: prop_fy = props[ej][ei].fy;
680: /* initialise element stiffness matrix */
681: PetscCall(PetscMemzero(Fe, sizeof(Fe)));
683: /* form element stiffness matrix */
684: FormMomentumRhsQ1(Fe, el_coords, prop_fx, prop_fy);
686: /* insert element matrix into global matrix */
687: PetscCall(DMDAGetElementEqnums_u(u_eqn, ei, ej));
689: PetscCall(DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, Fe));
690: }
691: }
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));
698: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
700: PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
701: PetscCall(DMRestoreLocalVector(properties_da, &local_properties));
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
705: static PetscErrorCode solve_elasticity_2d(PetscInt mx, PetscInt my)
706: {
707: DM elas_da, da_prop;
708: PetscInt u_dof, dof, stencil_width;
709: Mat A;
710: PetscInt mxl, myl;
711: DM prop_cda, vel_cda;
712: Vec prop_coords, vel_coords;
713: PetscInt si, sj, nx, ny, i, j, p;
714: Vec f, X;
715: PetscInt prop_dof, prop_stencil_width;
716: Vec properties, l_properties;
717: MatNullSpace matnull;
718: PetscReal dx, dy;
719: PetscInt M, N;
720: DMDACoor2d **_prop_coords, **_vel_coords;
721: GaussPointCoefficients **element_props;
722: KSP ksp_E;
723: PetscInt coefficient_structure = 0;
724: PetscInt cpu_x, cpu_y, *lx = NULL, *ly = NULL;
725: PetscBool use_gp_coords = PETSC_FALSE;
726: PetscBool use_nonsymbc = PETSC_FALSE;
727: PetscBool no_view = PETSC_FALSE;
728: PetscBool flg;
730: PetscFunctionBeginUser;
731: /* Generate the da for velocity and pressure */
732: /*
733: We use Q1 elements for the temperature.
734: FEM has a 9-point stencil (BOX) or connectivity pattern
735: Num nodes in each direction is mx+1, my+1
736: */
737: u_dof = U_DOFS; /* Vx, Vy - velocities */
738: dof = u_dof;
739: stencil_width = 1;
740: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, &elas_da));
742: PetscCall(DMSetMatType(elas_da, MATAIJ));
743: PetscCall(DMSetFromOptions(elas_da));
744: PetscCall(DMSetUp(elas_da));
746: PetscCall(DMDASetFieldName(elas_da, 0, "Ux"));
747: PetscCall(DMDASetFieldName(elas_da, 1, "Uy"));
749: /* unit box [0,1] x [0,1] */
750: PetscCall(DMDASetUniformCoordinates(elas_da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
752: /* Generate element properties, we will assume all material properties are constant over the element */
753: /* local number of elements */
754: PetscCall(DMDAGetElementsSizes(elas_da, &mxl, &myl, NULL));
756: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
757: PetscCall(DMDAGetInfo(elas_da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
758: PetscCall(DMDAGetElementOwnershipRanges2d(elas_da, &lx, &ly));
760: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
761: prop_stencil_width = 0;
762: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, cpu_x, cpu_y, prop_dof, prop_stencil_width, lx, ly, &da_prop));
763: PetscCall(DMSetFromOptions(da_prop));
764: PetscCall(DMSetUp(da_prop));
766: PetscCall(PetscFree(lx));
767: PetscCall(PetscFree(ly));
769: /* define centroid positions */
770: PetscCall(DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
771: dx = 1.0 / ((PetscReal)(M));
772: dy = 1.0 / ((PetscReal)(N));
774: PetscCall(DMDASetUniformCoordinates(da_prop, 0.0 + 0.5 * dx, 1.0 - 0.5 * dx, 0.0 + 0.5 * dy, 1.0 - 0.5 * dy, 0.0, 1.0));
776: /* define coefficients */
777: PetscCall(PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL));
779: PetscCall(DMCreateGlobalVector(da_prop, &properties));
780: PetscCall(DMCreateLocalVector(da_prop, &l_properties));
781: PetscCall(DMDAVecGetArray(da_prop, l_properties, &element_props));
783: PetscCall(DMGetCoordinateDM(da_prop, &prop_cda));
784: PetscCall(DMGetCoordinatesLocal(da_prop, &prop_coords));
785: PetscCall(DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords));
787: PetscCall(DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0));
789: PetscCall(DMGetCoordinateDM(elas_da, &vel_cda));
790: PetscCall(DMGetCoordinatesLocal(elas_da, &vel_coords));
791: PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));
793: /* interpolate the coordinates */
794: for (j = sj; j < sj + ny; j++) {
795: for (i = si; i < si + nx; i++) {
796: PetscInt ngp;
797: PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
798: PetscScalar el_coords[8];
800: PetscCall(GetElementCoords(_vel_coords, i, j, el_coords));
801: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
803: for (p = 0; p < GAUSS_POINTS; p++) {
804: PetscScalar gp_x, gp_y;
805: PetscInt n;
806: PetscScalar xi_p[2], Ni_p[4];
808: xi_p[0] = gp_xi[p][0];
809: xi_p[1] = gp_xi[p][1];
810: ConstructQ12D_Ni(xi_p, Ni_p);
812: gp_x = 0.0;
813: gp_y = 0.0;
814: for (n = 0; n < NODES_PER_EL; n++) {
815: gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
816: gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
817: }
818: element_props[j][i].gp_coords[2 * p] = gp_x;
819: element_props[j][i].gp_coords[2 * p + 1] = gp_y;
820: }
821: }
822: }
824: /* define the coefficients */
825: PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, &flg));
827: for (j = sj; j < sj + ny; j++) {
828: for (i = si; i < si + nx; i++) {
829: PetscScalar centroid_x = _prop_coords[j][i].x; /* centroids of cell */
830: PetscScalar centroid_y = _prop_coords[j][i].y;
831: PETSC_UNUSED PetscScalar coord_x, coord_y;
833: if (coefficient_structure == 0) { /* isotropic */
834: PetscScalar opts_E, opts_nu;
836: opts_E = 1.0;
837: opts_nu = 0.33;
838: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-iso_E", &opts_E, &flg));
839: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-iso_nu", &opts_nu, &flg));
841: for (p = 0; p < GAUSS_POINTS; p++) {
842: element_props[j][i].E[p] = opts_E;
843: element_props[j][i].nu[p] = opts_nu;
845: element_props[j][i].fx[p] = 0.0;
846: element_props[j][i].fy[p] = 0.0;
847: }
848: } else if (coefficient_structure == 1) { /* step */
849: PetscScalar opts_E0, opts_nu0, opts_xc;
850: PetscScalar opts_E1, opts_nu1;
852: opts_E0 = opts_E1 = 1.0;
853: opts_nu0 = opts_nu1 = 0.333;
854: opts_xc = 0.5;
855: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_E0", &opts_E0, &flg));
856: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_nu0", &opts_nu0, &flg));
857: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_E1", &opts_E1, &flg));
858: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_nu1", &opts_nu1, &flg));
859: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_xc", &opts_xc, &flg));
861: for (p = 0; p < GAUSS_POINTS; p++) {
862: coord_x = centroid_x;
863: coord_y = centroid_y;
864: if (use_gp_coords) {
865: coord_x = element_props[j][i].gp_coords[2 * p];
866: coord_y = element_props[j][i].gp_coords[2 * p + 1];
867: }
869: element_props[j][i].E[p] = opts_E0;
870: element_props[j][i].nu[p] = opts_nu0;
871: if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
872: element_props[j][i].E[p] = opts_E1;
873: element_props[j][i].nu[p] = opts_nu1;
874: }
876: element_props[j][i].fx[p] = 0.0;
877: element_props[j][i].fy[p] = 0.0;
878: }
879: } else if (coefficient_structure == 2) { /* brick */
880: PetscReal values_E[10];
881: PetscReal values_nu[10];
882: PetscInt nbricks, maxnbricks;
883: PetscInt index, span;
884: PetscInt jj;
886: flg = PETSC_FALSE;
887: maxnbricks = 10;
888: PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-brick_E", values_E, &maxnbricks, &flg));
889: nbricks = maxnbricks;
890: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply a list of E values for each brick");
892: flg = PETSC_FALSE;
893: maxnbricks = 10;
894: PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-brick_nu", values_nu, &maxnbricks, &flg));
895: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply a list of nu values for each brick");
896: PetscCheck(maxnbricks == nbricks, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply equal numbers of values for E and nu");
898: span = 1;
899: PetscCall(PetscOptionsGetInt(NULL, NULL, "-brick_span", &span, &flg));
901: /* cycle through the indices so that no two material properties are repeated in lines of x or y */
902: jj = (j / span) % nbricks;
903: index = (jj + i / span) % nbricks;
904: /*printf("j=%d: index = %d \n", j,index); */
906: for (p = 0; p < GAUSS_POINTS; p++) {
907: element_props[j][i].E[p] = values_E[index];
908: element_props[j][i].nu[p] = values_nu[index];
909: }
910: } else if (coefficient_structure == 3) { /* sponge */
911: PetscScalar opts_E0, opts_nu0;
912: PetscScalar opts_E1, opts_nu1;
913: PetscInt opts_t, opts_w;
914: PetscInt ii, jj, ci, cj;
916: opts_E0 = opts_E1 = 1.0;
917: opts_nu0 = opts_nu1 = 0.333;
918: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_E0", &opts_E0, &flg));
919: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_nu0", &opts_nu0, &flg));
920: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_E1", &opts_E1, &flg));
921: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_nu1", &opts_nu1, &flg));
923: opts_t = opts_w = 1;
924: PetscCall(PetscOptionsGetInt(NULL, NULL, "-sponge_t", &opts_t, &flg));
925: PetscCall(PetscOptionsGetInt(NULL, NULL, "-sponge_w", &opts_w, &flg));
927: ii = (i) / (opts_t + opts_w + opts_t);
928: jj = (j) / (opts_t + opts_w + opts_t);
930: ci = i - ii * (opts_t + opts_w + opts_t);
931: cj = j - jj * (opts_t + opts_w + opts_t);
933: for (p = 0; p < GAUSS_POINTS; p++) {
934: element_props[j][i].E[p] = opts_E0;
935: element_props[j][i].nu[p] = opts_nu0;
936: }
937: if ((ci >= opts_t) && (ci < opts_t + opts_w)) {
938: if ((cj >= opts_t) && (cj < opts_t + opts_w)) {
939: for (p = 0; p < GAUSS_POINTS; p++) {
940: element_props[j][i].E[p] = opts_E1;
941: element_props[j][i].nu[p] = opts_nu1;
942: }
943: }
944: }
945: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
946: }
947: }
948: PetscCall(DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords));
950: PetscCall(DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords));
952: PetscCall(DMDAVecRestoreArray(da_prop, l_properties, &element_props));
953: PetscCall(DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties));
954: PetscCall(DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties));
956: PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL));
957: if (!no_view) {
958: PetscCall(DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for elasticity eqn.", "properties"));
959: PetscCall(DMDACoordViewGnuplot2d(elas_da, "mesh"));
960: }
962: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
963: PetscCall(DMCreateMatrix(elas_da, &A));
964: PetscCall(DMGetCoordinates(elas_da, &vel_coords));
965: PetscCall(MatNullSpaceCreateRigidBody(vel_coords, &matnull));
966: PetscCall(MatSetNearNullSpace(A, matnull));
967: PetscCall(MatNullSpaceDestroy(&matnull));
968: PetscCall(MatCreateVecs(A, &f, &X));
970: /* assemble A11 */
971: PetscCall(MatZeroEntries(A));
972: PetscCall(VecZeroEntries(f));
974: PetscCall(AssembleA_Elasticity(A, elas_da, da_prop, properties));
975: /* build force vector */
976: PetscCall(AssembleF_Elasticity(f, elas_da, da_prop, properties));
978: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_E));
979: PetscCall(KSPSetOptionsPrefix(ksp_E, "elas_")); /* elasticity */
981: PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_nonsymbc", &use_nonsymbc, &flg));
982: /* solve */
983: if (!use_nonsymbc) {
984: Mat AA;
985: Vec ff, XX;
986: IS is;
987: VecScatter scat;
989: PetscCall(DMDABCApplySymmetricCompression(elas_da, A, f, &is, &AA, &ff));
990: PetscCall(VecDuplicate(ff, &XX));
992: PetscCall(KSPSetOperators(ksp_E, AA, AA));
993: PetscCall(KSPSetFromOptions(ksp_E));
995: PetscCall(KSPSolve(ksp_E, ff, XX));
997: /* push XX back into X */
998: PetscCall(DMDABCApplyCompression(elas_da, NULL, X));
1000: PetscCall(VecScatterCreate(XX, NULL, X, is, &scat));
1001: PetscCall(VecScatterBegin(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD));
1002: PetscCall(VecScatterEnd(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD));
1003: PetscCall(VecScatterDestroy(&scat));
1005: PetscCall(MatDestroy(&AA));
1006: PetscCall(VecDestroy(&ff));
1007: PetscCall(VecDestroy(&XX));
1008: PetscCall(ISDestroy(&is));
1009: } else {
1010: PetscCall(DMDABCApplyCompression(elas_da, A, f));
1012: PetscCall(KSPSetOperators(ksp_E, A, A));
1013: PetscCall(KSPSetFromOptions(ksp_E));
1015: PetscCall(KSPSolve(ksp_E, f, X));
1016: }
1018: if (!no_view) PetscCall(DMDAViewGnuplot2d(elas_da, X, "Displacement solution for elasticity eqn.", "X"));
1019: PetscCall(KSPDestroy(&ksp_E));
1021: PetscCall(VecDestroy(&X));
1022: PetscCall(VecDestroy(&f));
1023: PetscCall(MatDestroy(&A));
1025: PetscCall(DMDestroy(&elas_da));
1026: PetscCall(DMDestroy(&da_prop));
1028: PetscCall(VecDestroy(&properties));
1029: PetscCall(VecDestroy(&l_properties));
1030: PetscFunctionReturn(PETSC_SUCCESS);
1031: }
1033: int main(int argc, char **args)
1034: {
1035: PetscInt mx, my;
1037: PetscFunctionBeginUser;
1038: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
1039: mx = my = 10;
1040: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1041: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1042: PetscCall(solve_elasticity_2d(mx, my));
1043: PetscCall(PetscFinalize());
1044: return 0;
1045: }
1047: /* -------------------------- helpers for boundary conditions -------------------------------- */
1049: static PetscErrorCode BCApply_EAST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
1050: {
1051: DM cda;
1052: Vec coords;
1053: PetscInt si, sj, nx, ny, i, j;
1054: PetscInt M, N;
1055: DMDACoor2d **_coords;
1056: const PetscInt *g_idx;
1057: PetscInt *bc_global_ids;
1058: PetscScalar *bc_vals;
1059: PetscInt nbcs;
1060: PetscInt n_dofs;
1061: ISLocalToGlobalMapping ltogm;
1063: PetscFunctionBeginUser;
1064: /* enforce bc's */
1065: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1066: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1068: PetscCall(DMGetCoordinateDM(da, &cda));
1069: PetscCall(DMGetCoordinatesLocal(da, &coords));
1070: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1071: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1072: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1074: /* --- */
1076: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1077: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1079: /* init the entries to -1 so VecSetValues will ignore them */
1080: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1082: i = nx - 1;
1083: for (j = 0; j < ny; j++) {
1084: PetscInt local_id;
1085: PETSC_UNUSED PetscScalar coordx, coordy;
1087: local_id = i + j * nx;
1089: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1091: coordx = _coords[j + sj][i + si].x;
1092: coordy = _coords[j + sj][i + si].y;
1094: bc_vals[j] = bc_val;
1095: }
1096: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1097: nbcs = 0;
1098: if ((si + nx) == (M)) nbcs = ny;
1100: if (b) {
1101: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1102: PetscCall(VecAssemblyBegin(b));
1103: PetscCall(VecAssemblyEnd(b));
1104: }
1105: if (A) PetscCall(MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0));
1107: PetscCall(PetscFree(bc_vals));
1108: PetscCall(PetscFree(bc_global_ids));
1110: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1111: PetscFunctionReturn(PETSC_SUCCESS);
1112: }
1114: static PetscErrorCode BCApply_WEST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
1115: {
1116: DM cda;
1117: Vec coords;
1118: PetscInt si, sj, nx, ny, i, j;
1119: PetscInt M, N;
1120: DMDACoor2d **_coords;
1121: const PetscInt *g_idx;
1122: PetscInt *bc_global_ids;
1123: PetscScalar *bc_vals;
1124: PetscInt nbcs;
1125: PetscInt n_dofs;
1126: ISLocalToGlobalMapping ltogm;
1128: PetscFunctionBeginUser;
1129: /* enforce bc's */
1130: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1131: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1133: PetscCall(DMGetCoordinateDM(da, &cda));
1134: PetscCall(DMGetCoordinatesLocal(da, &coords));
1135: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1136: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1137: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1139: /* --- */
1141: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1142: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1144: /* init the entries to -1 so VecSetValues will ignore them */
1145: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1147: i = 0;
1148: for (j = 0; j < ny; j++) {
1149: PetscInt local_id;
1150: PETSC_UNUSED PetscScalar coordx, coordy;
1152: local_id = i + j * nx;
1154: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1156: coordx = _coords[j + sj][i + si].x;
1157: coordy = _coords[j + sj][i + si].y;
1159: bc_vals[j] = bc_val;
1160: }
1161: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1162: nbcs = 0;
1163: if (si == 0) nbcs = ny;
1165: if (b) {
1166: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1167: PetscCall(VecAssemblyBegin(b));
1168: PetscCall(VecAssemblyEnd(b));
1169: }
1170: if (A) PetscCall(MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0));
1172: PetscCall(PetscFree(bc_vals));
1173: PetscCall(PetscFree(bc_global_ids));
1175: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1176: PetscFunctionReturn(PETSC_SUCCESS);
1177: }
1179: static PetscErrorCode DMDABCApplyCompression(DM elas_da, Mat A, Vec f)
1180: {
1181: PetscFunctionBeginUser;
1182: PetscCall(BCApply_EAST(elas_da, 0, -1.0, A, f));
1183: PetscCall(BCApply_EAST(elas_da, 1, 0.0, A, f));
1184: PetscCall(BCApply_WEST(elas_da, 0, 1.0, A, f));
1185: PetscCall(BCApply_WEST(elas_da, 1, 0.0, A, f));
1186: PetscFunctionReturn(PETSC_SUCCESS);
1187: }
1189: static PetscErrorCode Orthogonalize(PetscInt n, Vec *vecs)
1190: {
1191: PetscInt i, j;
1192: PetscScalar dot;
1194: PetscFunctionBegin;
1195: for (i = 0; i < n; i++) {
1196: PetscCall(VecNormalize(vecs[i], NULL));
1197: for (j = i + 1; j < n; j++) {
1198: PetscCall(VecDot(vecs[i], vecs[j], &dot));
1199: PetscCall(VecAXPY(vecs[j], -dot, vecs[i]));
1200: }
1201: }
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da, Mat A, Vec f, IS *dofs, Mat *AA, Vec *ff)
1206: {
1207: PetscInt start, end, m;
1208: PetscInt *unconstrained;
1209: PetscInt cnt, i;
1210: Vec x;
1211: PetscScalar *_x;
1212: IS is;
1213: VecScatter scat;
1215: PetscFunctionBeginUser;
1216: /* push bc's into f and A */
1217: PetscCall(VecDuplicate(f, &x));
1218: PetscCall(BCApply_EAST(elas_da, 0, -1.0, A, x));
1219: PetscCall(BCApply_EAST(elas_da, 1, 0.0, A, x));
1220: PetscCall(BCApply_WEST(elas_da, 0, 1.0, A, x));
1221: PetscCall(BCApply_WEST(elas_da, 1, 0.0, A, x));
1223: /* define which dofs are not constrained */
1224: PetscCall(VecGetLocalSize(x, &m));
1225: PetscCall(PetscMalloc1(m, &unconstrained));
1226: PetscCall(VecGetOwnershipRange(x, &start, &end));
1227: PetscCall(VecGetArray(x, &_x));
1228: cnt = 0;
1229: for (i = 0; i < m; i += 2) {
1230: PetscReal val1, val2;
1232: val1 = PetscRealPart(_x[i]);
1233: val2 = PetscRealPart(_x[i + 1]);
1234: if (PetscAbs(val1) < 0.1 && PetscAbs(val2) < 0.1) {
1235: unconstrained[cnt] = start + i;
1236: cnt++;
1237: unconstrained[cnt] = start + i + 1;
1238: cnt++;
1239: }
1240: }
1241: PetscCall(VecRestoreArray(x, &_x));
1243: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, cnt, unconstrained, PETSC_COPY_VALUES, &is));
1244: PetscCall(PetscFree(unconstrained));
1245: PetscCall(ISSetBlockSize(is, 2));
1247: /* define correction for dirichlet in the rhs */
1248: PetscCall(MatMult(A, x, f));
1249: PetscCall(VecScale(f, -1.0));
1251: /* get new matrix */
1252: PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, AA));
1253: /* get new vector */
1254: PetscCall(MatCreateVecs(*AA, NULL, ff));
1256: PetscCall(VecScatterCreate(f, is, *ff, NULL, &scat));
1257: PetscCall(VecScatterBegin(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD));
1258: PetscCall(VecScatterEnd(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD));
1260: { /* Constrain near-null space */
1261: PetscInt nvecs;
1262: const Vec *vecs;
1263: Vec *uvecs;
1264: PetscBool has_const;
1265: MatNullSpace mnull, unull;
1267: PetscCall(MatGetNearNullSpace(A, &mnull));
1268: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvecs, &vecs));
1269: PetscCall(VecDuplicateVecs(*ff, nvecs, &uvecs));
1270: for (i = 0; i < nvecs; i++) {
1271: PetscCall(VecScatterBegin(scat, vecs[i], uvecs[i], INSERT_VALUES, SCATTER_FORWARD));
1272: PetscCall(VecScatterEnd(scat, vecs[i], uvecs[i], INSERT_VALUES, SCATTER_FORWARD));
1273: }
1274: PetscCall(Orthogonalize(nvecs, uvecs));
1275: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)A), PETSC_FALSE, nvecs, uvecs, &unull));
1276: PetscCall(MatSetNearNullSpace(*AA, unull));
1277: PetscCall(MatNullSpaceDestroy(&unull));
1278: PetscCall(VecDestroyVecs(nvecs, &uvecs));
1279: }
1281: PetscCall(VecScatterDestroy(&scat));
1283: *dofs = is;
1284: PetscCall(VecDestroy(&x));
1285: PetscFunctionReturn(PETSC_SUCCESS);
1286: }
1288: /*TEST
1290: build:
1291: requires: !complex !single
1293: test:
1294: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_rtol 5e-3 -elas_ksp_view
1295: output_file: output/ex49_1.out
1297: test:
1298: suffix: 2
1299: nsize: 4
1300: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type gcr -elas_pc_type asm -elas_sub_pc_type lu -elas_ksp_rtol 5e-3
1302: test:
1303: suffix: 3
1304: nsize: 4
1305: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,10,1000,100 -brick_nu 0.4,0.2,0.3,0.1 -brick_span 3 -elas_pc_type asm -elas_sub_pc_type lu -elas_ksp_rtol 5e-3
1307: test:
1308: suffix: 4
1309: nsize: 4
1310: args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type unpreconditioned -mx 40 -my 40 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_mg_levels_ksp_type chebyshev -elas_pc_type ml -elas_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -elas_mg_levels_pc_type pbjacobi -elas_mg_levels_ksp_max_it 3 -use_nonsymbc -elas_pc_ml_nullspace user
1311: requires: ml
1313: test:
1314: suffix: 5
1315: nsize: 3
1316: args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_pc_type gamg -elas_mg_levels_ksp_type chebyshev -elas_mg_levels_ksp_max_it 1 -elas_mg_levels_ksp_chebyshev_esteig 0.2,1.1 -elas_mg_levels_pc_type jacobi -elas_pc_gamg_esteig_ksp_type cg
1318: test:
1319: suffix: 6
1320: nsize: 4
1321: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipegcr -elas_pc_type asm -elas_sub_pc_type lu
1323: test:
1324: suffix: 7
1325: nsize: 4
1326: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipegcr -elas_pc_type asm -elas_sub_pc_type ksp -elas_sub_ksp_ksp_type cg -elas_sub_ksp_ksp_max_it 15
1328: test:
1329: suffix: 8
1330: nsize: 4
1331: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipefgmres -elas_pc_type asm -elas_sub_pc_type ksp -elas_sub_ksp_ksp_type cg -elas_sub_ksp_ksp_max_it 15
1333: test:
1334: suffix: hypre_nullspace
1335: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
1336: args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_pc_type hypre -elas_pc_hypre_boomeramg_nodal_coarsen 6 -elas_pc_hypre_boomeramg_vec_interp_variant 3 -elas_pc_hypre_boomeramg_interp_type ext+i -elas_ksp_view
1338: test:
1339: nsize: 4
1340: suffix: bddc
1341: args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic
1343: test:
1344: nsize: 4
1345: suffix: bddc_unsym
1346: args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic -use_nonsymbc -elas_pc_bddc_symmetric 0
1348: test:
1349: nsize: 4
1350: suffix: bddc_unsym_deluxe
1351: args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic -use_nonsymbc -elas_pc_bddc_symmetric 0 -elas_pc_bddc_use_deluxe_scaling -elas_sub_schurs_symmetric 0
1353: test:
1354: nsize: 4
1355: suffix: fetidp_unsym_deluxe
1356: args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type fetidp -elas_fetidp_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_fetidp_bddc_pc_bddc_monolithic -use_nonsymbc -elas_fetidp_bddc_pc_bddc_use_deluxe_scaling -elas_fetidp_bddc_sub_schurs_symmetric 0 -elas_fetidp_bddc_pc_bddc_deluxe_singlemat
1358: test:
1359: nsize: 4
1360: suffix: bddc_layerjump
1361: args: -mx 40 -my 40 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_ksp_norm_type natural
1363: test:
1364: nsize: 4
1365: suffix: bddc_subdomainjump
1366: args: -mx 40 -my 40 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,1000 -brick_nu 0.4,0.2 -brick_span 20 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_pc_is_use_stiffness_scaling -elas_ksp_norm_type natural
1368: test:
1369: nsize: 9
1370: suffix: bddc_subdomainjump_deluxe
1371: args: -mx 30 -my 30 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,1000 -brick_nu 0.4,0.2 -brick_span 10 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_pc_bddc_use_deluxe_scaling -elas_ksp_norm_type natural -elas_pc_bddc_schur_layers 1
1372: TEST*/