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;
209: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
211: 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: PetscMalloc1(cpu_x, &LX);
217: PetscMalloc1(cpu_y, &LY);
219: DMDAGetElementsSizes(da, &local_mx, &local_my, NULL);
220: VecCreate(PETSC_COMM_WORLD, &vlx);
221: VecSetSizes(vlx, PETSC_DECIDE, cpu_x);
222: VecSetFromOptions(vlx);
224: VecCreate(PETSC_COMM_WORLD, &vly);
225: VecSetSizes(vly, PETSC_DECIDE, cpu_y);
226: VecSetFromOptions(vly);
228: VecSetValue(vlx, proc_I, (PetscScalar)(local_mx + 1.0e-9), INSERT_VALUES);
229: VecSetValue(vly, proc_J, (PetscScalar)(local_my + 1.0e-9), INSERT_VALUES);
230: VecAssemblyBegin(vlx); VecAssemblyEnd(vlx);
231: VecAssemblyBegin(vly); VecAssemblyEnd(vly);
233: VecScatterCreateToAll(vlx, &ctx, &V_SEQ);
234: VecScatterBegin(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD);
235: VecScatterEnd(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD);
236: VecGetArray(V_SEQ, &_a);
237: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
238: VecRestoreArray(V_SEQ, &_a);
239: VecScatterDestroy(&ctx);
240: VecDestroy(&V_SEQ);
242: VecScatterCreateToAll(vly, &ctx, &V_SEQ);
243: VecScatterBegin(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD);
244: VecScatterEnd(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD);
245: VecGetArray(V_SEQ, &_a);
246: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
247: VecRestoreArray(V_SEQ, &_a);
248: VecScatterDestroy(&ctx);
249: VecDestroy(&V_SEQ);
251: *_lx = LX;
252: *_ly = LY;
254: VecDestroy(&vlx);
255: VecDestroy(&vly);
256: return 0;
257: }
259: static PetscErrorCode DMDACoordViewGnuplot2d(DM da, const char prefix[])
260: {
261: DM cda;
262: Vec coords;
263: DMDACoor2d **_coords;
264: PetscInt si, sj, nx, ny, i, j;
265: FILE *fp;
266: char fname[PETSC_MAX_PATH_LEN];
267: PetscMPIInt rank;
270: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
271: PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank);
272: PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp);
274: PetscFPrintf(PETSC_COMM_SELF, fp, "### Element geometry for processor %1.4d ### \n", rank);
276: DMGetCoordinateDM(da, &cda);
277: DMGetCoordinatesLocal(da, &coords);
278: DMDAVecGetArray(cda, coords, &_coords);
279: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
280: for (j = sj; j < sj + ny - 1; j++) {
281: for (i = si; i < si + nx - 1; i++) {
282: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y));
283: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i].x), (double)PetscRealPart(_coords[j + 1][i].y));
284: 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));
285: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i + 1].x), (double)PetscRealPart(_coords[j][i + 1].y));
286: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n\n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y));
287: }
288: }
289: DMDAVecRestoreArray(cda, coords, &_coords);
291: PetscFClose(PETSC_COMM_SELF, fp);
292: return 0;
293: }
295: static PetscErrorCode DMDAViewGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
296: {
297: DM cda;
298: Vec coords, local_fields;
299: DMDACoor2d **_coords;
300: FILE *fp;
301: char fname[PETSC_MAX_PATH_LEN];
302: const char *field_name;
303: PetscMPIInt rank;
304: PetscInt si, sj, nx, ny, i, j;
305: PetscInt n_dofs, d;
306: PetscScalar *_fields;
309: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
310: PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank);
311: PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp);
314: PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank);
315: DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);
316: PetscFPrintf(PETSC_COMM_SELF, fp, "### x y ");
317: for (d = 0; d < n_dofs; d++) {
318: DMDAGetFieldName(da, d, &field_name);
319: PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name);
320: }
321: PetscFPrintf(PETSC_COMM_SELF, fp, "###\n");
323: DMGetCoordinateDM(da, &cda);
324: DMGetCoordinatesLocal(da, &coords);
325: DMDAVecGetArray(cda, coords, &_coords);
326: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
328: DMCreateLocalVector(da, &local_fields);
329: DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields);
330: DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields);
331: VecGetArray(local_fields, &_fields);
333: for (j = sj; j < sj + ny; j++) {
334: for (i = si; i < si + nx; i++) {
335: PetscScalar coord_x, coord_y;
336: PetscScalar field_d;
338: coord_x = _coords[j][i].x;
339: coord_y = _coords[j][i].y;
341: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y));
342: for (d = 0; d < n_dofs; d++) {
343: field_d = _fields[n_dofs * ((i - si) + (j - sj) * (nx)) + d];
344: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e ", (double)PetscRealPart(field_d));
345: }
346: PetscFPrintf(PETSC_COMM_SELF, fp, "\n");
347: }
348: }
349: VecRestoreArray(local_fields, &_fields);
350: VecDestroy(&local_fields);
352: DMDAVecRestoreArray(cda, coords, &_coords);
354: PetscFClose(PETSC_COMM_SELF, fp);
355: return 0;
356: }
358: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
359: {
360: DM cda;
361: Vec local_fields;
362: FILE *fp;
363: char fname[PETSC_MAX_PATH_LEN];
364: const char *field_name;
365: PetscMPIInt rank;
366: PetscInt si, sj, nx, ny, i, j, p;
367: PetscInt n_dofs, d;
368: GaussPointCoefficients **_coefficients;
371: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
372: PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank);
373: PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp);
376: PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank);
377: DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);
378: PetscFPrintf(PETSC_COMM_SELF, fp, "### x y ");
379: for (d = 0; d < n_dofs; d++) {
380: DMDAGetFieldName(da, d, &field_name);
381: PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name);
382: }
383: PetscFPrintf(PETSC_COMM_SELF, fp, "###\n");
385: DMGetCoordinateDM(da, &cda);
386: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
388: DMCreateLocalVector(da, &local_fields);
389: DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields);
390: DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields);
391: DMDAVecGetArray(da, local_fields, &_coefficients);
393: for (j = sj; j < sj + ny; j++) {
394: for (i = si; i < si + nx; i++) {
395: PetscScalar coord_x, coord_y;
397: for (p = 0; p < GAUSS_POINTS; p++) {
398: coord_x = _coefficients[j][i].gp_coords[2 * p];
399: coord_y = _coefficients[j][i].gp_coords[2 * p + 1];
401: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y));
403: 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]),
404: (double)PetscRealPart(_coefficients[j][i].fy[p])));
405: }
406: }
407: }
408: DMDAVecRestoreArray(da, local_fields, &_coefficients);
409: VecDestroy(&local_fields);
411: PetscFClose(PETSC_COMM_SELF, fp);
412: return 0;
413: }
415: static void FormStressOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar E[], PetscScalar nu[])
416: {
417: PetscInt ngp;
418: PetscScalar gp_xi[GAUSS_POINTS][2];
419: PetscScalar gp_weight[GAUSS_POINTS];
420: PetscInt p, i, j, k, l;
421: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
422: PetscScalar J_p;
423: PetscScalar B[3][U_DOFS * NODES_PER_EL];
424: PetscScalar prop_E, prop_nu, factor, constit_D[3][3];
426: /* define quadrature rule */
427: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
429: /* evaluate integral */
430: for (p = 0; p < ngp; p++) {
431: ConstructQ12D_GNi(gp_xi[p], GNi_p);
432: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
434: for (i = 0; i < NODES_PER_EL; i++) {
435: PetscScalar d_dx_i = GNx_p[0][i];
436: PetscScalar d_dy_i = GNx_p[1][i];
438: B[0][2 * i] = d_dx_i;
439: B[0][2 * i + 1] = 0.0;
440: B[1][2 * i] = 0.0;
441: B[1][2 * i + 1] = d_dy_i;
442: B[2][2 * i] = d_dy_i;
443: B[2][2 * i + 1] = d_dx_i;
444: }
446: /* form D for the quadrature point */
447: prop_E = E[p];
448: prop_nu = nu[p];
449: factor = prop_E / ((1.0 + prop_nu) * (1.0 - 2.0 * prop_nu));
450: constit_D[0][0] = 1.0 - prop_nu;
451: constit_D[0][1] = prop_nu;
452: constit_D[0][2] = 0.0;
453: constit_D[1][0] = prop_nu;
454: constit_D[1][1] = 1.0 - prop_nu;
455: constit_D[1][2] = 0.0;
456: constit_D[2][0] = 0.0;
457: constit_D[2][1] = 0.0;
458: constit_D[2][2] = 0.5 * (1.0 - 2.0 * prop_nu);
459: for (i = 0; i < 3; i++) {
460: for (j = 0; j < 3; j++) constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
461: }
463: /* form Bt tildeD B */
464: /*
465: Ke_ij = Bt_ik . D_kl . B_lj
466: = B_ki . D_kl . B_lj
467: */
468: for (i = 0; i < 8; i++) {
469: for (j = 0; j < 8; j++) {
470: for (k = 0; k < 3; k++) {
471: for (l = 0; l < 3; l++) Ke[8 * i + j] = Ke[8 * i + j] + B[k][i] * constit_D[k][l] * B[l][j];
472: }
473: }
474: }
476: } /* end quadrature */
477: }
479: static void FormMomentumRhsQ1(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
480: {
481: PetscInt ngp;
482: PetscScalar gp_xi[GAUSS_POINTS][2];
483: PetscScalar gp_weight[GAUSS_POINTS];
484: PetscInt p, i;
485: PetscScalar Ni_p[NODES_PER_EL];
486: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
487: PetscScalar J_p, fac;
489: /* define quadrature rule */
490: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
492: /* evaluate integral */
493: for (p = 0; p < ngp; p++) {
494: ConstructQ12D_Ni(gp_xi[p], Ni_p);
495: ConstructQ12D_GNi(gp_xi[p], GNi_p);
496: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
497: fac = gp_weight[p] * J_p;
499: for (i = 0; i < NODES_PER_EL; i++) {
500: Fe[NSD * i] += fac * Ni_p[i] * fx[p];
501: Fe[NSD * i + 1] += fac * Ni_p[i] * fy[p];
502: }
503: }
504: }
506: /*
507: i,j are the element indices
508: The unknown is a vector quantity.
509: The s[].c is used to indicate the degree of freedom.
510: */
511: static PetscErrorCode DMDAGetElementEqnums_u(MatStencil s_u[], PetscInt i, PetscInt j)
512: {
514: /* displacement */
515: /* node 0 */
516: s_u[0].i = i;
517: s_u[0].j = j;
518: s_u[0].c = 0; /* Ux0 */
519: s_u[1].i = i;
520: s_u[1].j = j;
521: s_u[1].c = 1; /* Uy0 */
523: /* node 1 */
524: s_u[2].i = i;
525: s_u[2].j = j + 1;
526: s_u[2].c = 0; /* Ux1 */
527: s_u[3].i = i;
528: s_u[3].j = j + 1;
529: s_u[3].c = 1; /* Uy1 */
531: /* node 2 */
532: s_u[4].i = i + 1;
533: s_u[4].j = j + 1;
534: s_u[4].c = 0; /* Ux2 */
535: s_u[5].i = i + 1;
536: s_u[5].j = j + 1;
537: s_u[5].c = 1; /* Uy2 */
539: /* node 3 */
540: s_u[6].i = i + 1;
541: s_u[6].j = j;
542: s_u[6].c = 0; /* Ux3 */
543: s_u[7].i = i + 1;
544: s_u[7].j = j;
545: s_u[7].c = 1; /* Uy3 */
546: return 0;
547: }
549: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords, PetscInt ei, PetscInt ej, PetscScalar el_coords[])
550: {
552: /* get coords for the element */
553: el_coords[NSD * 0 + 0] = _coords[ej][ei].x;
554: el_coords[NSD * 0 + 1] = _coords[ej][ei].y;
555: el_coords[NSD * 1 + 0] = _coords[ej + 1][ei].x;
556: el_coords[NSD * 1 + 1] = _coords[ej + 1][ei].y;
557: el_coords[NSD * 2 + 0] = _coords[ej + 1][ei + 1].x;
558: el_coords[NSD * 2 + 1] = _coords[ej + 1][ei + 1].y;
559: el_coords[NSD * 3 + 0] = _coords[ej][ei + 1].x;
560: el_coords[NSD * 3 + 1] = _coords[ej][ei + 1].y;
561: return 0;
562: }
564: static PetscErrorCode AssembleA_Elasticity(Mat A, DM elas_da, DM properties_da, Vec properties)
565: {
566: DM cda;
567: Vec coords;
568: DMDACoor2d **_coords;
569: MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
570: PetscInt sex, sey, mx, my;
571: PetscInt ei, ej;
572: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
573: PetscScalar el_coords[NODES_PER_EL * NSD];
574: Vec local_properties;
575: GaussPointCoefficients **props;
576: PetscScalar *prop_E, *prop_nu;
579: /* setup for coords */
580: DMGetCoordinateDM(elas_da, &cda);
581: DMGetCoordinatesLocal(elas_da, &coords);
582: DMDAVecGetArray(cda, coords, &_coords);
584: /* setup for coefficients */
585: DMCreateLocalVector(properties_da, &local_properties);
586: DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties);
587: DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties);
588: DMDAVecGetArray(properties_da, local_properties, &props);
590: DMDAGetElementsCorners(elas_da, &sex, &sey, 0);
591: DMDAGetElementsSizes(elas_da, &mx, &my, 0);
592: for (ej = sey; ej < sey + my; ej++) {
593: for (ei = sex; ei < sex + mx; ei++) {
594: /* get coords for the element */
595: GetElementCoords(_coords, ei, ej, el_coords);
597: /* get coefficients for the element */
598: prop_E = props[ej][ei].E;
599: prop_nu = props[ej][ei].nu;
601: /* initialise element stiffness matrix */
602: PetscMemzero(Ae, sizeof(Ae));
604: /* form element stiffness matrix */
605: FormStressOperatorQ1(Ae, el_coords, prop_E, prop_nu);
607: /* insert element matrix into global matrix */
608: DMDAGetElementEqnums_u(u_eqn, ei, ej);
609: MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES);
610: }
611: }
612: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
613: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
615: DMDAVecRestoreArray(cda, coords, &_coords);
617: DMDAVecRestoreArray(properties_da, local_properties, &props);
618: VecDestroy(&local_properties);
619: return 0;
620: }
622: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F, MatStencil u_eqn[], PetscScalar Fe_u[])
623: {
624: PetscInt n;
627: for (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: return 0;
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: PetscInt ei, ej;
642: PetscScalar Fe[NODES_PER_EL * U_DOFS];
643: PetscScalar el_coords[NODES_PER_EL * NSD];
644: Vec local_properties;
645: GaussPointCoefficients **props;
646: PetscScalar *prop_fx, *prop_fy;
647: Vec local_F;
648: ElasticityDOF **ff;
651: /* setup for coords */
652: DMGetCoordinateDM(elas_da, &cda);
653: DMGetCoordinatesLocal(elas_da, &coords);
654: DMDAVecGetArray(cda, coords, &_coords);
656: /* setup for coefficients */
657: DMGetLocalVector(properties_da, &local_properties);
658: DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties);
659: DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties);
660: DMDAVecGetArray(properties_da, local_properties, &props);
662: /* get access to the vector */
663: DMGetLocalVector(elas_da, &local_F);
664: VecZeroEntries(local_F);
665: DMDAVecGetArray(elas_da, local_F, &ff);
667: DMDAGetElementsCorners(elas_da, &sex, &sey, 0);
668: DMDAGetElementsSizes(elas_da, &mx, &my, 0);
669: for (ej = sey; ej < sey + my; ej++) {
670: for (ei = sex; ei < sex + mx; ei++) {
671: /* get coords for the element */
672: GetElementCoords(_coords, ei, ej, el_coords);
674: /* get coefficients for the element */
675: prop_fx = props[ej][ei].fx;
676: prop_fy = props[ej][ei].fy;
678: /* initialise element stiffness matrix */
679: PetscMemzero(Fe, sizeof(Fe));
681: /* form element stiffness matrix */
682: FormMomentumRhsQ1(Fe, el_coords, prop_fx, prop_fy);
684: /* insert element matrix into global matrix */
685: DMDAGetElementEqnums_u(u_eqn, ei, ej);
687: DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, Fe);
688: }
689: }
691: DMDAVecRestoreArray(elas_da, local_F, &ff);
692: DMLocalToGlobalBegin(elas_da, local_F, ADD_VALUES, F);
693: DMLocalToGlobalEnd(elas_da, local_F, ADD_VALUES, F);
694: DMRestoreLocalVector(elas_da, &local_F);
696: DMDAVecRestoreArray(cda, coords, &_coords);
698: DMDAVecRestoreArray(properties_da, local_properties, &props);
699: DMRestoreLocalVector(properties_da, &local_properties);
700: return 0;
701: }
703: static PetscErrorCode solve_elasticity_2d(PetscInt mx, PetscInt my)
704: {
705: DM elas_da, da_prop;
706: PetscInt u_dof, dof, stencil_width;
707: Mat A;
708: PetscInt mxl, myl;
709: DM prop_cda, vel_cda;
710: Vec prop_coords, vel_coords;
711: PetscInt si, sj, nx, ny, i, j, p;
712: Vec f, X;
713: PetscInt prop_dof, prop_stencil_width;
714: Vec properties, l_properties;
715: MatNullSpace matnull;
716: PetscReal dx, dy;
717: PetscInt M, N;
718: DMDACoor2d **_prop_coords, **_vel_coords;
719: GaussPointCoefficients **element_props;
720: KSP ksp_E;
721: PetscInt coefficient_structure = 0;
722: PetscInt cpu_x, cpu_y, *lx = NULL, *ly = NULL;
723: PetscBool use_gp_coords = PETSC_FALSE;
724: PetscBool use_nonsymbc = PETSC_FALSE;
725: PetscBool no_view = PETSC_FALSE;
726: PetscBool flg;
729: /* Generate the da for velocity and pressure */
730: /*
731: We use Q1 elements for the temperature.
732: FEM has a 9-point stencil (BOX) or connectivity pattern
733: Num nodes in each direction is mx+1, my+1
734: */
735: u_dof = U_DOFS; /* Vx, Vy - velocities */
736: dof = u_dof;
737: stencil_width = 1;
738: 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);
740: DMSetMatType(elas_da, MATAIJ);
741: DMSetFromOptions(elas_da);
742: DMSetUp(elas_da);
744: DMDASetFieldName(elas_da, 0, "Ux");
745: DMDASetFieldName(elas_da, 1, "Uy");
747: /* unit box [0,1] x [0,1] */
748: DMDASetUniformCoordinates(elas_da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
750: /* Generate element properties, we will assume all material properties are constant over the element */
751: /* local number of elements */
752: DMDAGetElementsSizes(elas_da, &mxl, &myl, NULL);
754: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
755: DMDAGetInfo(elas_da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0);
756: DMDAGetElementOwnershipRanges2d(elas_da, &lx, &ly);
758: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
759: prop_stencil_width = 0;
760: 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);
761: DMSetFromOptions(da_prop);
762: DMSetUp(da_prop);
764: PetscFree(lx);
765: PetscFree(ly);
767: /* define centroid positions */
768: DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
769: dx = 1.0 / ((PetscReal)(M));
770: dy = 1.0 / ((PetscReal)(N));
772: 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);
774: /* define coefficients */
775: PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL);
777: DMCreateGlobalVector(da_prop, &properties);
778: DMCreateLocalVector(da_prop, &l_properties);
779: DMDAVecGetArray(da_prop, l_properties, &element_props);
781: DMGetCoordinateDM(da_prop, &prop_cda);
782: DMGetCoordinatesLocal(da_prop, &prop_coords);
783: DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords);
785: DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0);
787: DMGetCoordinateDM(elas_da, &vel_cda);
788: DMGetCoordinatesLocal(elas_da, &vel_coords);
789: DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords);
791: /* interpolate the coordinates */
792: for (j = sj; j < sj + ny; j++) {
793: for (i = si; i < si + nx; i++) {
794: PetscInt ngp;
795: PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
796: PetscScalar el_coords[8];
798: GetElementCoords(_vel_coords, i, j, el_coords);
799: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
801: for (p = 0; p < GAUSS_POINTS; p++) {
802: PetscScalar gp_x, gp_y;
803: PetscInt n;
804: PetscScalar xi_p[2], Ni_p[4];
806: xi_p[0] = gp_xi[p][0];
807: xi_p[1] = gp_xi[p][1];
808: ConstructQ12D_Ni(xi_p, Ni_p);
810: gp_x = 0.0;
811: gp_y = 0.0;
812: for (n = 0; n < NODES_PER_EL; n++) {
813: gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
814: gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
815: }
816: element_props[j][i].gp_coords[2 * p] = gp_x;
817: element_props[j][i].gp_coords[2 * p + 1] = gp_y;
818: }
819: }
820: }
822: /* define the coefficients */
823: PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, &flg);
825: for (j = sj; j < sj + ny; j++) {
826: for (i = si; i < si + nx; i++) {
827: PetscScalar centroid_x = _prop_coords[j][i].x; /* centroids of cell */
828: PetscScalar centroid_y = _prop_coords[j][i].y;
829: PETSC_UNUSED PetscScalar coord_x, coord_y;
831: if (coefficient_structure == 0) { /* isotropic */
832: PetscScalar opts_E, opts_nu;
834: opts_E = 1.0;
835: opts_nu = 0.33;
836: PetscOptionsGetScalar(NULL, NULL, "-iso_E", &opts_E, &flg);
837: PetscOptionsGetScalar(NULL, NULL, "-iso_nu", &opts_nu, &flg);
839: for (p = 0; p < GAUSS_POINTS; p++) {
840: element_props[j][i].E[p] = opts_E;
841: element_props[j][i].nu[p] = opts_nu;
843: element_props[j][i].fx[p] = 0.0;
844: element_props[j][i].fy[p] = 0.0;
845: }
846: } else if (coefficient_structure == 1) { /* step */
847: PetscScalar opts_E0, opts_nu0, opts_xc;
848: PetscScalar opts_E1, opts_nu1;
850: opts_E0 = opts_E1 = 1.0;
851: opts_nu0 = opts_nu1 = 0.333;
852: opts_xc = 0.5;
853: PetscOptionsGetScalar(NULL, NULL, "-step_E0", &opts_E0, &flg);
854: PetscOptionsGetScalar(NULL, NULL, "-step_nu0", &opts_nu0, &flg);
855: PetscOptionsGetScalar(NULL, NULL, "-step_E1", &opts_E1, &flg);
856: PetscOptionsGetScalar(NULL, NULL, "-step_nu1", &opts_nu1, &flg);
857: PetscOptionsGetScalar(NULL, NULL, "-step_xc", &opts_xc, &flg);
859: for (p = 0; p < GAUSS_POINTS; p++) {
860: coord_x = centroid_x;
861: coord_y = centroid_y;
862: if (use_gp_coords) {
863: coord_x = element_props[j][i].gp_coords[2 * p];
864: coord_y = element_props[j][i].gp_coords[2 * p + 1];
865: }
867: element_props[j][i].E[p] = opts_E0;
868: element_props[j][i].nu[p] = opts_nu0;
869: if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
870: element_props[j][i].E[p] = opts_E1;
871: element_props[j][i].nu[p] = opts_nu1;
872: }
874: element_props[j][i].fx[p] = 0.0;
875: element_props[j][i].fy[p] = 0.0;
876: }
877: } else if (coefficient_structure == 2) { /* brick */
878: PetscReal values_E[10];
879: PetscReal values_nu[10];
880: PetscInt nbricks, maxnbricks;
881: PetscInt index, span;
882: PetscInt jj;
884: flg = PETSC_FALSE;
885: maxnbricks = 10;
886: PetscOptionsGetRealArray(NULL, NULL, "-brick_E", values_E, &maxnbricks, &flg);
887: nbricks = maxnbricks;
890: flg = PETSC_FALSE;
891: maxnbricks = 10;
892: PetscOptionsGetRealArray(NULL, NULL, "-brick_nu", values_nu, &maxnbricks, &flg);
896: span = 1;
897: PetscOptionsGetInt(NULL, NULL, "-brick_span", &span, &flg);
899: /* cycle through the indices so that no two material properties are repeated in lines of x or y */
900: jj = (j / span) % nbricks;
901: index = (jj + i / span) % nbricks;
902: /*printf("j=%d: index = %d \n", j,index); */
904: for (p = 0; p < GAUSS_POINTS; p++) {
905: element_props[j][i].E[p] = values_E[index];
906: element_props[j][i].nu[p] = values_nu[index];
907: }
908: } else if (coefficient_structure == 3) { /* sponge */
909: PetscScalar opts_E0, opts_nu0;
910: PetscScalar opts_E1, opts_nu1;
911: PetscInt opts_t, opts_w;
912: PetscInt ii, jj, ci, cj;
914: opts_E0 = opts_E1 = 1.0;
915: opts_nu0 = opts_nu1 = 0.333;
916: PetscOptionsGetScalar(NULL, NULL, "-sponge_E0", &opts_E0, &flg);
917: PetscOptionsGetScalar(NULL, NULL, "-sponge_nu0", &opts_nu0, &flg);
918: PetscOptionsGetScalar(NULL, NULL, "-sponge_E1", &opts_E1, &flg);
919: PetscOptionsGetScalar(NULL, NULL, "-sponge_nu1", &opts_nu1, &flg);
921: opts_t = opts_w = 1;
922: PetscOptionsGetInt(NULL, NULL, "-sponge_t", &opts_t, &flg);
923: PetscOptionsGetInt(NULL, NULL, "-sponge_w", &opts_w, &flg);
925: ii = (i) / (opts_t + opts_w + opts_t);
926: jj = (j) / (opts_t + opts_w + opts_t);
928: ci = i - ii * (opts_t + opts_w + opts_t);
929: cj = j - jj * (opts_t + opts_w + opts_t);
931: for (p = 0; p < GAUSS_POINTS; p++) {
932: element_props[j][i].E[p] = opts_E0;
933: element_props[j][i].nu[p] = opts_nu0;
934: }
935: if ((ci >= opts_t) && (ci < opts_t + opts_w)) {
936: if ((cj >= opts_t) && (cj < opts_t + opts_w)) {
937: for (p = 0; p < GAUSS_POINTS; p++) {
938: element_props[j][i].E[p] = opts_E1;
939: element_props[j][i].nu[p] = opts_nu1;
940: }
941: }
942: }
943: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
944: }
945: }
946: DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords);
948: DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords);
950: DMDAVecRestoreArray(da_prop, l_properties, &element_props);
951: DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties);
952: DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties);
954: PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL);
955: if (!no_view) {
956: DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for elasticity eqn.", "properties");
957: DMDACoordViewGnuplot2d(elas_da, "mesh");
958: }
960: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
961: DMCreateMatrix(elas_da, &A);
962: DMGetCoordinates(elas_da, &vel_coords);
963: MatNullSpaceCreateRigidBody(vel_coords, &matnull);
964: MatSetNearNullSpace(A, matnull);
965: MatNullSpaceDestroy(&matnull);
966: MatCreateVecs(A, &f, &X);
968: /* assemble A11 */
969: MatZeroEntries(A);
970: VecZeroEntries(f);
972: AssembleA_Elasticity(A, elas_da, da_prop, properties);
973: /* build force vector */
974: AssembleF_Elasticity(f, elas_da, da_prop, properties);
976: KSPCreate(PETSC_COMM_WORLD, &ksp_E);
977: KSPSetOptionsPrefix(ksp_E, "elas_"); /* elasticity */
979: PetscOptionsGetBool(NULL, NULL, "-use_nonsymbc", &use_nonsymbc, &flg);
980: /* solve */
981: if (!use_nonsymbc) {
982: Mat AA;
983: Vec ff, XX;
984: IS is;
985: VecScatter scat;
987: DMDABCApplySymmetricCompression(elas_da, A, f, &is, &AA, &ff);
988: VecDuplicate(ff, &XX);
990: KSPSetOperators(ksp_E, AA, AA);
991: KSPSetFromOptions(ksp_E);
993: KSPSolve(ksp_E, ff, XX);
995: /* push XX back into X */
996: DMDABCApplyCompression(elas_da, NULL, X);
998: VecScatterCreate(XX, NULL, X, is, &scat);
999: VecScatterBegin(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD);
1000: VecScatterEnd(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD);
1001: VecScatterDestroy(&scat);
1003: MatDestroy(&AA);
1004: VecDestroy(&ff);
1005: VecDestroy(&XX);
1006: ISDestroy(&is);
1007: } else {
1008: DMDABCApplyCompression(elas_da, A, f);
1010: KSPSetOperators(ksp_E, A, A);
1011: KSPSetFromOptions(ksp_E);
1013: KSPSolve(ksp_E, f, X);
1014: }
1016: if (!no_view) DMDAViewGnuplot2d(elas_da, X, "Displacement solution for elasticity eqn.", "X");
1017: KSPDestroy(&ksp_E);
1019: VecDestroy(&X);
1020: VecDestroy(&f);
1021: MatDestroy(&A);
1023: DMDestroy(&elas_da);
1024: DMDestroy(&da_prop);
1026: VecDestroy(&properties);
1027: VecDestroy(&l_properties);
1028: return 0;
1029: }
1031: int main(int argc, char **args)
1032: {
1033: PetscInt mx, my;
1036: PetscInitialize(&argc, &args, (char *)0, help);
1037: mx = my = 10;
1038: PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL);
1039: PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL);
1040: solve_elasticity_2d(mx, my);
1041: PetscFinalize();
1042: return 0;
1043: }
1045: /* -------------------------- helpers for boundary conditions -------------------------------- */
1047: static PetscErrorCode BCApply_EAST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
1048: {
1049: DM cda;
1050: Vec coords;
1051: PetscInt si, sj, nx, ny, i, j;
1052: PetscInt M, N;
1053: DMDACoor2d **_coords;
1054: const PetscInt *g_idx;
1055: PetscInt *bc_global_ids;
1056: PetscScalar *bc_vals;
1057: PetscInt nbcs;
1058: PetscInt n_dofs;
1059: ISLocalToGlobalMapping ltogm;
1062: /* enforce bc's */
1063: DMGetLocalToGlobalMapping(da, <ogm);
1064: ISLocalToGlobalMappingGetIndices(ltogm, &g_idx);
1066: DMGetCoordinateDM(da, &cda);
1067: DMGetCoordinatesLocal(da, &coords);
1068: DMDAVecGetArray(cda, coords, &_coords);
1069: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
1070: DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);
1072: /* --- */
1074: PetscMalloc1(ny * n_dofs, &bc_global_ids);
1075: PetscMalloc1(ny * n_dofs, &bc_vals);
1077: /* init the entries to -1 so VecSetValues will ignore them */
1078: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1080: i = nx - 1;
1081: for (j = 0; j < ny; j++) {
1082: PetscInt local_id;
1083: PETSC_UNUSED PetscScalar coordx, coordy;
1085: local_id = i + j * nx;
1087: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1089: coordx = _coords[j + sj][i + si].x;
1090: coordy = _coords[j + sj][i + si].y;
1092: bc_vals[j] = bc_val;
1093: }
1094: ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1095: nbcs = 0;
1096: if ((si + nx) == (M)) nbcs = ny;
1098: if (b) {
1099: VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES);
1100: VecAssemblyBegin(b);
1101: VecAssemblyEnd(b);
1102: }
1103: if (A) MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0);
1105: PetscFree(bc_vals);
1106: PetscFree(bc_global_ids);
1108: DMDAVecRestoreArray(cda, coords, &_coords);
1109: return 0;
1110: }
1112: static PetscErrorCode BCApply_WEST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
1113: {
1114: DM cda;
1115: Vec coords;
1116: PetscInt si, sj, nx, ny, i, j;
1117: PetscInt M, N;
1118: DMDACoor2d **_coords;
1119: const PetscInt *g_idx;
1120: PetscInt *bc_global_ids;
1121: PetscScalar *bc_vals;
1122: PetscInt nbcs;
1123: PetscInt n_dofs;
1124: ISLocalToGlobalMapping ltogm;
1127: /* enforce bc's */
1128: DMGetLocalToGlobalMapping(da, <ogm);
1129: ISLocalToGlobalMappingGetIndices(ltogm, &g_idx);
1131: DMGetCoordinateDM(da, &cda);
1132: DMGetCoordinatesLocal(da, &coords);
1133: DMDAVecGetArray(cda, coords, &_coords);
1134: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
1135: DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);
1137: /* --- */
1139: PetscMalloc1(ny * n_dofs, &bc_global_ids);
1140: PetscMalloc1(ny * n_dofs, &bc_vals);
1142: /* init the entries to -1 so VecSetValues will ignore them */
1143: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1145: i = 0;
1146: for (j = 0; j < ny; j++) {
1147: PetscInt local_id;
1148: PETSC_UNUSED PetscScalar coordx, coordy;
1150: local_id = i + j * nx;
1152: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1154: coordx = _coords[j + sj][i + si].x;
1155: coordy = _coords[j + sj][i + si].y;
1157: bc_vals[j] = bc_val;
1158: }
1159: ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1160: nbcs = 0;
1161: if (si == 0) nbcs = ny;
1163: if (b) {
1164: VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES);
1165: VecAssemblyBegin(b);
1166: VecAssemblyEnd(b);
1167: }
1168: if (A) MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0);
1170: PetscFree(bc_vals);
1171: PetscFree(bc_global_ids);
1173: DMDAVecRestoreArray(cda, coords, &_coords);
1174: return 0;
1175: }
1177: static PetscErrorCode DMDABCApplyCompression(DM elas_da, Mat A, Vec f)
1178: {
1180: BCApply_EAST(elas_da, 0, -1.0, A, f);
1181: BCApply_EAST(elas_da, 1, 0.0, A, f);
1182: BCApply_WEST(elas_da, 0, 1.0, A, f);
1183: BCApply_WEST(elas_da, 1, 0.0, A, f);
1184: return 0;
1185: }
1187: static PetscErrorCode Orthogonalize(PetscInt n, Vec *vecs)
1188: {
1189: PetscInt i, j;
1190: PetscScalar dot;
1192: for (i = 0; i < n; i++) {
1193: VecNormalize(vecs[i], NULL);
1194: for (j = i + 1; j < n; j++) {
1195: VecDot(vecs[i], vecs[j], &dot);
1196: VecAXPY(vecs[j], -dot, vecs[i]);
1197: }
1198: }
1199: return 0;
1200: }
1202: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da, Mat A, Vec f, IS *dofs, Mat *AA, Vec *ff)
1203: {
1204: PetscInt start, end, m;
1205: PetscInt *unconstrained;
1206: PetscInt cnt, i;
1207: Vec x;
1208: PetscScalar *_x;
1209: IS is;
1210: VecScatter scat;
1213: /* push bc's into f and A */
1214: VecDuplicate(f, &x);
1215: BCApply_EAST(elas_da, 0, -1.0, A, x);
1216: BCApply_EAST(elas_da, 1, 0.0, A, x);
1217: BCApply_WEST(elas_da, 0, 1.0, A, x);
1218: BCApply_WEST(elas_da, 1, 0.0, A, x);
1220: /* define which dofs are not constrained */
1221: VecGetLocalSize(x, &m);
1222: PetscMalloc1(m, &unconstrained);
1223: VecGetOwnershipRange(x, &start, &end);
1224: VecGetArray(x, &_x);
1225: cnt = 0;
1226: for (i = 0; i < m; i += 2) {
1227: PetscReal val1, val2;
1229: val1 = PetscRealPart(_x[i]);
1230: val2 = PetscRealPart(_x[i + 1]);
1231: if (PetscAbs(val1) < 0.1 && PetscAbs(val2) < 0.1) {
1232: unconstrained[cnt] = start + i;
1233: cnt++;
1234: unconstrained[cnt] = start + i + 1;
1235: cnt++;
1236: }
1237: }
1238: VecRestoreArray(x, &_x);
1240: ISCreateGeneral(PETSC_COMM_WORLD, cnt, unconstrained, PETSC_COPY_VALUES, &is);
1241: PetscFree(unconstrained);
1242: ISSetBlockSize(is, 2);
1244: /* define correction for dirichlet in the rhs */
1245: MatMult(A, x, f);
1246: VecScale(f, -1.0);
1248: /* get new matrix */
1249: MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, AA);
1250: /* get new vector */
1251: MatCreateVecs(*AA, NULL, ff);
1253: VecScatterCreate(f, is, *ff, NULL, &scat);
1254: VecScatterBegin(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD);
1255: VecScatterEnd(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD);
1257: { /* Constrain near-null space */
1258: PetscInt nvecs;
1259: const Vec *vecs;
1260: Vec *uvecs;
1261: PetscBool has_const;
1262: MatNullSpace mnull, unull;
1264: MatGetNearNullSpace(A, &mnull);
1265: MatNullSpaceGetVecs(mnull, &has_const, &nvecs, &vecs);
1266: VecDuplicateVecs(*ff, nvecs, &uvecs);
1267: for (i = 0; i < nvecs; i++) {
1268: VecScatterBegin(scat, vecs[i], uvecs[i], INSERT_VALUES, SCATTER_FORWARD);
1269: VecScatterEnd(scat, vecs[i], uvecs[i], INSERT_VALUES, SCATTER_FORWARD);
1270: }
1271: Orthogonalize(nvecs, uvecs);
1272: MatNullSpaceCreate(PetscObjectComm((PetscObject)A), PETSC_FALSE, nvecs, uvecs, &unull);
1273: MatSetNearNullSpace(*AA, unull);
1274: MatNullSpaceDestroy(&unull);
1275: VecDestroyVecs(nvecs, &uvecs);
1276: }
1278: VecScatterDestroy(&scat);
1280: *dofs = is;
1281: VecDestroy(&x);
1282: return 0;
1283: }
1285: /*TEST
1287: build:
1288: requires: !complex !single
1290: test:
1291: 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
1292: output_file: output/ex49_1.out
1294: test:
1295: suffix: 2
1296: nsize: 4
1297: 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
1299: test:
1300: suffix: 3
1301: nsize: 4
1302: 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
1304: test:
1305: suffix: 4
1306: nsize: 4
1307: 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
1308: requires: ml
1310: test:
1311: suffix: 5
1312: nsize: 3
1313: 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
1315: test:
1316: suffix: 6
1317: nsize: 4
1318: 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
1320: test:
1321: suffix: 7
1322: nsize: 4
1323: 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
1325: test:
1326: suffix: 8
1327: nsize: 4
1328: 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
1330: test:
1331: suffix: hypre_nullspace
1332: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
1333: 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
1335: test:
1336: nsize: 4
1337: suffix: bddc
1338: 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
1340: test:
1341: nsize: 4
1342: suffix: bddc_unsym
1343: 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
1345: test:
1346: nsize: 4
1347: suffix: bddc_unsym_deluxe
1348: 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
1350: test:
1351: nsize: 4
1352: suffix: fetidp_unsym_deluxe
1353: 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
1355: test:
1356: nsize: 4
1357: suffix: bddc_layerjump
1358: 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
1360: test:
1361: nsize: 4
1362: suffix: bddc_subdomainjump
1363: 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
1365: test:
1366: nsize: 9
1367: suffix: bddc_subdomainjump_deluxe
1368: 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
1369: TEST*/