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 evaulated 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);
78: #define NSD 2 /* number of spatial dimensions */
79: #define NODES_PER_EL 4 /* nodes per element */
80: #define U_DOFS 2 /* degrees of freedom per displacement node */
81: #define GAUSS_POINTS 4
83: /* cell based evaluation */
84: typedef struct {
85: PetscScalar E,nu,fx,fy;
86: } Coefficients;
88: /* Gauss point based evaluation 8+4+4+4 = 20 */
89: typedef struct {
90: PetscScalar gp_coords[2*GAUSS_POINTS];
91: PetscScalar E[GAUSS_POINTS];
92: PetscScalar nu[GAUSS_POINTS];
93: PetscScalar fx[GAUSS_POINTS];
94: PetscScalar fy[GAUSS_POINTS];
95: } GaussPointCoefficients;
97: typedef struct {
98: PetscScalar ux_dof;
99: PetscScalar uy_dof;
100: } ElasticityDOF;
103: /*
105: D = E/((1+nu)(1-2nu)) * [ 1-nu nu 0 ]
106: [ nu 1-nu 0 ]
107: [ 0 0 0.5*(1-2nu) ]
109: B = [ d_dx 0 ]
110: [ 0 d_dy ]
111: [ d_dy d_dx ]
113: */
115: /* FEM routines */
116: /*
117: Element: Local basis function ordering
118: 1-----2
119: | |
120: | |
121: 0-----3
122: */
123: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
124: {
125: PetscScalar xi = _xi[0];
126: PetscScalar eta = _xi[1];
128: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
129: Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
130: Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
131: Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
132: }
134: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
135: {
136: PetscScalar xi = _xi[0];
137: PetscScalar eta = _xi[1];
139: GNi[0][0] = -0.25*(1.0-eta);
140: GNi[0][1] = -0.25*(1.0+eta);
141: GNi[0][2] = 0.25*(1.0+eta);
142: GNi[0][3] = 0.25*(1.0-eta);
144: GNi[1][0] = -0.25*(1.0-xi);
145: GNi[1][1] = 0.25*(1.0-xi);
146: GNi[1][2] = 0.25*(1.0+xi);
147: GNi[1][3] = -0.25*(1.0+xi);
148: }
150: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
151: {
152: PetscScalar J00,J01,J10,J11,J;
153: PetscScalar iJ00,iJ01,iJ10,iJ11;
154: PetscInt i;
156: J00 = J01 = J10 = J11 = 0.0;
157: for (i = 0; i < NODES_PER_EL; i++) {
158: PetscScalar cx = coords[2*i+0];
159: PetscScalar cy = coords[2*i+1];
161: J00 = J00+GNi[0][i]*cx; /* J_xx = dx/dxi */
162: J01 = J01+GNi[0][i]*cy; /* J_xy = dy/dxi */
163: J10 = J10+GNi[1][i]*cx; /* J_yx = dx/deta */
164: J11 = J11+GNi[1][i]*cy; /* J_yy = dy/deta */
165: }
166: J = (J00*J11)-(J01*J10);
168: iJ00 = J11/J;
169: iJ01 = -J01/J;
170: iJ10 = -J10/J;
171: iJ11 = J00/J;
174: for (i = 0; i < NODES_PER_EL; i++) {
175: GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
176: GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
177: }
179: if (det_J) *det_J = J;
180: }
182: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
183: {
184: *ngp = 4;
185: gp_xi[0][0] = -0.57735026919;gp_xi[0][1] = -0.57735026919;
186: gp_xi[1][0] = -0.57735026919;gp_xi[1][1] = 0.57735026919;
187: gp_xi[2][0] = 0.57735026919;gp_xi[2][1] = 0.57735026919;
188: gp_xi[3][0] = 0.57735026919;gp_xi[3][1] = -0.57735026919;
189: gp_weight[0] = 1.0;
190: gp_weight[1] = 1.0;
191: gp_weight[2] = 1.0;
192: gp_weight[3] = 1.0;
193: }
195: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
196: {
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;
271: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
272: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
273: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
274: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
275: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
277: DMGetCoordinateDM(da,&cda);
278: DMGetCoordinatesLocal(da,&coords);
279: DMDAVecGetArray(cda,coords,&_coords);
280: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
281: for (j = sj; j < sj+ny-1; j++) {
282: for (i = si; i < si+nx-1; i++) {
283: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
284: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i].x),PetscRealPart(_coords[j+1][i].y));
285: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i+1].x),PetscRealPart(_coords[j+1][i+1].y));
286: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i+1].x),PetscRealPart(_coords[j][i+1].y));
287: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
288: }
289: }
290: DMDAVecRestoreArray(cda,coords,&_coords);
292: PetscFClose(PETSC_COMM_SELF,fp);
293: return(0);
294: }
296: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
297: {
298: DM cda;
299: Vec coords,local_fields;
300: DMDACoor2d **_coords;
301: FILE *fp;
302: char fname[PETSC_MAX_PATH_LEN];
303: const char *field_name;
304: PetscMPIInt rank;
305: PetscInt si,sj,nx,ny,i,j;
306: PetscInt n_dofs,d;
307: PetscScalar *_fields;
311: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
312: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
313: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
314: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
316: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
317: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
318: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
319: for (d = 0; d < n_dofs; d++) {
320: DMDAGetFieldName(da,d,&field_name);
321: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
322: }
323: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
326: DMGetCoordinateDM(da,&cda);
327: DMGetCoordinatesLocal(da,&coords);
328: DMDAVecGetArray(cda,coords,&_coords);
329: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
331: DMCreateLocalVector(da,&local_fields);
332: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
333: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
334: VecGetArray(local_fields,&_fields);
336: for (j = sj; j < sj+ny; j++) {
337: for (i = si; i < si+nx; i++) {
338: PetscScalar coord_x,coord_y;
339: PetscScalar field_d;
341: coord_x = _coords[j][i].x;
342: coord_y = _coords[j][i].y;
344: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
345: for (d = 0; d < n_dofs; d++) {
346: field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
347: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",PetscRealPart(field_d));
348: }
349: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
350: }
351: }
352: VecRestoreArray(local_fields,&_fields);
353: VecDestroy(&local_fields);
355: DMDAVecRestoreArray(cda,coords,&_coords);
357: PetscFClose(PETSC_COMM_SELF,fp);
358: return(0);
359: }
361: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
362: {
363: DM cda;
364: Vec local_fields;
365: FILE *fp;
366: char fname[PETSC_MAX_PATH_LEN];
367: const char *field_name;
368: PetscMPIInt rank;
369: PetscInt si,sj,nx,ny,i,j,p;
370: PetscInt n_dofs,d;
371: GaussPointCoefficients **_coefficients;
372: PetscErrorCode ierr;
375: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
376: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
377: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
378: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
380: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
381: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
382: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
383: for (d = 0; d < n_dofs; d++) {
384: DMDAGetFieldName(da,d,&field_name);
385: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
386: }
387: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
390: DMGetCoordinateDM(da,&cda);
391: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
393: DMCreateLocalVector(da,&local_fields);
394: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
395: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
396: DMDAVecGetArray(da,local_fields,&_coefficients);
398: for (j = sj; j < sj+ny; j++) {
399: for (i = si; i < si+nx; i++) {
400: PetscScalar coord_x,coord_y;
402: for (p = 0; p < GAUSS_POINTS; p++) {
403: coord_x = _coefficients[j][i].gp_coords[2*p];
404: coord_y = _coefficients[j][i].gp_coords[2*p+1];
406: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
408: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e %1.6e",
409: PetscRealPart(_coefficients[j][i].E[p]),PetscRealPart(_coefficients[j][i].nu[p]),
410: PetscRealPart(_coefficients[j][i].fx[p]),PetscRealPart(_coefficients[j][i].fy[p]));
411: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
412: }
413: }
414: }
415: DMDAVecRestoreArray(da,local_fields,&_coefficients);
416: VecDestroy(&local_fields);
418: PetscFClose(PETSC_COMM_SELF,fp);
419: return(0);
420: }
422: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar E[],PetscScalar nu[])
423: {
424: PetscInt ngp;
425: PetscScalar gp_xi[GAUSS_POINTS][2];
426: PetscScalar gp_weight[GAUSS_POINTS];
427: PetscInt p,i,j,k,l;
428: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
429: PetscScalar J_p;
430: PetscScalar B[3][U_DOFS*NODES_PER_EL];
431: PetscScalar prop_E,prop_nu,factor,constit_D[3][3];
433: /* define quadrature rule */
434: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
436: /* evaluate integral */
437: for (p = 0; p < ngp; p++) {
438: ConstructQ12D_GNi(gp_xi[p],GNi_p);
439: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
441: for (i = 0; i < NODES_PER_EL; i++) {
442: PetscScalar d_dx_i = GNx_p[0][i];
443: PetscScalar d_dy_i = GNx_p[1][i];
445: B[0][2*i] = d_dx_i; B[0][2*i+1] = 0.0;
446: B[1][2*i] = 0.0; B[1][2*i+1] = d_dy_i;
447: B[2][2*i] = d_dy_i; B[2][2*i+1] = d_dx_i;
448: }
450: /* form D for the quadrature point */
451: prop_E = E[p];
452: prop_nu = nu[p];
453: factor = prop_E / ((1.0+prop_nu)*(1.0-2.0*prop_nu));
454: constit_D[0][0] = 1.0-prop_nu; constit_D[0][1] = prop_nu; constit_D[0][2] = 0.0;
455: constit_D[1][0] = prop_nu; constit_D[1][1] = 1.0-prop_nu; constit_D[1][2] = 0.0;
456: constit_D[2][0] = 0.0; constit_D[2][1] = 0.0; constit_D[2][2] = 0.5*(1.0-2.0*prop_nu);
457: for (i = 0; i < 3; i++) {
458: for (j = 0; j < 3; j++) {
459: constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
460: }
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++) {
472: Ke[8*i+j] = Ke[8*i+j] + B[k][i] * constit_D[k][l] * B[l][j];
473: }
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: {
516: /* displacement */
517: /* node 0 */
518: s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0; /* Ux0 */
519: s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1; /* Uy0 */
521: /* node 1 */
522: s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0; /* Ux1 */
523: s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1; /* Uy1 */
525: /* node 2 */
526: s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0; /* Ux2 */
527: s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1; /* Uy2 */
529: /* node 3 */
530: s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0; /* Ux3 */
531: s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1; /* Uy3 */
532: return(0);
533: }
535: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
536: {
538: /* get coords for the element */
539: el_coords[NSD*0+0] = _coords[ej][ei].x; el_coords[NSD*0+1] = _coords[ej][ei].y;
540: el_coords[NSD*1+0] = _coords[ej+1][ei].x; el_coords[NSD*1+1] = _coords[ej+1][ei].y;
541: el_coords[NSD*2+0] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
542: el_coords[NSD*3+0] = _coords[ej][ei+1].x; el_coords[NSD*3+1] = _coords[ej][ei+1].y;
543: return(0);
544: }
546: static PetscErrorCode AssembleA_Elasticity(Mat A,DM elas_da,DM properties_da,Vec properties)
547: {
548: DM cda;
549: Vec coords;
550: DMDACoor2d **_coords;
551: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
552: PetscInt sex,sey,mx,my;
553: PetscInt ei,ej;
554: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
555: PetscScalar el_coords[NODES_PER_EL*NSD];
556: Vec local_properties;
557: GaussPointCoefficients **props;
558: PetscScalar *prop_E,*prop_nu;
559: PetscErrorCode ierr;
562: /* setup for coords */
563: DMGetCoordinateDM(elas_da,&cda);
564: DMGetCoordinatesLocal(elas_da,&coords);
565: DMDAVecGetArray(cda,coords,&_coords);
567: /* setup for coefficients */
568: DMCreateLocalVector(properties_da,&local_properties);
569: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
570: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
571: DMDAVecGetArray(properties_da,local_properties,&props);
573: DMDAGetElementsCorners(elas_da,&sex,&sey,0);
574: DMDAGetElementsSizes(elas_da,&mx,&my,0);
575: for (ej = sey; ej < sey+my; ej++) {
576: for (ei = sex; ei < sex+mx; ei++) {
577: /* get coords for the element */
578: GetElementCoords(_coords,ei,ej,el_coords);
580: /* get coefficients for the element */
581: prop_E = props[ej][ei].E;
582: prop_nu = props[ej][ei].nu;
584: /* initialise element stiffness matrix */
585: PetscMemzero(Ae,sizeof(Ae));
587: /* form element stiffness matrix */
588: FormStressOperatorQ1(Ae,el_coords,prop_E,prop_nu);
590: /* insert element matrix into global matrix */
591: DMDAGetElementEqnums_u(u_eqn,ei,ej);
592: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
593: }
594: }
595: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
596: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
598: DMDAVecRestoreArray(cda,coords,&_coords);
600: DMDAVecRestoreArray(properties_da,local_properties,&props);
601: VecDestroy(&local_properties);
602: return(0);
603: }
606: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F,MatStencil u_eqn[],PetscScalar Fe_u[])
607: {
608: PetscInt n;
611: for (n = 0; n < 4; n++) {
612: 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];
613: 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];
614: }
615: return(0);
616: }
618: static PetscErrorCode AssembleF_Elasticity(Vec F,DM elas_da,DM properties_da,Vec properties)
619: {
620: DM cda;
621: Vec coords;
622: DMDACoor2d **_coords;
623: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
624: PetscInt sex,sey,mx,my;
625: PetscInt ei,ej;
626: PetscScalar Fe[NODES_PER_EL*U_DOFS];
627: PetscScalar el_coords[NODES_PER_EL*NSD];
628: Vec local_properties;
629: GaussPointCoefficients **props;
630: PetscScalar *prop_fx,*prop_fy;
631: Vec local_F;
632: ElasticityDOF **ff;
633: PetscErrorCode ierr;
636: /* setup for coords */
637: DMGetCoordinateDM(elas_da,&cda);
638: DMGetCoordinatesLocal(elas_da,&coords);
639: DMDAVecGetArray(cda,coords,&_coords);
641: /* setup for coefficients */
642: DMGetLocalVector(properties_da,&local_properties);
643: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
644: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
645: DMDAVecGetArray(properties_da,local_properties,&props);
647: /* get acces to the vector */
648: DMGetLocalVector(elas_da,&local_F);
649: VecZeroEntries(local_F);
650: DMDAVecGetArray(elas_da,local_F,&ff);
652: DMDAGetElementsCorners(elas_da,&sex,&sey,0);
653: DMDAGetElementsSizes(elas_da,&mx,&my,0);
654: for (ej = sey; ej < sey+my; ej++) {
655: for (ei = sex; ei < sex+mx; ei++) {
656: /* get coords for the element */
657: GetElementCoords(_coords,ei,ej,el_coords);
659: /* get coefficients for the element */
660: prop_fx = props[ej][ei].fx;
661: prop_fy = props[ej][ei].fy;
663: /* initialise element stiffness matrix */
664: PetscMemzero(Fe,sizeof(Fe));
666: /* form element stiffness matrix */
667: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
669: /* insert element matrix into global matrix */
670: DMDAGetElementEqnums_u(u_eqn,ei,ej);
672: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,Fe);
673: }
674: }
676: DMDAVecRestoreArray(elas_da,local_F,&ff);
677: DMLocalToGlobalBegin(elas_da,local_F,ADD_VALUES,F);
678: DMLocalToGlobalEnd(elas_da,local_F,ADD_VALUES,F);
679: DMRestoreLocalVector(elas_da,&local_F);
681: DMDAVecRestoreArray(cda,coords,&_coords);
683: DMDAVecRestoreArray(properties_da,local_properties,&props);
684: DMRestoreLocalVector(properties_da,&local_properties);
685: return(0);
686: }
688: static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
689: {
690: DM elas_da,da_prop;
691: PetscInt u_dof,dof,stencil_width;
692: Mat A;
693: PetscInt mxl,myl;
694: DM prop_cda,vel_cda;
695: Vec prop_coords,vel_coords;
696: PetscInt si,sj,nx,ny,i,j,p;
697: Vec f,X;
698: PetscInt prop_dof,prop_stencil_width;
699: Vec properties,l_properties;
700: MatNullSpace matnull;
701: PetscReal dx,dy;
702: PetscInt M,N;
703: DMDACoor2d **_prop_coords,**_vel_coords;
704: GaussPointCoefficients **element_props;
705: KSP ksp_E;
706: PetscInt coefficient_structure = 0;
707: PetscInt cpu_x,cpu_y,*lx = NULL,*ly = NULL;
708: PetscBool use_gp_coords = PETSC_FALSE;
709: PetscBool use_nonsymbc = PETSC_FALSE;
710: PetscBool no_view = PETSC_FALSE;
711: PetscBool flg;
712: PetscErrorCode ierr;
715: /* Generate the da for velocity and pressure */
716: /*
717: We use Q1 elements for the temperature.
718: FEM has a 9-point stencil (BOX) or connectivity pattern
719: Num nodes in each direction is mx+1, my+1
720: */
721: u_dof = U_DOFS; /* Vx, Vy - velocities */
722: dof = u_dof;
723: stencil_width = 1;
724: 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);
726: DMSetMatType(elas_da,MATAIJ);
727: DMSetFromOptions(elas_da);
728: DMSetUp(elas_da);
730: DMDASetFieldName(elas_da,0,"Ux");
731: DMDASetFieldName(elas_da,1,"Uy");
733: /* unit box [0,1] x [0,1] */
734: DMDASetUniformCoordinates(elas_da,0.0,1.0,0.0,1.0,0.0,1.0);
736: /* Generate element properties, we will assume all material properties are constant over the element */
737: /* local number of elements */
738: DMDAGetElementsSizes(elas_da,&mxl,&myl,NULL);
740: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
741: DMDAGetInfo(elas_da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
742: DMDAGetElementOwnershipRanges2d(elas_da,&lx,&ly);
744: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
745: prop_stencil_width = 0;
746: 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);
747: DMSetFromOptions(da_prop);
748: DMSetUp(da_prop);
750: PetscFree(lx);
751: PetscFree(ly);
753: /* define centroid positions */
754: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
755: dx = 1.0/((PetscReal)(M));
756: dy = 1.0/((PetscReal)(N));
758: 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);
760: /* define coefficients */
761: PetscOptionsGetInt(NULL,NULL,"-c_str",&coefficient_structure,NULL);
763: DMCreateGlobalVector(da_prop,&properties);
764: DMCreateLocalVector(da_prop,&l_properties);
765: DMDAVecGetArray(da_prop,l_properties,&element_props);
767: DMGetCoordinateDM(da_prop,&prop_cda);
768: DMGetCoordinatesLocal(da_prop,&prop_coords);
769: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
771: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
773: DMGetCoordinateDM(elas_da,&vel_cda);
774: DMGetCoordinatesLocal(elas_da,&vel_coords);
775: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
777: /* interpolate the coordinates */
778: for (j = sj; j < sj+ny; j++) {
779: for (i = si; i < si+nx; i++) {
780: PetscInt ngp;
781: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
782: PetscScalar el_coords[8];
784: GetElementCoords(_vel_coords,i,j,el_coords);
785: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
787: for (p = 0; p < GAUSS_POINTS; p++) {
788: PetscScalar gp_x,gp_y;
789: PetscInt n;
790: PetscScalar xi_p[2],Ni_p[4];
792: xi_p[0] = gp_xi[p][0];
793: xi_p[1] = gp_xi[p][1];
794: ConstructQ12D_Ni(xi_p,Ni_p);
796: gp_x = 0.0;
797: gp_y = 0.0;
798: for (n = 0; n < NODES_PER_EL; n++) {
799: gp_x = gp_x+Ni_p[n]*el_coords[2*n];
800: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
801: }
802: element_props[j][i].gp_coords[2*p] = gp_x;
803: element_props[j][i].gp_coords[2*p+1] = gp_y;
804: }
805: }
806: }
808: /* define the coefficients */
809: PetscOptionsGetBool(NULL,NULL,"-use_gp_coords",&use_gp_coords,&flg);
811: for (j = sj; j < sj+ny; j++) {
812: for (i = si; i < si+nx; i++) {
813: PetscScalar centroid_x = _prop_coords[j][i].x; /* centroids of cell */
814: PetscScalar centroid_y = _prop_coords[j][i].y;
815: PETSC_UNUSED PetscScalar coord_x,coord_y;
817: if (coefficient_structure == 0) { /* isotropic */
818: PetscScalar opts_E,opts_nu;
820: opts_E = 1.0;
821: opts_nu = 0.33;
822: PetscOptionsGetScalar(NULL,NULL,"-iso_E",&opts_E,&flg);
823: PetscOptionsGetScalar(NULL,NULL,"-iso_nu",&opts_nu,&flg);
825: for (p = 0; p < GAUSS_POINTS; p++) {
826: element_props[j][i].E[p] = opts_E;
827: element_props[j][i].nu[p] = opts_nu;
829: element_props[j][i].fx[p] = 0.0;
830: element_props[j][i].fy[p] = 0.0;
831: }
832: } else if (coefficient_structure == 1) { /* step */
833: PetscScalar opts_E0,opts_nu0,opts_xc;
834: PetscScalar opts_E1,opts_nu1;
836: opts_E0 = opts_E1 = 1.0;
837: opts_nu0 = opts_nu1 = 0.333;
838: opts_xc = 0.5;
839: PetscOptionsGetScalar(NULL,NULL,"-step_E0",&opts_E0,&flg);
840: PetscOptionsGetScalar(NULL,NULL,"-step_nu0",&opts_nu0,&flg);
841: PetscOptionsGetScalar(NULL,NULL,"-step_E1",&opts_E1,&flg);
842: PetscOptionsGetScalar(NULL,NULL,"-step_nu1",&opts_nu1,&flg);
843: PetscOptionsGetScalar(NULL,NULL,"-step_xc",&opts_xc,&flg);
845: for (p = 0; p < GAUSS_POINTS; p++) {
846: coord_x = centroid_x;
847: coord_y = centroid_y;
848: if (use_gp_coords) {
849: coord_x = element_props[j][i].gp_coords[2*p];
850: coord_y = element_props[j][i].gp_coords[2*p+1];
851: }
853: element_props[j][i].E[p] = opts_E0;
854: element_props[j][i].nu[p] = opts_nu0;
855: if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
856: element_props[j][i].E[p] = opts_E1;
857: element_props[j][i].nu[p] = opts_nu1;
858: }
860: element_props[j][i].fx[p] = 0.0;
861: element_props[j][i].fy[p] = 0.0;
862: }
863: } else if (coefficient_structure == 2) { /* brick */
864: PetscReal values_E[10];
865: PetscReal values_nu[10];
866: PetscInt nbricks,maxnbricks;
867: PetscInt index,span;
868: PetscInt jj;
870: flg = PETSC_FALSE;
871: maxnbricks = 10;
872: PetscOptionsGetRealArray(NULL,NULL, "-brick_E",values_E,&maxnbricks,&flg);
873: nbricks = maxnbricks;
874: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of E values for each brick");
876: flg = PETSC_FALSE;
877: maxnbricks = 10;
878: PetscOptionsGetRealArray(NULL,NULL, "-brick_nu",values_nu,&maxnbricks,&flg);
879: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of nu values for each brick");
880: if (maxnbricks != nbricks) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply equal numbers of values for E and nu");
882: span = 1;
883: PetscOptionsGetInt(NULL,NULL,"-brick_span",&span,&flg);
885: /* cycle through the indices so that no two material properties are repeated in lines of x or y */
886: jj = (j/span)%nbricks;
887: index = (jj+i/span)%nbricks;
888: /*printf("j=%d: index = %d \n", j,index); */
890: for (p = 0; p < GAUSS_POINTS; p++) {
891: element_props[j][i].E[p] = values_E[index];
892: element_props[j][i].nu[p] = values_nu[index];
893: }
894: } else if (coefficient_structure == 3) { /* sponge */
895: PetscScalar opts_E0,opts_nu0;
896: PetscScalar opts_E1,opts_nu1;
897: PetscInt opts_t,opts_w;
898: PetscInt ii,jj,ci,cj;
900: opts_E0 = opts_E1 = 1.0;
901: opts_nu0 = opts_nu1 = 0.333;
902: PetscOptionsGetScalar(NULL,NULL,"-sponge_E0",&opts_E0,&flg);
903: PetscOptionsGetScalar(NULL,NULL,"-sponge_nu0",&opts_nu0,&flg);
904: PetscOptionsGetScalar(NULL,NULL,"-sponge_E1",&opts_E1,&flg);
905: PetscOptionsGetScalar(NULL,NULL,"-sponge_nu1",&opts_nu1,&flg);
907: opts_t = opts_w = 1;
908: PetscOptionsGetInt(NULL,NULL,"-sponge_t",&opts_t,&flg);
909: PetscOptionsGetInt(NULL,NULL,"-sponge_w",&opts_w,&flg);
911: ii = (i)/(opts_t+opts_w+opts_t);
912: jj = (j)/(opts_t+opts_w+opts_t);
914: ci = i - ii*(opts_t+opts_w+opts_t);
915: cj = j - jj*(opts_t+opts_w+opts_t);
917: for (p = 0; p < GAUSS_POINTS; p++) {
918: element_props[j][i].E[p] = opts_E0;
919: element_props[j][i].nu[p] = opts_nu0;
920: }
921: if ((ci >= opts_t) && (ci < opts_t+opts_w)) {
922: if ((cj >= opts_t) && (cj < opts_t+opts_w)) {
923: for (p = 0; p < GAUSS_POINTS; p++) {
924: element_props[j][i].E[p] = opts_E1;
925: element_props[j][i].nu[p] = opts_nu1;
926: }
927: }
928: }
929: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
930: }
931: }
932: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
934: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
936: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
937: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
938: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
940: PetscOptionsGetBool(NULL,NULL,"-no_view",&no_view,NULL);
941: if (!no_view) {
942: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for elasticity eqn.","properties");
943: DMDACoordViewGnuplot2d(elas_da,"mesh");
944: }
946: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
947: DMCreateMatrix(elas_da,&A);
948: DMGetCoordinates(elas_da,&vel_coords);
949: MatNullSpaceCreateRigidBody(vel_coords,&matnull);
950: MatSetNearNullSpace(A,matnull);
951: MatNullSpaceDestroy(&matnull);
952: MatCreateVecs(A,&f,&X);
954: /* assemble A11 */
955: MatZeroEntries(A);
956: VecZeroEntries(f);
958: AssembleA_Elasticity(A,elas_da,da_prop,properties);
959: /* build force vector */
960: AssembleF_Elasticity(f,elas_da,da_prop,properties);
962: KSPCreate(PETSC_COMM_WORLD,&ksp_E);
963: KSPSetOptionsPrefix(ksp_E,"elas_"); /* elasticity */
965: PetscOptionsGetBool(NULL,NULL,"-use_nonsymbc",&use_nonsymbc,&flg);
966: /* solve */
967: if (!use_nonsymbc) {
968: Mat AA;
969: Vec ff,XX;
970: IS is;
971: VecScatter scat;
973: DMDABCApplySymmetricCompression(elas_da,A,f,&is,&AA,&ff);
974: VecDuplicate(ff,&XX);
976: KSPSetOperators(ksp_E,AA,AA);
977: KSPSetFromOptions(ksp_E);
979: KSPSolve(ksp_E,ff,XX);
981: /* push XX back into X */
982: DMDABCApplyCompression(elas_da,NULL,X);
984: VecScatterCreate(XX,NULL,X,is,&scat);
985: VecScatterBegin(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
986: VecScatterEnd(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
987: VecScatterDestroy(&scat);
989: MatDestroy(&AA);
990: VecDestroy(&ff);
991: VecDestroy(&XX);
992: ISDestroy(&is);
993: } else {
994: DMDABCApplyCompression(elas_da,A,f);
996: KSPSetOperators(ksp_E,A,A);
997: KSPSetFromOptions(ksp_E);
999: KSPSolve(ksp_E,f,X);
1000: }
1002: if (!no_view) {DMDAViewGnuplot2d(elas_da,X,"Displacement solution for elasticity eqn.","X");}
1003: KSPDestroy(&ksp_E);
1005: VecDestroy(&X);
1006: VecDestroy(&f);
1007: MatDestroy(&A);
1009: DMDestroy(&elas_da);
1010: DMDestroy(&da_prop);
1012: VecDestroy(&properties);
1013: VecDestroy(&l_properties);
1014: return(0);
1015: }
1017: int main(int argc,char **args)
1018: {
1020: PetscInt mx,my;
1022: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
1023: mx = my = 10;
1024: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1025: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1026: solve_elasticity_2d(mx,my);
1027: PetscFinalize();
1028: return ierr;
1029: }
1031: /* -------------------------- helpers for boundary conditions -------------------------------- */
1033: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1034: {
1035: DM cda;
1036: Vec coords;
1037: PetscInt si,sj,nx,ny,i,j;
1038: PetscInt M,N;
1039: DMDACoor2d **_coords;
1040: const PetscInt *g_idx;
1041: PetscInt *bc_global_ids;
1042: PetscScalar *bc_vals;
1043: PetscInt nbcs;
1044: PetscInt n_dofs;
1045: PetscErrorCode ierr;
1046: ISLocalToGlobalMapping ltogm;
1049: /* enforce bc's */
1050: DMGetLocalToGlobalMapping(da,<ogm);
1051: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1053: DMGetCoordinateDM(da,&cda);
1054: DMGetCoordinatesLocal(da,&coords);
1055: DMDAVecGetArray(cda,coords,&_coords);
1056: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1057: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1059: /* --- */
1061: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1062: PetscMalloc1(ny*n_dofs,&bc_vals);
1064: /* init the entries to -1 so VecSetValues will ignore them */
1065: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1067: i = nx-1;
1068: for (j = 0; j < ny; j++) {
1069: PetscInt local_id;
1070: PETSC_UNUSED PetscScalar coordx,coordy;
1072: local_id = i+j*nx;
1074: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1076: coordx = _coords[j+sj][i+si].x;
1077: coordy = _coords[j+sj][i+si].y;
1079: bc_vals[j] = bc_val;
1080: }
1081: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1082: nbcs = 0;
1083: if ((si+nx) == (M)) nbcs = ny;
1085: if (b) {
1086: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1087: VecAssemblyBegin(b);
1088: VecAssemblyEnd(b);
1089: }
1090: if (A) {
1091: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1092: }
1094: PetscFree(bc_vals);
1095: PetscFree(bc_global_ids);
1097: DMDAVecRestoreArray(cda,coords,&_coords);
1098: return(0);
1099: }
1101: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1102: {
1103: DM cda;
1104: Vec coords;
1105: PetscInt si,sj,nx,ny,i,j;
1106: PetscInt M,N;
1107: DMDACoor2d **_coords;
1108: const PetscInt *g_idx;
1109: PetscInt *bc_global_ids;
1110: PetscScalar *bc_vals;
1111: PetscInt nbcs;
1112: PetscInt n_dofs;
1113: PetscErrorCode ierr;
1114: ISLocalToGlobalMapping ltogm;
1117: /* enforce bc's */
1118: DMGetLocalToGlobalMapping(da,<ogm);
1119: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1121: DMGetCoordinateDM(da,&cda);
1122: DMGetCoordinatesLocal(da,&coords);
1123: DMDAVecGetArray(cda,coords,&_coords);
1124: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1125: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1127: /* --- */
1129: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1130: PetscMalloc1(ny*n_dofs,&bc_vals);
1132: /* init the entries to -1 so VecSetValues will ignore them */
1133: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1135: i = 0;
1136: for (j = 0; j < ny; j++) {
1137: PetscInt local_id;
1138: PETSC_UNUSED PetscScalar coordx,coordy;
1140: local_id = i+j*nx;
1142: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1144: coordx = _coords[j+sj][i+si].x;
1145: coordy = _coords[j+sj][i+si].y;
1147: bc_vals[j] = bc_val;
1148: }
1149: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1150: nbcs = 0;
1151: if (si == 0) nbcs = ny;
1153: if (b) {
1154: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1155: VecAssemblyBegin(b);
1156: VecAssemblyEnd(b);
1157: }
1158: if (A) {
1159: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1160: }
1162: PetscFree(bc_vals);
1163: PetscFree(bc_global_ids);
1165: DMDAVecRestoreArray(cda,coords,&_coords);
1166: return(0);
1167: }
1169: static PetscErrorCode DMDABCApplyCompression(DM elas_da,Mat A,Vec f)
1170: {
1174: BCApply_EAST(elas_da,0,-1.0,A,f);
1175: BCApply_EAST(elas_da,1, 0.0,A,f);
1176: BCApply_WEST(elas_da,0,1.0,A,f);
1177: BCApply_WEST(elas_da,1,0.0,A,f);
1178: return(0);
1179: }
1181: static PetscErrorCode Orthogonalize(PetscInt n,Vec *vecs)
1182: {
1183: PetscInt i,j;
1184: PetscScalar dot;
1188: for (i=0; i<n; i++) {
1189: VecNormalize(vecs[i],NULL);
1190: for (j=i+1; j<n; j++) {
1191: VecDot(vecs[i],vecs[j],&dot);
1192: VecAXPY(vecs[j],-dot,vecs[i]);
1193: }
1194: }
1195: return(0);
1196: }
1198: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff)
1199: {
1201: PetscInt start,end,m;
1202: PetscInt *unconstrained;
1203: PetscInt cnt,i;
1204: Vec x;
1205: PetscScalar *_x;
1206: IS is;
1207: VecScatter scat;
1210: /* push bc's into f and A */
1211: VecDuplicate(f,&x);
1212: BCApply_EAST(elas_da,0,-1.0,A,x);
1213: BCApply_EAST(elas_da,1, 0.0,A,x);
1214: BCApply_WEST(elas_da,0,1.0,A,x);
1215: BCApply_WEST(elas_da,1,0.0,A,x);
1217: /* define which dofs are not constrained */
1218: VecGetLocalSize(x,&m);
1219: PetscMalloc1(m,&unconstrained);
1220: VecGetOwnershipRange(x,&start,&end);
1221: VecGetArray(x,&_x);
1222: cnt = 0;
1223: for (i = 0; i < m; i+=2) {
1224: PetscReal val1,val2;
1226: val1 = PetscRealPart(_x[i]);
1227: val2 = PetscRealPart(_x[i+1]);
1228: if (PetscAbs(val1) < 0.1 && PetscAbs(val2) < 0.1) {
1229: unconstrained[cnt] = start + i;
1230: cnt++;
1231: unconstrained[cnt] = start + i + 1;
1232: cnt++;
1233: }
1234: }
1235: VecRestoreArray(x,&_x);
1237: ISCreateGeneral(PETSC_COMM_WORLD,cnt,unconstrained,PETSC_COPY_VALUES,&is);
1238: PetscFree(unconstrained);
1239: ISSetBlockSize(is,2);
1241: /* define correction for dirichlet in the rhs */
1242: MatMult(A,x,f);
1243: VecScale(f,-1.0);
1245: /* get new matrix */
1246: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,AA);
1247: /* get new vector */
1248: MatCreateVecs(*AA,NULL,ff);
1250: VecScatterCreate(f,is,*ff,NULL,&scat);
1251: VecScatterBegin(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1252: VecScatterEnd(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1254: { /* Constrain near-null space */
1255: PetscInt nvecs;
1256: const Vec *vecs;
1257: Vec *uvecs;
1258: PetscBool has_const;
1259: MatNullSpace mnull,unull;
1261: MatGetNearNullSpace(A,&mnull);
1262: MatNullSpaceGetVecs(mnull,&has_const,&nvecs,&vecs);
1263: VecDuplicateVecs(*ff,nvecs,&uvecs);
1264: for (i=0; i<nvecs; i++) {
1265: VecScatterBegin(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1266: VecScatterEnd(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1267: }
1268: Orthogonalize(nvecs,uvecs);
1269: MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,nvecs,uvecs,&unull);
1270: MatSetNearNullSpace(*AA,unull);
1271: MatNullSpaceDestroy(&unull);
1272: VecDestroyVecs(nvecs,&uvecs);
1273: }
1275: VecScatterDestroy(&scat);
1277: *dofs = is;
1278: VecDestroy(&x);
1279: return(0);
1280: }
1283: /*TEST
1285: build:
1286: requires: !complex !single
1288: test:
1289: 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
1290: output_file: output/ex49_1.out
1292: test:
1293: suffix: 2
1294: nsize: 4
1295: 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
1297: test:
1298: suffix: 3
1299: nsize: 4
1300: 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
1302: test:
1303: suffix: 4
1304: nsize: 4
1305: 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
1306: requires: ml
1308: test:
1309: suffix: 5
1310: nsize: 3
1311: 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
1313: test:
1314: suffix: 6
1315: nsize: 4
1316: 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
1318: test:
1319: suffix: 7
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 ksp -elas_sub_ksp_ksp_type cg -elas_sub_ksp_ksp_max_it 15
1323: test:
1324: suffix: 8
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 pipefgmres -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: hypre_nullspace
1330: requires: hypre
1331: 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_ksp_view
1333: test:
1334: nsize: 4
1335: suffix: bddc
1336: 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
1338: test:
1339: nsize: 4
1340: suffix: bddc_unsym
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 -use_nonsymbc -elas_pc_bddc_symmetric 0
1343: test:
1344: nsize: 4
1345: suffix: bddc_unsym_deluxe
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 -elas_pc_bddc_use_deluxe_scaling -elas_sub_schurs_symmetric 0
1348: test:
1349: nsize: 4
1350: suffix: fetidp_unsym_deluxe
1351: 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
1353: test:
1354: nsize: 4
1355: suffix: bddc_layerjump
1356: 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
1358: test:
1359: nsize: 4
1360: suffix: bddc_subdomainjump
1361: 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
1363: test:
1364: nsize: 9
1365: suffix: bddc_subdomainjump_deluxe
1366: 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
1367: TEST*/