Actual source code: ex49.c
petsc-3.5.4 2015-05-23
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 moduls) and nu (Poisson ratio) may vary as a function of space. \n\
3: The model utilisse boundary conditions which produce compression in the x direction. \n\
4: Options: \n\
5: -mx : number elements in x-direciton \n\
6: -my : number elements in y-direciton \n\
7: -c_str : indicates the structure of the coefficients to use. \n\
8: -c_str 0 => Setup for an isotropic material with constant coefficients. \n\
9: Parameters: \n\
10: -iso_E : Youngs modulus \n\
11: -iso_nu : Poisson ratio \n\
12: -c_str 1 => Setup for a step function in the material properties in x. \n\
13: Parameters: \n\
14: -step_E0 : Youngs modulus to the left of the step \n\
15: -step_nu0 : Poisson ratio to the left of the step \n\
16: -step_E1 : Youngs modulus to the right of the step \n\
17: -step_n1 : Poisson ratio to the right of the step \n\
18: -step_xc : x coordinate of the step \n\
19: -c_str 2 => Setup for a checkerboard material with alternating properties. \n\
20: Repeats the following pattern throughout the domain. For example with 4 materials specified, we would heve \n\
21: -------------------------\n\
22: | D | A | B | C |\n\
23: ------|-----|-----|------\n\
24: | C | D | A | B |\n\
25: ------|-----|-----|------\n\
26: | B | C | D | A |\n\
27: ------|-----|-----|------\n\
28: | A | B | C | D |\n\
29: -------------------------\n\
30: \n\
31: Parameters: \n\
32: -brick_E : a comma seperated list of Young's modulii \n\
33: -brick_nu : a comma seperated list of Poisson ratio's \n\
34: -brick_span : the number of elements in x and y each brick will span \n\
35: -c_str 3 => Setup for a sponge-like material with alternating properties. \n\
36: Repeats the following pattern throughout the domain \n\
37: -----------------------------\n\
38: | [background] |\n\
39: | E0,nu0 |\n\
40: | ----------------- |\n\
41: | | [inclusion] | |\n\
42: | | E1,nu1 | |\n\
43: | | | |\n\
44: | | <---- w ----> | |\n\
45: | | | |\n\
46: | | | |\n\
47: | ----------------- |\n\
48: | |\n\
49: | |\n\
50: -----------------------------\n\
51: <-------- t + w + t ------->\n\
52: \n\
53: Parameters: \n\
54: -sponge_E0 : Youngs moduls of the surrounding material \n\
55: -sponge_E1 : Youngs moduls of the inclusio \n\
56: -sponge_nu0 : Poisson ratio of the surrounding material \n\
57: -sponge_nu1 : Poisson ratio of the inclusio \n\
58: -sponge_t : the number of elements defining the border around each inclusion \n\
59: -sponge_w : the number of elements in x and y each inclusion will span\n\
60: -use_gp_coords : Evaluate the Youngs modulus, Poisson ratio and the body force at the global coordinates of the quadrature points.\n\
61: By default, E, nu and the body force are evaulated at the element center and applied as a constant over the entire element.\n\
62: -use_nonsymbc : Option to use non-symmetric boundary condition imposition. This choice will use less memory.";
64: /* Contributed by Dave May */
66: #include <petscksp.h>
67: #include <petscdm.h>
68: #include <petscdmda.h>
70: static PetscErrorCode DMDABCApplyCompression(DM,Mat,Vec);
71: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff);
74: #define NSD 2 /* number of spatial dimensions */
75: #define NODES_PER_EL 4 /* nodes per element */
76: #define U_DOFS 2 /* degrees of freedom per displacement node */
77: #define GAUSS_POINTS 4
79: /* cell based evaluation */
80: typedef struct {
81: PetscScalar E,nu,fx,fy;
82: } Coefficients;
84: /* Gauss point based evaluation 8+4+4+4 = 20 */
85: typedef struct {
86: PetscScalar gp_coords[2*GAUSS_POINTS];
87: PetscScalar E[GAUSS_POINTS];
88: PetscScalar nu[GAUSS_POINTS];
89: PetscScalar fx[GAUSS_POINTS];
90: PetscScalar fy[GAUSS_POINTS];
91: } GaussPointCoefficients;
93: typedef struct {
94: PetscScalar ux_dof;
95: PetscScalar uy_dof;
96: } ElasticityDOF;
99: /*
101: D = E/((1+nu)(1-2nu)) * [ 1-nu nu 0 ]
102: [ nu 1-nu 0 ]
103: [ 0 0 0.5*(1-2nu) ]
105: B = [ d_dx 0 ]
106: [ 0 d_dy ]
107: [ d_dy d_dx ]
109: */
111: /* FEM routines */
112: /*
113: Element: Local basis function ordering
114: 1-----2
115: | |
116: | |
117: 0-----3
118: */
119: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
120: {
121: PetscScalar xi = _xi[0];
122: PetscScalar eta = _xi[1];
124: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
125: Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
126: Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
127: Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
128: }
130: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
131: {
132: PetscScalar xi = _xi[0];
133: PetscScalar eta = _xi[1];
135: GNi[0][0] = -0.25*(1.0-eta);
136: GNi[0][1] = -0.25*(1.0+eta);
137: GNi[0][2] = 0.25*(1.0+eta);
138: GNi[0][3] = 0.25*(1.0-eta);
140: GNi[1][0] = -0.25*(1.0-xi);
141: GNi[1][1] = 0.25*(1.0-xi);
142: GNi[1][2] = 0.25*(1.0+xi);
143: GNi[1][3] = -0.25*(1.0+xi);
144: }
146: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
147: {
148: PetscScalar J00,J01,J10,J11,J;
149: PetscScalar iJ00,iJ01,iJ10,iJ11;
150: PetscInt i;
152: J00 = J01 = J10 = J11 = 0.0;
153: for (i = 0; i < NODES_PER_EL; i++) {
154: PetscScalar cx = coords[2*i+0];
155: PetscScalar cy = coords[2*i+1];
157: J00 = J00+GNi[0][i]*cx; /* J_xx = dx/dxi */
158: J01 = J01+GNi[0][i]*cy; /* J_xy = dy/dxi */
159: J10 = J10+GNi[1][i]*cx; /* J_yx = dx/deta */
160: J11 = J11+GNi[1][i]*cy; /* J_yy = dy/deta */
161: }
162: J = (J00*J11)-(J01*J10);
164: iJ00 = J11/J;
165: iJ01 = -J01/J;
166: iJ10 = -J10/J;
167: iJ11 = J00/J;
170: for (i = 0; i < NODES_PER_EL; i++) {
171: GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
172: GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
173: }
175: if (det_J != NULL) *det_J = J;
176: }
178: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
179: {
180: *ngp = 4;
181: gp_xi[0][0] = -0.57735026919;gp_xi[0][1] = -0.57735026919;
182: gp_xi[1][0] = -0.57735026919;gp_xi[1][1] = 0.57735026919;
183: gp_xi[2][0] = 0.57735026919;gp_xi[2][1] = 0.57735026919;
184: gp_xi[3][0] = 0.57735026919;gp_xi[3][1] = -0.57735026919;
185: gp_weight[0] = 1.0;
186: gp_weight[1] = 1.0;
187: gp_weight[2] = 1.0;
188: gp_weight[3] = 1.0;
189: }
192: /* procs to the left claim the ghost node as their element */
195: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
196: {
198: PetscInt m,n,p,M,N,P;
199: PetscInt sx,sy,sz;
202: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
203: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
205: if (mxl != NULL) {
206: *mxl = m;
207: if ((sx+m) == M) *mxl = m-1; /* last proc */
208: }
209: if (myl != NULL) {
210: *myl = n;
211: if ((sy+n) == N) *myl = n-1; /* last proc */
212: }
213: if (mzl != NULL) {
214: *mzl = p;
215: if ((sz+p) == P) *mzl = p-1; /* last proc */
216: }
217: return(0);
218: }
222: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
223: {
225: PetscInt si,sj,sk;
228: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
230: if (sx) {
231: *sx = si;
232: if (si != 0) *sx = si+1;
233: }
234: if (sy) {
235: *sy = sj;
236: if (sj != 0) *sy = sj+1;
237: }
239: if (sk) {
240: *sz = sk;
241: if (sk != 0) *sz = sk+1;
242: }
244: DMDAGetLocalElementSize(da,mx,my,mz);
245: return(0);
246: }
250: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
251: {
253: PetscMPIInt rank;
254: PetscInt proc_I,proc_J;
255: PetscInt cpu_x,cpu_y;
256: PetscInt local_mx,local_my;
257: Vec vlx,vly;
258: PetscInt *LX,*LY,i;
259: PetscScalar *_a;
260: Vec V_SEQ;
261: VecScatter ctx;
264: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
266: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
268: proc_J = rank/cpu_x;
269: proc_I = rank-cpu_x*proc_J;
271: PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);
272: PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);
274: DMDAGetLocalElementSize(da,&local_mx,&local_my,NULL);
275: VecCreate(PETSC_COMM_WORLD,&vlx);
276: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
277: VecSetFromOptions(vlx);
279: VecCreate(PETSC_COMM_WORLD,&vly);
280: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
281: VecSetFromOptions(vly);
283: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
284: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
285: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
286: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
290: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
291: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
292: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
293: VecGetArray(V_SEQ,&_a);
294: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
295: VecRestoreArray(V_SEQ,&_a);
296: VecScatterDestroy(&ctx);
297: VecDestroy(&V_SEQ);
299: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
300: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
301: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
302: VecGetArray(V_SEQ,&_a);
303: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
304: VecRestoreArray(V_SEQ,&_a);
305: VecScatterDestroy(&ctx);
306: VecDestroy(&V_SEQ);
308: *_lx = LX;
309: *_ly = LY;
311: VecDestroy(&vlx);
312: VecDestroy(&vly);
313: return(0);
314: }
318: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
319: {
320: DM cda;
321: Vec coords;
322: DMDACoor2d **_coords;
323: PetscInt si,sj,nx,ny,i,j;
324: FILE *fp;
325: char fname[PETSC_MAX_PATH_LEN];
326: PetscMPIInt rank;
330: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
331: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
332: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
333: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
334: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
336: DMGetCoordinateDM(da,&cda);
337: DMGetCoordinatesLocal(da,&coords);
338: DMDAVecGetArray(cda,coords,&_coords);
339: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
340: for (j = sj; j < sj+ny-1; j++) {
341: for (i = si; i < si+nx-1; i++) {
342: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
343: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i].x),PetscRealPart(_coords[j+1][i].y));
344: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i+1].x),PetscRealPart(_coords[j+1][i+1].y));
345: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i+1].x),PetscRealPart(_coords[j][i+1].y));
346: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
347: }
348: }
349: DMDAVecRestoreArray(cda,coords,&_coords);
351: PetscFClose(PETSC_COMM_SELF,fp);
352: return(0);
353: }
357: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
358: {
359: DM cda;
360: Vec coords,local_fields;
361: DMDACoor2d **_coords;
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;
367: PetscInt n_dofs,d;
368: PetscScalar *_fields;
372: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
373: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
374: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
375: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
377: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
378: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
379: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
380: for (d = 0; d < n_dofs; d++) {
381: DMDAGetFieldName(da,d,&field_name);
382: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
383: }
384: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
387: DMGetCoordinateDM(da,&cda);
388: DMGetCoordinatesLocal(da,&coords);
389: DMDAVecGetArray(cda,coords,&_coords);
390: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
392: DMCreateLocalVector(da,&local_fields);
393: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
394: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
395: VecGetArray(local_fields,&_fields);
397: for (j = sj; j < sj+ny; j++) {
398: for (i = si; i < si+nx; i++) {
399: PetscScalar coord_x,coord_y;
400: PetscScalar field_d;
402: coord_x = _coords[j][i].x;
403: coord_y = _coords[j][i].y;
405: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
406: for (d = 0; d < n_dofs; d++) {
407: field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
408: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",PetscRealPart(field_d));
409: }
410: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
411: }
412: }
413: VecRestoreArray(local_fields,&_fields);
414: VecDestroy(&local_fields);
416: DMDAVecRestoreArray(cda,coords,&_coords);
418: PetscFClose(PETSC_COMM_SELF,fp);
419: return(0);
420: }
424: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
425: {
426: DM cda;
427: Vec local_fields;
428: FILE *fp;
429: char fname[PETSC_MAX_PATH_LEN];
430: const char *field_name;
431: PetscMPIInt rank;
432: PetscInt si,sj,nx,ny,i,j,p;
433: PetscInt n_dofs,d;
434: GaussPointCoefficients **_coefficients;
435: PetscErrorCode ierr;
438: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
439: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
440: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
441: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
443: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
444: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
445: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
446: for (d = 0; d < n_dofs; d++) {
447: DMDAGetFieldName(da,d,&field_name);
448: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
449: }
450: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
453: DMGetCoordinateDM(da,&cda);
454: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
456: DMCreateLocalVector(da,&local_fields);
457: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
458: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
459: DMDAVecGetArray(da,local_fields,&_coefficients);
461: for (j = sj; j < sj+ny; j++) {
462: for (i = si; i < si+nx; i++) {
463: PetscScalar coord_x,coord_y;
465: for (p = 0; p < GAUSS_POINTS; p++) {
466: coord_x = _coefficients[j][i].gp_coords[2*p];
467: coord_y = _coefficients[j][i].gp_coords[2*p+1];
469: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
471: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e %1.6e",
472: PetscRealPart(_coefficients[j][i].E[p]),PetscRealPart(_coefficients[j][i].nu[p]),
473: PetscRealPart(_coefficients[j][i].fx[p]),PetscRealPart(_coefficients[j][i].fy[p]));
474: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
475: }
476: }
477: }
478: DMDAVecRestoreArray(da,local_fields,&_coefficients);
479: VecDestroy(&local_fields);
481: PetscFClose(PETSC_COMM_SELF,fp);
482: return(0);
483: }
485: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar E[],PetscScalar nu[])
486: {
487: PetscInt ngp;
488: PetscScalar gp_xi[GAUSS_POINTS][2];
489: PetscScalar gp_weight[GAUSS_POINTS];
490: PetscInt p,i,j,k,l;
491: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
492: PetscScalar J_p;
493: PetscScalar B[3][U_DOFS*NODES_PER_EL];
494: PetscScalar prop_E,prop_nu,factor,constit_D[3][3];
496: /* define quadrature rule */
497: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
499: /* evaluate integral */
500: for (p = 0; p < ngp; p++) {
501: ConstructQ12D_GNi(gp_xi[p],GNi_p);
502: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
504: for (i = 0; i < NODES_PER_EL; i++) {
505: PetscScalar d_dx_i = GNx_p[0][i];
506: PetscScalar d_dy_i = GNx_p[1][i];
508: B[0][2*i] = d_dx_i; B[0][2*i+1] = 0.0;
509: B[1][2*i] = 0.0; B[1][2*i+1] = d_dy_i;
510: B[2][2*i] = d_dy_i; B[2][2*i+1] = d_dx_i;
511: }
513: /* form D for the quadrature point */
514: prop_E = E[p];
515: prop_nu = nu[p];
516: factor = prop_E / ((1.0+prop_nu)*(1.0-2.0*prop_nu));
517: constit_D[0][0] = 1.0-prop_nu; constit_D[0][1] = prop_nu; constit_D[0][2] = 0.0;
518: constit_D[1][0] = prop_nu; constit_D[1][1] = 1.0-prop_nu; constit_D[1][2] = 0.0;
519: constit_D[2][0] = 0.0; constit_D[2][1] = 0.0; constit_D[2][2] = 0.5*(1.0-2.0*prop_nu);
520: for (i = 0; i < 3; i++) {
521: for (j = 0; j < 3; j++) {
522: constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
523: }
524: }
526: /* form Bt tildeD B */
527: /*
528: Ke_ij = Bt_ik . D_kl . B_lj
529: = B_ki . D_kl . B_lj
530: */
531: for (i = 0; i < 8; i++) {
532: for (j = 0; j < 8; j++) {
533: for (k = 0; k < 3; k++) {
534: for (l = 0; l < 3; l++) {
535: Ke[8*i+j] = Ke[8*i+j] + B[k][i] * constit_D[k][l] * B[l][j];
536: }
537: }
538: }
539: }
541: } /* end quadrature */
542: }
544: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
545: {
546: PetscInt ngp;
547: PetscScalar gp_xi[GAUSS_POINTS][2];
548: PetscScalar gp_weight[GAUSS_POINTS];
549: PetscInt p,i;
550: PetscScalar Ni_p[NODES_PER_EL];
551: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
552: PetscScalar J_p,fac;
554: /* define quadrature rule */
555: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
557: /* evaluate integral */
558: for (p = 0; p < ngp; p++) {
559: ConstructQ12D_Ni(gp_xi[p],Ni_p);
560: ConstructQ12D_GNi(gp_xi[p],GNi_p);
561: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
562: fac = gp_weight[p]*J_p;
564: for (i = 0; i < NODES_PER_EL; i++) {
565: Fe[NSD*i] += fac*Ni_p[i]*fx[p];
566: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
567: }
568: }
569: }
571: /*
572: i,j are the element indices
573: The unknown is a vector quantity.
574: The s[].c is used to indicate the degree of freedom.
575: */
578: static PetscErrorCode DMDAGetElementEqnums_u(MatStencil s_u[],PetscInt i,PetscInt j)
579: {
581: /* displacement */
582: /* node 0 */
583: s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0; /* Ux0 */
584: s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1; /* Uy0 */
586: /* node 1 */
587: s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0; /* Ux1 */
588: s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1; /* Uy1 */
590: /* node 2 */
591: s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0; /* Ux2 */
592: s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1; /* Uy2 */
594: /* node 3 */
595: s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0; /* Ux3 */
596: s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1; /* Uy3 */
597: return(0);
598: }
602: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
603: {
605: /* get coords for the element */
606: el_coords[NSD*0+0] = _coords[ej][ei].x; el_coords[NSD*0+1] = _coords[ej][ei].y;
607: el_coords[NSD*1+0] = _coords[ej+1][ei].x; el_coords[NSD*1+1] = _coords[ej+1][ei].y;
608: el_coords[NSD*2+0] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
609: el_coords[NSD*3+0] = _coords[ej][ei+1].x; el_coords[NSD*3+1] = _coords[ej][ei+1].y;
610: return(0);
611: }
615: static PetscErrorCode AssembleA_Elasticity(Mat A,DM elas_da,DM properties_da,Vec properties)
616: {
617: DM cda;
618: Vec coords;
619: DMDACoor2d **_coords;
620: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
621: PetscInt sex,sey,mx,my;
622: PetscInt ei,ej;
623: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
624: PetscScalar el_coords[NODES_PER_EL*NSD];
625: Vec local_properties;
626: GaussPointCoefficients **props;
627: PetscScalar *prop_E,*prop_nu;
628: PetscErrorCode ierr;
631: /* setup for coords */
632: DMGetCoordinateDM(elas_da,&cda);
633: DMGetCoordinatesLocal(elas_da,&coords);
634: DMDAVecGetArray(cda,coords,&_coords);
636: /* setup for coefficients */
637: DMCreateLocalVector(properties_da,&local_properties);
638: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
639: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
640: DMDAVecGetArray(properties_da,local_properties,&props);
642: DMDAGetElementCorners(elas_da,&sex,&sey,0,&mx,&my,0);
643: for (ej = sey; ej < sey+my; ej++) {
644: for (ei = sex; ei < sex+mx; ei++) {
645: /* get coords for the element */
646: GetElementCoords(_coords,ei,ej,el_coords);
648: /* get coefficients for the element */
649: prop_E = props[ej][ei].E;
650: prop_nu = props[ej][ei].nu;
652: /* initialise element stiffness matrix */
653: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
655: /* form element stiffness matrix */
656: FormStressOperatorQ1(Ae,el_coords,prop_E,prop_nu);
658: /* insert element matrix into global matrix */
659: DMDAGetElementEqnums_u(u_eqn,ei,ej);
660: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
661: }
662: }
663: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
664: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
666: DMDAVecRestoreArray(cda,coords,&_coords);
668: DMDAVecRestoreArray(properties_da,local_properties,&props);
669: VecDestroy(&local_properties);
670: return(0);
671: }
676: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F,MatStencil u_eqn[],PetscScalar Fe_u[])
677: {
678: PetscInt n;
681: for (n = 0; n < 4; n++) {
682: 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];
683: 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];
684: }
685: return(0);
686: }
690: static PetscErrorCode AssembleF_Elasticity(Vec F,DM elas_da,DM properties_da,Vec properties)
691: {
692: DM cda;
693: Vec coords;
694: DMDACoor2d **_coords;
695: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
696: PetscInt sex,sey,mx,my;
697: PetscInt ei,ej;
698: PetscScalar Fe[NODES_PER_EL*U_DOFS];
699: PetscScalar el_coords[NODES_PER_EL*NSD];
700: Vec local_properties;
701: GaussPointCoefficients **props;
702: PetscScalar *prop_fx,*prop_fy;
703: Vec local_F;
704: ElasticityDOF **ff;
705: PetscErrorCode ierr;
708: /* setup for coords */
709: DMGetCoordinateDM(elas_da,&cda);
710: DMGetCoordinatesLocal(elas_da,&coords);
711: DMDAVecGetArray(cda,coords,&_coords);
713: /* setup for coefficients */
714: DMGetLocalVector(properties_da,&local_properties);
715: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
716: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
717: DMDAVecGetArray(properties_da,local_properties,&props);
719: /* get acces to the vector */
720: DMGetLocalVector(elas_da,&local_F);
721: VecZeroEntries(local_F);
722: DMDAVecGetArray(elas_da,local_F,&ff);
725: DMDAGetElementCorners(elas_da,&sex,&sey,0,&mx,&my,0);
726: for (ej = sey; ej < sey+my; ej++) {
727: for (ei = sex; ei < sex+mx; ei++) {
728: /* get coords for the element */
729: GetElementCoords(_coords,ei,ej,el_coords);
731: /* get coefficients for the element */
732: prop_fx = props[ej][ei].fx;
733: prop_fy = props[ej][ei].fy;
735: /* initialise element stiffness matrix */
736: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
738: /* form element stiffness matrix */
739: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
741: /* insert element matrix into global matrix */
742: DMDAGetElementEqnums_u(u_eqn,ei,ej);
744: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,Fe);
745: }
746: }
748: DMDAVecRestoreArray(elas_da,local_F,&ff);
749: DMLocalToGlobalBegin(elas_da,local_F,ADD_VALUES,F);
750: DMLocalToGlobalEnd(elas_da,local_F,ADD_VALUES,F);
751: DMRestoreLocalVector(elas_da,&local_F);
753: DMDAVecRestoreArray(cda,coords,&_coords);
755: DMDAVecRestoreArray(properties_da,local_properties,&props);
756: DMRestoreLocalVector(properties_da,&local_properties);
757: return(0);
758: }
762: static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
763: {
764: DM elas_da,da_prop;
765: PetscInt u_dof,dof,stencil_width;
766: Mat A;
767: PetscInt mxl,myl;
768: DM prop_cda,vel_cda;
769: Vec prop_coords,vel_coords;
770: PetscInt si,sj,nx,ny,i,j,p;
771: Vec f,X;
772: PetscInt prop_dof,prop_stencil_width;
773: Vec properties,l_properties;
774: MatNullSpace matnull;
775: PetscReal dx,dy;
776: PetscInt M,N;
777: DMDACoor2d **_prop_coords,**_vel_coords;
778: GaussPointCoefficients **element_props;
779: KSP ksp_E;
780: PetscInt coefficient_structure = 0;
781: PetscInt cpu_x,cpu_y,*lx = NULL,*ly = NULL;
782: PetscBool use_gp_coords = PETSC_FALSE;
783: PetscBool use_nonsymbc = PETSC_FALSE;
784: PetscBool no_view = PETSC_FALSE;
785: PetscBool flg;
786: PetscErrorCode ierr;
789: /* Generate the da for velocity and pressure */
790: /*
791: We use Q1 elements for the temperature.
792: FEM has a 9-point stencil (BOX) or connectivity pattern
793: Num nodes in each direction is mx+1, my+1
794: */
795: u_dof = U_DOFS; /* Vx, Vy - velocities */
796: dof = u_dof;
797: stencil_width = 1;
798: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
799: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&elas_da);
800: DMDASetFieldName(elas_da,0,"Ux");
801: DMDASetFieldName(elas_da,1,"Uy");
803: /* unit box [0,1] x [0,1] */
804: DMDASetUniformCoordinates(elas_da,0.0,1.0,0.0,1.0,0.0,1.0);
807: /* Generate element properties, we will assume all material properties are constant over the element */
808: /* local number of elements */
809: DMDAGetLocalElementSize(elas_da,&mxl,&myl,NULL);
811: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
812: DMDAGetInfo(elas_da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
813: DMDAGetElementOwnershipRanges2d(elas_da,&lx,&ly);
815: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
816: prop_stencil_width = 0;
817: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
818: mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
819: PetscFree(lx);
820: PetscFree(ly);
822: /* define centroid positions */
823: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
824: dx = 1.0/((PetscReal)(M));
825: dy = 1.0/((PetscReal)(N));
827: 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);
829: /* define coefficients */
830: PetscOptionsGetInt(NULL,"-c_str",&coefficient_structure,NULL);
832: DMCreateGlobalVector(da_prop,&properties);
833: DMCreateLocalVector(da_prop,&l_properties);
834: DMDAVecGetArray(da_prop,l_properties,&element_props);
836: DMGetCoordinateDM(da_prop,&prop_cda);
837: DMGetCoordinatesLocal(da_prop,&prop_coords);
838: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
840: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
842: DMGetCoordinateDM(elas_da,&vel_cda);
843: DMGetCoordinatesLocal(elas_da,&vel_coords);
844: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
847: /* interpolate the coordinates */
848: for (j = sj; j < sj+ny; j++) {
849: for (i = si; i < si+nx; i++) {
850: PetscInt ngp;
851: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
852: PetscScalar el_coords[8];
854: GetElementCoords(_vel_coords,i,j,el_coords);
855: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
857: for (p = 0; p < GAUSS_POINTS; p++) {
858: PetscScalar gp_x,gp_y;
859: PetscInt n;
860: PetscScalar xi_p[2],Ni_p[4];
862: xi_p[0] = gp_xi[p][0];
863: xi_p[1] = gp_xi[p][1];
864: ConstructQ12D_Ni(xi_p,Ni_p);
866: gp_x = 0.0;
867: gp_y = 0.0;
868: for (n = 0; n < NODES_PER_EL; n++) {
869: gp_x = gp_x+Ni_p[n]*el_coords[2*n];
870: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
871: }
872: element_props[j][i].gp_coords[2*p] = gp_x;
873: element_props[j][i].gp_coords[2*p+1] = gp_y;
874: }
875: }
876: }
878: /* define the coefficients */
879: PetscOptionsGetBool(NULL,"-use_gp_coords",&use_gp_coords,&flg);
881: for (j = sj; j < sj+ny; j++) {
882: for (i = si; i < si+nx; i++) {
883: PetscScalar centroid_x = _prop_coords[j][i].x; /* centroids of cell */
884: PetscScalar centroid_y = _prop_coords[j][i].y;
885: PETSC_UNUSED PetscScalar coord_x,coord_y;
888: if (coefficient_structure == 0) { /* isotropic */
889: PetscScalar opts_E,opts_nu;
891: opts_E = 1.0;
892: opts_nu = 0.33;
893: PetscOptionsGetScalar(NULL,"-iso_E",&opts_E,&flg);
894: PetscOptionsGetScalar(NULL,"-iso_nu",&opts_nu,&flg);
896: for (p = 0; p < GAUSS_POINTS; p++) {
897: element_props[j][i].E[p] = opts_E;
898: element_props[j][i].nu[p] = opts_nu;
900: element_props[j][i].fx[p] = 0.0;
901: element_props[j][i].fy[p] = 0.0;
902: }
903: } else if (coefficient_structure == 1) { /* step */
904: PetscScalar opts_E0,opts_nu0,opts_xc;
905: PetscScalar opts_E1,opts_nu1;
907: opts_E0 = opts_E1 = 1.0;
908: opts_nu0 = opts_nu1 = 0.333;
909: opts_xc = 0.5;
910: PetscOptionsGetScalar(NULL,"-step_E0",&opts_E0,&flg);
911: PetscOptionsGetScalar(NULL,"-step_nu0",&opts_nu0,&flg);
912: PetscOptionsGetScalar(NULL,"-step_E1",&opts_E1,&flg);
913: PetscOptionsGetScalar(NULL,"-step_nu1",&opts_nu1,&flg);
914: PetscOptionsGetScalar(NULL,"-step_xc",&opts_xc,&flg);
916: for (p = 0; p < GAUSS_POINTS; p++) {
917: coord_x = centroid_x;
918: coord_y = centroid_y;
919: if (use_gp_coords) {
920: coord_x = element_props[j][i].gp_coords[2*p];
921: coord_y = element_props[j][i].gp_coords[2*p+1];
922: }
924: element_props[j][i].E[p] = opts_E0;
925: element_props[j][i].nu[p] = opts_nu0;
926: if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
927: element_props[j][i].E[p] = opts_E1;
928: element_props[j][i].nu[p] = opts_nu1;
929: }
931: element_props[j][i].fx[p] = 0.0;
932: element_props[j][i].fy[p] = 0.0;
933: }
934: } else if (coefficient_structure == 2) { /* brick */
935: PetscReal values_E[10];
936: PetscReal values_nu[10];
937: PetscInt nbricks,maxnbricks;
938: PetscInt index,span;
939: PetscInt jj;
941: flg = PETSC_FALSE;
942: maxnbricks = 10;
943: PetscOptionsGetRealArray(NULL, "-brick_E",values_E,&maxnbricks,&flg);
944: nbricks = maxnbricks;
945: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of E values for each brick");
947: flg = PETSC_FALSE;
948: maxnbricks = 10;
949: PetscOptionsGetRealArray(NULL, "-brick_nu",values_nu,&maxnbricks,&flg);
950: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of nu values for each brick");
951: if (maxnbricks != nbricks) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply equal numbers of values for E and nu");
953: span = 1;
954: PetscOptionsGetInt(NULL,"-brick_span",&span,&flg);
956: /* cycle through the indices so that no two material properties are repeated in lines of x or y */
957: jj = (j/span)%nbricks;
958: index = (jj+i/span)%nbricks;
959: /*printf("j=%d: index = %d \n", j,index); */
961: for (p = 0; p < GAUSS_POINTS; p++) {
962: element_props[j][i].E[p] = values_E[index];
963: element_props[j][i].nu[p] = values_nu[index];
964: }
965: } else if (coefficient_structure == 3) { /* sponge */
966: PetscScalar opts_E0,opts_nu0;
967: PetscScalar opts_E1,opts_nu1;
968: PetscInt opts_t,opts_w;
969: PetscInt ii,jj,ci,cj;
971: opts_E0 = opts_E1 = 1.0;
972: opts_nu0 = opts_nu1 = 0.333;
973: PetscOptionsGetScalar(NULL,"-sponge_E0",&opts_E0,&flg);
974: PetscOptionsGetScalar(NULL,"-sponge_nu0",&opts_nu0,&flg);
975: PetscOptionsGetScalar(NULL,"-sponge_E1",&opts_E1,&flg);
976: PetscOptionsGetScalar(NULL,"-sponge_nu1",&opts_nu1,&flg);
978: opts_t = opts_w = 1;
979: PetscOptionsGetInt(NULL,"-sponge_t",&opts_t,&flg);
980: PetscOptionsGetInt(NULL,"-sponge_w",&opts_w,&flg);
982: ii = (i)/(opts_t+opts_w+opts_t);
983: jj = (j)/(opts_t+opts_w+opts_t);
985: ci = i - ii*(opts_t+opts_w+opts_t);
986: cj = j - jj*(opts_t+opts_w+opts_t);
988: for (p = 0; p < GAUSS_POINTS; p++) {
989: element_props[j][i].E[p] = opts_E0;
990: element_props[j][i].nu[p] = opts_nu0;
991: }
992: if ((ci >= opts_t) && (ci < opts_t+opts_w)) {
993: if ((cj >= opts_t) && (cj < opts_t+opts_w)) {
994: for (p = 0; p < GAUSS_POINTS; p++) {
995: element_props[j][i].E[p] = opts_E1;
996: element_props[j][i].nu[p] = opts_nu1;
997: }
998: }
999: }
1001: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1002: }
1003: }
1004: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
1006: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1008: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1009: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1010: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
1012: PetscOptionsGetBool(NULL,"-no_view",&no_view,NULL);
1013: if (!no_view) {
1014: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for elasticity eqn.","properties");
1015: DMDACoordViewGnuplot2d(elas_da,"mesh");
1016: }
1018: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1019: DMSetMatType(elas_da,MATAIJ);
1020: DMCreateMatrix(elas_da,&A);
1021: DMGetCoordinates(elas_da,&vel_coords);
1022: MatNullSpaceCreateRigidBody(vel_coords,&matnull);
1023: MatSetNearNullSpace(A,matnull);
1024: MatNullSpaceDestroy(&matnull);
1025: MatGetVecs(A,&f,&X);
1027: /* assemble A11 */
1028: MatZeroEntries(A);
1029: VecZeroEntries(f);
1031: AssembleA_Elasticity(A,elas_da,da_prop,properties);
1032: /* build force vector */
1033: AssembleF_Elasticity(f,elas_da,da_prop,properties);
1036: KSPCreate(PETSC_COMM_WORLD,&ksp_E);
1037: KSPSetOptionsPrefix(ksp_E,"elas_"); /* elasticity */
1039: PetscOptionsGetBool(NULL,"-use_nonsymbc",&use_nonsymbc,&flg);
1040: /* solve */
1041: if (!use_nonsymbc) {
1042: Mat AA;
1043: Vec ff,XX;
1044: IS is;
1045: VecScatter scat;
1047: DMDABCApplySymmetricCompression(elas_da,A,f,&is,&AA,&ff);
1048: VecDuplicate(ff,&XX);
1050: KSPSetOperators(ksp_E,AA,AA);
1051: KSPSetFromOptions(ksp_E);
1053: KSPSolve(ksp_E,ff,XX);
1055: /* push XX back into X */
1056: DMDABCApplyCompression(elas_da,NULL,X);
1058: VecScatterCreate(XX,NULL,X,is,&scat);
1059: VecScatterBegin(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
1060: VecScatterEnd(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
1061: VecScatterDestroy(&scat);
1063: MatDestroy(&AA);
1064: VecDestroy(&ff);
1065: VecDestroy(&XX);
1066: ISDestroy(&is);
1067: } else {
1068: DMDABCApplyCompression(elas_da,A,f);
1070: KSPSetOperators(ksp_E,A,A);
1071: KSPSetFromOptions(ksp_E);
1073: KSPSolve(ksp_E,f,X);
1074: }
1076: if (!no_view) {DMDAViewGnuplot2d(elas_da,X,"Displacement solution for elasticity eqn.","X");}
1077: KSPDestroy(&ksp_E);
1079: VecDestroy(&X);
1080: VecDestroy(&f);
1081: MatDestroy(&A);
1083: DMDestroy(&elas_da);
1084: DMDestroy(&da_prop);
1086: VecDestroy(&properties);
1087: VecDestroy(&l_properties);
1088: return(0);
1089: }
1093: int main(int argc,char **args)
1094: {
1096: PetscInt mx,my;
1098: PetscInitialize(&argc,&args,(char*)0,help);
1100: mx = my = 10;
1101: PetscOptionsGetInt(NULL,"-mx",&mx,NULL);
1102: PetscOptionsGetInt(NULL,"-my",&my,NULL);
1104: solve_elasticity_2d(mx,my);
1106: PetscFinalize();
1107: return 0;
1108: }
1110: /* -------------------------- helpers for boundary conditions -------------------------------- */
1114: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1115: {
1116: DM cda;
1117: Vec coords;
1118: PetscInt si,sj,nx,ny,i,j;
1119: PetscInt M,N;
1120: DMDACoor2d **_coords;
1121: const PetscInt *g_idx;
1122: PetscInt *bc_global_ids;
1123: PetscScalar *bc_vals;
1124: PetscInt nbcs;
1125: PetscInt n_dofs;
1126: PetscErrorCode ierr;
1127: ISLocalToGlobalMapping ltogm;
1130: /* enforce bc's */
1131: DMGetLocalToGlobalMapping(da,<ogm);
1132: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1134: DMGetCoordinateDM(da,&cda);
1135: DMGetCoordinatesLocal(da,&coords);
1136: DMDAVecGetArray(cda,coords,&_coords);
1137: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1138: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1140: /* --- */
1142: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1143: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1145: /* init the entries to -1 so VecSetValues will ignore them */
1146: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1148: i = nx-1;
1149: for (j = 0; j < ny; j++) {
1150: PetscInt local_id;
1151: PETSC_UNUSED PetscScalar coordx,coordy;
1153: local_id = i+j*nx;
1155: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1157: coordx = _coords[j+sj][i+si].x;
1158: coordy = _coords[j+sj][i+si].y;
1160: bc_vals[j] = bc_val;
1161: }
1162: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1163: nbcs = 0;
1164: if ((si+nx) == (M)) nbcs = ny;
1166: if (b) {
1167: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1168: VecAssemblyBegin(b);
1169: VecAssemblyEnd(b);
1170: }
1171: if (A) {
1172: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1173: }
1175: PetscFree(bc_vals);
1176: PetscFree(bc_global_ids);
1178: DMDAVecRestoreArray(cda,coords,&_coords);
1179: return(0);
1180: }
1184: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1185: {
1186: DM cda;
1187: Vec coords;
1188: PetscInt si,sj,nx,ny,i,j;
1189: PetscInt M,N;
1190: DMDACoor2d **_coords;
1191: const PetscInt *g_idx;
1192: PetscInt *bc_global_ids;
1193: PetscScalar *bc_vals;
1194: PetscInt nbcs;
1195: PetscInt n_dofs;
1196: PetscErrorCode ierr;
1197: ISLocalToGlobalMapping ltogm;
1200: /* enforce bc's */
1201: DMGetLocalToGlobalMapping(da,<ogm);
1202: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1204: DMGetCoordinateDM(da,&cda);
1205: DMGetCoordinatesLocal(da,&coords);
1206: DMDAVecGetArray(cda,coords,&_coords);
1207: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1208: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1210: /* --- */
1212: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1213: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1215: /* init the entries to -1 so VecSetValues will ignore them */
1216: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1218: i = 0;
1219: for (j = 0; j < ny; j++) {
1220: PetscInt local_id;
1221: PETSC_UNUSED PetscScalar coordx,coordy;
1223: local_id = i+j*nx;
1225: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1227: coordx = _coords[j+sj][i+si].x;
1228: coordy = _coords[j+sj][i+si].y;
1230: bc_vals[j] = bc_val;
1231: }
1232: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1233: nbcs = 0;
1234: if (si == 0) nbcs = ny;
1236: if (b) {
1237: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1238: VecAssemblyBegin(b);
1239: VecAssemblyEnd(b);
1240: }
1241: if (A) {
1242: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1243: }
1245: PetscFree(bc_vals);
1246: PetscFree(bc_global_ids);
1248: DMDAVecRestoreArray(cda,coords,&_coords);
1249: return(0);
1250: }
1254: static PetscErrorCode DMDABCApplyCompression(DM elas_da,Mat A,Vec f)
1255: {
1259: BCApply_EAST(elas_da,0,-1.0,A,f);
1260: BCApply_EAST(elas_da,1, 0.0,A,f);
1261: BCApply_WEST(elas_da,0,1.0,A,f);
1262: BCApply_WEST(elas_da,1,0.0,A,f);
1263: return(0);
1264: }
1268: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff)
1269: {
1271: PetscInt start,end,m;
1272: PetscInt *unconstrained;
1273: PetscInt cnt,i;
1274: Vec x;
1275: PetscScalar *_x;
1276: IS is;
1277: VecScatter scat;
1280: /* push bc's into f and A */
1281: VecDuplicate(f,&x);
1282: BCApply_EAST(elas_da,0,-1.0,A,x);
1283: BCApply_EAST(elas_da,1, 0.0,A,x);
1284: BCApply_WEST(elas_da,0,1.0,A,x);
1285: BCApply_WEST(elas_da,1,0.0,A,x);
1287: /* define which dofs are not constrained */
1288: VecGetLocalSize(x,&m);
1289: PetscMalloc(sizeof(PetscInt)*m,&unconstrained);
1290: VecGetOwnershipRange(x,&start,&end);
1291: VecGetArray(x,&_x);
1292: cnt = 0;
1293: for (i = 0; i < m; i++) {
1294: PetscReal val;
1296: val = PetscRealPart(_x[i]);
1297: if (fabs(val) < 0.1) {
1298: unconstrained[cnt] = start + i;
1299: cnt++;
1300: }
1301: }
1302: VecRestoreArray(x,&_x);
1304: ISCreateGeneral(PETSC_COMM_WORLD,cnt,unconstrained,PETSC_COPY_VALUES,&is);
1305: PetscFree(unconstrained);
1307: /* define correction for dirichlet in the rhs */
1308: MatMult(A,x,f);
1309: VecScale(f,-1.0);
1311: /* get new matrix */
1312: MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,AA);
1313: /* get new vector */
1314: MatGetVecs(*AA,NULL,ff);
1316: VecScatterCreate(f,is,*ff,NULL,&scat);
1317: VecScatterBegin(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1318: VecScatterEnd(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1320: { /* Constrain near-null space */
1321: PetscInt nvecs;
1322: const Vec *vecs;
1323: Vec *uvecs;
1324: PetscBool has_const;
1325: MatNullSpace mnull,unull;
1326: MatGetNearNullSpace(A,&mnull);
1327: MatNullSpaceGetVecs(mnull,&has_const,&nvecs,&vecs);
1328: VecDuplicateVecs(*ff,nvecs,&uvecs);
1329: for (i=0; i<nvecs; i++) {
1330: VecScatterBegin(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1331: VecScatterEnd(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1332: }
1333: MatNullSpaceCreate(PetscObjectComm((PetscObject)A),has_const,nvecs,uvecs,&unull);
1334: MatSetNearNullSpace(*AA,unull);
1335: MatNullSpaceDestroy(&unull);
1336: for (i=0; i<nvecs; i++) {
1337: VecDestroy(&uvecs[i]);
1338: }
1339: PetscFree(uvecs);
1340: }
1342: VecScatterDestroy(&scat);
1344: *dofs = is;
1345: VecDestroy(&x);
1346: return(0);
1347: }