Actual source code: ex49.c
petsc-3.3-p7 2013-05-11
1: static char help[] =
2: " Solves the compressible plane strain elasticity equations in 2d on the unit domain using Q1 finite elements. \n\
3: Material properties E (Youngs moduls) and nu (Poisson ratio) may vary as a function of space. \n\
4: The model utilisse boundary conditions which produce compression in the x direction. \n\
5: Options: \n\
6: -mx : number elements in x-direciton \n\
7: -my : number elements in y-direciton \n\
8: -c_str : indicates the structure of the coefficients to use. \n\
9: -c_str 0 => Setup for an isotropic material with constant coefficients. \n\
10: Parameters: \n\
11: -iso_E : Youngs modulus \n\
12: -iso_nu : Poisson ratio \n\
13: -c_str 1 => Setup for a step function in the material properties in x. \n\
14: Parameters: \n\
15: -step_E0 : Youngs modulus to the left of the step \n\
16: -step_nu0 : Poisson ratio to the left of the step \n\
17: -step_E1 : Youngs modulus to the right of the step \n\
18: -step_n1 : Poisson ratio to the right of the step \n\
19: -step_xc : x coordinate of the step \n\
20: -c_str 2 => Setup for a checkerboard material with alternating properties. \n\
21: Repeats the following pattern throughout the domain. For example with 4 materials specified, we would heve \n\
22: -------------------------\n\
23: | D | A | B | C |\n\
24: ------|-----|-----|------\n\
25: | C | D | A | B |\n\
26: ------|-----|-----|------\n\
27: | B | C | D | A |\n\
28: ------|-----|-----|------\n\
29: | A | B | C | D |\n\
30: -------------------------\n\
31: \n\
32: Parameters: \n\
33: -brick_E : a comma seperated list of Young's modulii \n\
34: -brick_nu : a comma seperated list of Poisson ratio's \n\
35: -brick_span : the number of elements in x and y each brick will span \n\
36: -c_str 3 => Setup for a sponge-like material with alternating properties. \n\
37: Repeats the following pattern throughout the domain \n\
38: -----------------------------\n\
39: | [background] |\n\
40: | E0,nu0 |\n\
41: | ----------------- |\n\
42: | | [inclusion] | |\n\
43: | | E1,nu1 | |\n\
44: | | | |\n\
45: | | <---- w ----> | |\n\
46: | | | |\n\
47: | | | |\n\
48: | ----------------- |\n\
49: | |\n\
50: | |\n\
51: -----------------------------\n\
52: <-------- t + w + t ------->\n\
53: \n\
54: Parameters: \n\
55: -sponge_E0 : Youngs moduls of the surrounding material \n\
56: -sponge_E1 : Youngs moduls of the inclusio \n\
57: -sponge_nu0 : Poisson ratio of the surrounding material \n\
58: -sponge_nu1 : Poisson ratio of the inclusio \n\
59: -sponge_t : the number of elements defining the border around each inclusion \n\
60: -sponge_w : the number of elements in x and y each inclusion will span\n\
61: -use_gp_coords : Evaluate the Youngs modulus, Poisson ratio and the body force at the global coordinates of the quadrature points.\n\
62: By default, E, nu and the body force are evaulated at the element center and applied as a constant over the entire element.\n\
63: -use_nonsymbc : Option to use non-symmetric boundary condition imposition. This choice will use less memory.";
65: /* Contributed by Dave May */
67: #include <petscksp.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) {
176: *det_J = J;
177: }
178: }
180: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
181: {
182: *ngp = 4;
183: gp_xi[0][0] = -0.57735026919;gp_xi[0][1] = -0.57735026919;
184: gp_xi[1][0] = -0.57735026919;gp_xi[1][1] = 0.57735026919;
185: gp_xi[2][0] = 0.57735026919;gp_xi[2][1] = 0.57735026919;
186: gp_xi[3][0] = 0.57735026919;gp_xi[3][1] = -0.57735026919;
187: gp_weight[0] = 1.0;
188: gp_weight[1] = 1.0;
189: gp_weight[2] = 1.0;
190: gp_weight[3] = 1.0;
191: }
194: /* procs to the left claim the ghost node as their element */
197: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
198: {
200: PetscInt m,n,p,M,N,P;
201: PetscInt sx,sy,sz;
204: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
205: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
207: if (mxl != PETSC_NULL) {
208: *mxl = m;
209: if ((sx+m) == M) { /* last proc */
210: *mxl = m-1;
211: }
212: }
213: if (myl != PETSC_NULL) {
214: *myl = n;
215: if ((sy+n) == N) { /* last proc */
216: *myl = n-1;
217: }
218: }
219: if (mzl != PETSC_NULL) {
220: *mzl = p;
221: if ((sz+p) == P) { /* last proc */
222: *mzl = p-1;
223: }
224: }
226: return(0);
227: }
231: static PetscErrorCode DMDAGetElementCorners(DM da,
232: PetscInt *sx,PetscInt *sy,PetscInt *sz,
233: PetscInt *mx,PetscInt *my,PetscInt *mz)
234: {
236: PetscInt si,sj,sk;
239: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
241: if (sx != PETSC_NULL) {
242: *sx = si;
243: if (si != 0) {
244: *sx = si+1;
245: }
246: }
247: if (sy != PETSC_NULL) {
248: *sy = sj;
249: if (sj != 0) {
250: *sy = sj+1;
251: }
252: }
254: if (sk != PETSC_NULL) {
255: *sz = sk;
256: if (sk != 0) {
257: *sz = sk+1;
258: }
259: }
261: DMDAGetLocalElementSize(da,mx,my,mz);
263: return(0);
264: }
268: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
269: {
271: PetscMPIInt rank;
272: PetscInt proc_I,proc_J;
273: PetscInt cpu_x,cpu_y;
274: PetscInt local_mx,local_my;
275: Vec vlx,vly;
276: PetscInt *LX,*LY,i;
277: PetscScalar *_a;
278: Vec V_SEQ;
279: VecScatter ctx;
282: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
284: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
286: proc_J = rank/cpu_x;
287: proc_I = rank-cpu_x*proc_J;
289: PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);
290: PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);
292: DMDAGetLocalElementSize(da,&local_mx,&local_my,PETSC_NULL);
293: VecCreate(PETSC_COMM_WORLD,&vlx);
294: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
295: VecSetFromOptions(vlx);
297: VecCreate(PETSC_COMM_WORLD,&vly);
298: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
299: VecSetFromOptions(vly);
301: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
302: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
303: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
304: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
308: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
309: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
310: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
311: VecGetArray(V_SEQ,&_a);
312: for (i = 0; i < cpu_x; i++) {
313: LX[i] = (PetscInt)PetscRealPart(_a[i]);
314: }
315: VecRestoreArray(V_SEQ,&_a);
316: VecScatterDestroy(&ctx);
317: VecDestroy(&V_SEQ);
319: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
320: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
321: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
322: VecGetArray(V_SEQ,&_a);
323: for (i = 0; i < cpu_y; i++) {
324: LY[i] = (PetscInt)PetscRealPart(_a[i]);
325: }
326: VecRestoreArray(V_SEQ,&_a);
327: VecScatterDestroy(&ctx);
328: VecDestroy(&V_SEQ);
330: *_lx = LX;
331: *_ly = LY;
333: VecDestroy(&vlx);
334: VecDestroy(&vly);
336: return(0);
337: }
341: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
342: {
343: DM cda;
344: Vec coords;
345: DMDACoor2d **_coords;
346: PetscInt si,sj,nx,ny,i,j;
347: FILE *fp;
348: char fname[PETSC_MAX_PATH_LEN];
349: PetscMPIInt rank;
353: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
354: PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
355: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
356: if (fp == PETSC_NULL) {
357: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
358: }
360: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
362: DMDAGetCoordinateDA(da,&cda);
363: DMDAGetGhostedCoordinates(da,&coords);
364: DMDAVecGetArray(cda,coords,&_coords);
365: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
366: for (j = sj; j < sj+ny-1; j++) {
367: for (i = si; i < si+nx-1; i++) {
368: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
369: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i].x),PetscRealPart(_coords[j+1][i].y));
370: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i+1].x),PetscRealPart(_coords[j+1][i+1].y));
371: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i+1].x),PetscRealPart(_coords[j][i+1].y));
372: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
373: }
374: }
375: DMDAVecRestoreArray(cda,coords,&_coords);
377: PetscFClose(PETSC_COMM_SELF,fp);
378: return(0);
379: }
383: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
384: {
385: DM cda;
386: Vec coords,local_fields;
387: DMDACoor2d **_coords;
388: FILE *fp;
389: char fname[PETSC_MAX_PATH_LEN];
390: const char *field_name;
391: PetscMPIInt rank;
392: PetscInt si,sj,nx,ny,i,j;
393: PetscInt n_dofs,d;
394: PetscScalar *_fields;
398: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
399: PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
400: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
401: if (fp == PETSC_NULL) {
402: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
403: }
405: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
406: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
407: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
408: for (d = 0; d < n_dofs; d++) {
409: DMDAGetFieldName(da,d,&field_name);
410: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
411: }
412: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
415: DMDAGetCoordinateDA(da,&cda);
416: DMDAGetGhostedCoordinates(da,&coords);
417: DMDAVecGetArray(cda,coords,&_coords);
418: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
420: DMCreateLocalVector(da,&local_fields);
421: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
422: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
423: VecGetArray(local_fields,&_fields);
425: for (j = sj; j < sj+ny; j++) {
426: for (i = si; i < si+nx; i++) {
427: PetscScalar coord_x,coord_y;
428: PetscScalar field_d;
430: coord_x = _coords[j][i].x;
431: coord_y = _coords[j][i].y;
433: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
434: for (d = 0; d < n_dofs; d++) {
435: field_d = _fields[ n_dofs*((i-si)+(j-sj)*(nx))+d ];
436: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",PetscRealPart(field_d));
437: }
438: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
439: }
440: }
441: VecRestoreArray(local_fields,&_fields);
442: VecDestroy(&local_fields);
444: DMDAVecRestoreArray(cda,coords,&_coords);
446: PetscFClose(PETSC_COMM_SELF,fp);
447: return(0);
448: }
452: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
453: {
454: DM cda;
455: Vec local_fields;
456: FILE *fp;
457: char fname[PETSC_MAX_PATH_LEN];
458: const char *field_name;
459: PetscMPIInt rank;
460: PetscInt si,sj,nx,ny,i,j,p;
461: PetscInt n_dofs,d;
462: GaussPointCoefficients **_coefficients;
463: PetscErrorCode ierr;
466: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
467: PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
468: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
469: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
471: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
472: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
473: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
474: for (d = 0; d < n_dofs; d++) {
475: DMDAGetFieldName(da,d,&field_name);
476: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
477: }
478: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
481: DMDAGetCoordinateDA(da,&cda);
482: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
484: DMCreateLocalVector(da,&local_fields);
485: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
486: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
487: DMDAVecGetArray(da,local_fields,&_coefficients);
489: for (j = sj; j < sj+ny; j++) {
490: for (i = si; i < si+nx; i++) {
491: PetscScalar coord_x,coord_y;
493: for (p = 0; p < GAUSS_POINTS; p++) {
494: coord_x = _coefficients[j][i].gp_coords[2*p];
495: coord_y = _coefficients[j][i].gp_coords[2*p+1];
497: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
499: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e %1.6e",
500: PetscRealPart(_coefficients[j][i].E[p]),PetscRealPart(_coefficients[j][i].nu[p]),
501: PetscRealPart(_coefficients[j][i].fx[p]),PetscRealPart(_coefficients[j][i].fy[p]));
502: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
503: }
504: }
505: }
506: DMDAVecRestoreArray(da,local_fields,&_coefficients);
507: VecDestroy(&local_fields);
509: PetscFClose(PETSC_COMM_SELF,fp);
510: return(0);
511: }
513: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar E[],PetscScalar nu[])
514: {
515: PetscInt ngp;
516: PetscScalar gp_xi[GAUSS_POINTS][2];
517: PetscScalar gp_weight[GAUSS_POINTS];
518: PetscInt p,i,j,k,l;
519: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
520: PetscScalar J_p;
521: PetscScalar B[3][U_DOFS*NODES_PER_EL];
522: PetscScalar prop_E,prop_nu,factor,constit_D[3][3];
524: /* define quadrature rule */
525: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
527: /* evaluate integral */
528: for (p = 0; p < ngp; p++) {
529: ConstructQ12D_GNi(gp_xi[p],GNi_p);
530: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
532: for (i = 0; i < NODES_PER_EL; i++) {
533: PetscScalar d_dx_i = GNx_p[0][i];
534: PetscScalar d_dy_i = GNx_p[1][i];
536: B[0][2*i] = d_dx_i; B[0][2*i+1] = 0.0;
537: B[1][2*i] = 0.0; B[1][2*i+1] = d_dy_i;
538: B[2][2*i] = d_dy_i; B[2][2*i+1] = d_dx_i;
539: }
541: /* form D for the quadrature point */
542: prop_E = E[p];
543: prop_nu = nu[p];
544: factor = prop_E / ( (1.0+prop_nu)*(1.0-2.0*prop_nu) );
545: constit_D[0][0] = 1.0-prop_nu; constit_D[0][1] = prop_nu; constit_D[0][2] = 0.0;
546: constit_D[1][0] = prop_nu; constit_D[1][1] = 1.0-prop_nu; constit_D[1][2] = 0.0;
547: constit_D[2][0] = 0.0; constit_D[2][1] = 0.0; constit_D[2][2] = 0.5*(1.0-2.0*prop_nu);
548: for (i = 0; i < 3; i++) {
549: for (j = 0; j < 3; j++) {
550: constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
551: }
552: }
554: /* form Bt tildeD B */
555: /*
556: Ke_ij = Bt_ik . D_kl . B_lj
557: = B_ki . D_kl . B_lj
558: */
559: for (i = 0; i < 8; i++) {
560: for (j = 0; j < 8; j++) {
562: for (k = 0; k < 3; k++) {
563: for (l = 0; l < 3; l++) {
564: Ke[8*i+j] = Ke[8*i+j] + B[k][i] * constit_D[k][l] * B[l][j];
565: }
566: }
567: }
568: }
570: } /* end quadrature */
571: }
573: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
574: {
575: PetscInt ngp;
576: PetscScalar gp_xi[GAUSS_POINTS][2];
577: PetscScalar gp_weight[GAUSS_POINTS];
578: PetscInt p,i;
579: PetscScalar Ni_p[NODES_PER_EL];
580: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
581: PetscScalar J_p,fac;
583: /* define quadrature rule */
584: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
586: /* evaluate integral */
587: for (p = 0; p < ngp; p++) {
588: ConstructQ12D_Ni(gp_xi[p],Ni_p);
589: ConstructQ12D_GNi(gp_xi[p],GNi_p);
590: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
591: fac = gp_weight[p]*J_p;
593: for (i = 0; i < NODES_PER_EL; i++) {
594: Fe[NSD*i ] += fac*Ni_p[i]*fx[p];
595: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
596: }
597: }
598: }
600: /*
601: i,j are the element indices
602: The unknown is a vector quantity.
603: The s[].c is used to indicate the degree of freedom.
604: */
607: static PetscErrorCode DMDAGetElementEqnums_u(MatStencil s_u[],PetscInt i,PetscInt j)
608: {
610: /* displacement */
611: /* node 0 */
612: s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0; /* Ux0 */
613: s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1; /* Uy0 */
615: /* node 1 */
616: s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0; /* Ux1 */
617: s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1; /* Uy1 */
619: /* node 2 */
620: s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0; /* Ux2 */
621: s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1; /* Uy2 */
623: /* node 3 */
624: s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0; /* Ux3 */
625: s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1; /* Uy3 */
627: return(0);
628: }
632: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
633: {
635: /* get coords for the element */
636: el_coords[NSD*0+0] = _coords[ej ][ei ].x; el_coords[NSD*0+1] = _coords[ej ][ei ].y;
637: el_coords[NSD*1+0] = _coords[ej+1][ei ].x; el_coords[NSD*1+1] = _coords[ej+1][ei ].y;
638: el_coords[NSD*2+0] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
639: el_coords[NSD*3+0] = _coords[ej ][ei+1].x; el_coords[NSD*3+1] = _coords[ej ][ei+1].y;
641: return(0);
642: }
646: static PetscErrorCode AssembleA_Elasticity(Mat A,DM elas_da,DM properties_da,Vec properties)
647: {
648: DM cda;
649: Vec coords;
650: DMDACoor2d **_coords;
651: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
652: PetscInt sex,sey,mx,my;
653: PetscInt ei,ej;
654: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
655: PetscScalar el_coords[NODES_PER_EL*NSD];
656: Vec local_properties;
657: GaussPointCoefficients **props;
658: PetscScalar *prop_E,*prop_nu;
659: PetscErrorCode ierr;
662: /* setup for coords */
663: DMDAGetCoordinateDA(elas_da,&cda);
664: DMDAGetGhostedCoordinates(elas_da,&coords);
665: DMDAVecGetArray(cda,coords,&_coords);
667: /* setup for coefficients */
668: DMCreateLocalVector(properties_da,&local_properties);
669: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
670: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
671: DMDAVecGetArray(properties_da,local_properties,&props);
673: DMDAGetElementCorners(elas_da,&sex,&sey,0,&mx,&my,0);
674: for (ej = sey; ej < sey+my; ej++) {
675: for (ei = sex; ei < sex+mx; ei++) {
676: /* get coords for the element */
677: GetElementCoords(_coords,ei,ej,el_coords);
679: /* get coefficients for the element */
680: prop_E = props[ej][ei].E;
681: prop_nu = props[ej][ei].nu;
683: /* initialise element stiffness matrix */
684: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
686: /* form element stiffness matrix */
687: FormStressOperatorQ1(Ae,el_coords,prop_E,prop_nu);
689: /* insert element matrix into global matrix */
690: DMDAGetElementEqnums_u(u_eqn,ei,ej);
691: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
692: }
693: }
694: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
695: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
697: DMDAVecRestoreArray(cda,coords,&_coords);
699: DMDAVecRestoreArray(properties_da,local_properties,&props);
700: VecDestroy(&local_properties);
702: return(0);
703: }
708: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F,MatStencil u_eqn[],PetscScalar Fe_u[])
709: {
710: PetscInt n;
713: for (n = 0; n < 4; n++) {
714: 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 ];
715: 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];
716: }
717: return(0);
718: }
722: static PetscErrorCode AssembleF_Elasticity(Vec F,DM elas_da,DM properties_da,Vec properties)
723: {
724: DM cda;
725: Vec coords;
726: DMDACoor2d **_coords;
727: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
728: PetscInt sex,sey,mx,my;
729: PetscInt ei,ej;
730: PetscScalar Fe[NODES_PER_EL*U_DOFS];
731: PetscScalar el_coords[NODES_PER_EL*NSD];
732: Vec local_properties;
733: GaussPointCoefficients **props;
734: PetscScalar *prop_fx,*prop_fy;
735: Vec local_F;
736: ElasticityDOF **ff;
737: PetscErrorCode ierr;
740: /* setup for coords */
741: DMDAGetCoordinateDA(elas_da,&cda);
742: DMDAGetGhostedCoordinates(elas_da,&coords);
743: DMDAVecGetArray(cda,coords,&_coords);
745: /* setup for coefficients */
746: DMGetLocalVector(properties_da,&local_properties);
747: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
748: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
749: DMDAVecGetArray(properties_da,local_properties,&props);
751: /* get acces to the vector */
752: DMGetLocalVector(elas_da,&local_F);
753: VecZeroEntries(local_F);
754: DMDAVecGetArray(elas_da,local_F,&ff);
757: DMDAGetElementCorners(elas_da,&sex,&sey,0,&mx,&my,0);
758: for (ej = sey; ej < sey+my; ej++) {
759: for (ei = sex; ei < sex+mx; ei++) {
760: /* get coords for the element */
761: GetElementCoords(_coords,ei,ej,el_coords);
763: /* get coefficients for the element */
764: prop_fx = props[ej][ei].fx;
765: prop_fy = props[ej][ei].fy;
767: /* initialise element stiffness matrix */
768: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
770: /* form element stiffness matrix */
771: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
773: /* insert element matrix into global matrix */
774: DMDAGetElementEqnums_u(u_eqn,ei,ej);
776: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,Fe);
777: }
778: }
780: DMDAVecRestoreArray(elas_da,local_F,&ff);
781: DMLocalToGlobalBegin(elas_da,local_F,ADD_VALUES,F);
782: DMLocalToGlobalEnd(elas_da,local_F,ADD_VALUES,F);
783: DMRestoreLocalVector(elas_da,&local_F);
785: DMDAVecRestoreArray(cda,coords,&_coords);
787: DMDAVecRestoreArray(properties_da,local_properties,&props);
788: DMRestoreLocalVector(properties_da,&local_properties);
789: return(0);
790: }
794: static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
795: {
796: DM elas_da,da_prop;
797: PetscInt u_dof,dof,stencil_width;
798: Mat A;
799: PetscInt mxl,myl;
800: DM prop_cda,vel_cda;
801: Vec prop_coords,vel_coords;
802: PetscInt si,sj,nx,ny,i,j,p;
803: Vec f,X;
804: PetscInt prop_dof,prop_stencil_width;
805: Vec properties,l_properties;
806: MatNullSpace matnull;
807: PetscReal dx,dy;
808: PetscInt M,N;
809: DMDACoor2d **_prop_coords,**_vel_coords;
810: GaussPointCoefficients **element_props;
811: KSP ksp_E;
812: PetscInt coefficient_structure = 0;
813: PetscInt cpu_x,cpu_y,*lx = PETSC_NULL,*ly = PETSC_NULL;
814: PetscBool use_gp_coords = PETSC_FALSE;
815: PetscBool use_nonsymbc = PETSC_FALSE;
816: PetscBool no_view = PETSC_FALSE;
817: PetscBool flg;
818: PetscErrorCode ierr;
821: /* Generate the da for velocity and pressure */
822: /*
823: We use Q1 elements for the temperature.
824: FEM has a 9-point stencil (BOX) or connectivity pattern
825: Num nodes in each direction is mx+1, my+1
826: */
827: u_dof = U_DOFS; /* Vx, Vy - velocities */
828: dof = u_dof;
829: stencil_width = 1;
830: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
831: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,PETSC_NULL,PETSC_NULL,&elas_da);
832: DMDASetFieldName(elas_da,0,"Ux");
833: DMDASetFieldName(elas_da,1,"Uy");
835: /* unit box [0,1] x [0,1] */
836: DMDASetUniformCoordinates(elas_da,0.0,1.0,0.0,1.0,PETSC_NULL,PETSC_NULL);
839: /* Generate element properties, we will assume all material properties are constant over the element */
840: /* local number of elements */
841: DMDAGetLocalElementSize(elas_da,&mxl,&myl,PETSC_NULL);
843: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! // */
844: DMDAGetInfo(elas_da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
845: DMDAGetElementOwnershipRanges2d(elas_da,&lx,&ly);
847: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
848: prop_stencil_width = 0;
849: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
850: mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
851: PetscFree(lx);
852: PetscFree(ly);
854: /* define centroid positions */
855: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
856: dx = 1.0/((PetscReal)(M));
857: dy = 1.0/((PetscReal)(N));
859: DMDASetUniformCoordinates(da_prop,0.0+0.5*dx,1.0-0.5*dx,0.0+0.5*dy,1.0-0.5*dy,PETSC_NULL,PETSC_NULL);
861: /* define coefficients */
862: PetscOptionsGetInt(PETSC_NULL,"-c_str",&coefficient_structure,PETSC_NULL);
864: DMCreateGlobalVector(da_prop,&properties);
865: DMCreateLocalVector(da_prop,&l_properties);
866: DMDAVecGetArray(da_prop,l_properties,&element_props);
868: DMDAGetCoordinateDA(da_prop,&prop_cda);
869: DMDAGetGhostedCoordinates(da_prop,&prop_coords);
870: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
872: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
874: DMDAGetCoordinateDA(elas_da,&vel_cda);
875: DMDAGetGhostedCoordinates(elas_da,&vel_coords);
876: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
879: /* interpolate the coordinates */
880: for (j = sj; j < sj+ny; j++) {
881: for (i = si; i < si+nx; i++) {
882: PetscInt ngp;
883: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
884: PetscScalar el_coords[8];
886: GetElementCoords(_vel_coords,i,j,el_coords);
887: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
889: for (p = 0; p < GAUSS_POINTS; p++) {
890: PetscScalar gp_x,gp_y;
891: PetscInt n;
892: PetscScalar xi_p[2],Ni_p[4];
894: xi_p[0] = gp_xi[p][0];
895: xi_p[1] = gp_xi[p][1];
896: ConstructQ12D_Ni(xi_p,Ni_p);
898: gp_x = 0.0;
899: gp_y = 0.0;
900: for (n = 0; n < NODES_PER_EL; n++) {
901: gp_x = gp_x+Ni_p[n]*el_coords[2*n ];
902: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
903: }
904: element_props[j][i].gp_coords[2*p ] = gp_x;
905: element_props[j][i].gp_coords[2*p+1] = gp_y;
906: }
907: }
908: }
910: /* define the coefficients */
911: PetscOptionsGetBool(PETSC_NULL,"-use_gp_coords",&use_gp_coords,&flg);
913: for (j = sj; j < sj+ny; j++) {
914: for (i = si; i < si+nx; i++) {
915: PetscScalar centroid_x = _prop_coords[j][i].x; /* centroids of cell */
916: PetscScalar centroid_y = _prop_coords[j][i].y;
917: PETSC_UNUSED PetscScalar coord_x,coord_y;
920: if (coefficient_structure == 0) { /* isotropic */
921: PetscScalar opts_E,opts_nu;
923: opts_E = 1.0;
924: opts_nu = 0.33;
925: PetscOptionsGetScalar(PETSC_NULL,"-iso_E",&opts_E,&flg);
926: PetscOptionsGetScalar(PETSC_NULL,"-iso_nu",&opts_nu,&flg);
928: for (p = 0; p < GAUSS_POINTS; p++) {
929: element_props[j][i].E[p] = opts_E;
930: element_props[j][i].nu[p] = opts_nu;
932: element_props[j][i].fx[p] = 0.0;
933: element_props[j][i].fy[p] = 0.0;
934: }
935: } else if (coefficient_structure == 1) { /* step */
936: PetscScalar opts_E0,opts_nu0,opts_xc;
937: PetscScalar opts_E1,opts_nu1;
939: opts_E0 = opts_E1 = 1.0;
940: opts_nu0 = opts_nu1 = 0.333;
941: opts_xc = 0.5;
942: PetscOptionsGetScalar(PETSC_NULL,"-step_E0",&opts_E0,&flg);
943: PetscOptionsGetScalar(PETSC_NULL,"-step_nu0",&opts_nu0,&flg);
944: PetscOptionsGetScalar(PETSC_NULL,"-step_E1",&opts_E1,&flg);
945: PetscOptionsGetScalar(PETSC_NULL,"-step_nu1",&opts_nu1,&flg);
946: PetscOptionsGetScalar(PETSC_NULL,"-step_xc",&opts_xc,&flg);
948: for (p = 0; p < GAUSS_POINTS; p++) {
949: coord_x = centroid_x;
950: coord_y = centroid_y;
951: if (use_gp_coords == PETSC_TRUE) {
952: coord_x = element_props[j][i].gp_coords[2*p];
953: coord_y = element_props[j][i].gp_coords[2*p+1];
954: }
956: element_props[j][i].E[p] = opts_E0;
957: element_props[j][i].nu[p] = opts_nu0;
958: if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
959: element_props[j][i].E[p] = opts_E1;
960: element_props[j][i].nu[p] = opts_nu1;
961: }
963: element_props[j][i].fx[p] = 0.0;
964: element_props[j][i].fy[p] = 0.0;
965: }
966: } else if (coefficient_structure == 2) { /* brick */
967: PetscReal values_E[10];
968: PetscReal values_nu[10];
969: PetscInt nbricks,maxnbricks;
970: PetscInt index,span;
971: PetscInt jj;
973: flg = PETSC_FALSE;
974: maxnbricks = 10;
975: PetscOptionsGetRealArray( PETSC_NULL, "-brick_E",values_E,&maxnbricks,&flg);
976: nbricks = maxnbricks;
977: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of E values for each brick");
979: flg = PETSC_FALSE;
980: maxnbricks = 10;
981: PetscOptionsGetRealArray( PETSC_NULL, "-brick_nu",values_nu,&maxnbricks,&flg);
982: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of nu values for each brick");
983: if (maxnbricks != nbricks) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply equal numbers of values for E and nu");
985: span = 1;
986: PetscOptionsGetInt(PETSC_NULL,"-brick_span",&span,&flg);
988: /* cycle through the indices so that no two material properties are repeated in lines of x or y */
989: jj = (j/span)%nbricks;
990: index = (jj+i/span)%nbricks;
991: /*printf("j=%d: index = %d \n", j,index ); */
993: for (p = 0; p < GAUSS_POINTS; p++) {
994: element_props[j][i].E[p] = values_E[index];
995: element_props[j][i].nu[p] = values_nu[index];
996: }
997: } else if (coefficient_structure == 3) { /* sponge */
998: PetscScalar opts_E0,opts_nu0;
999: PetscScalar opts_E1,opts_nu1;
1000: PetscInt opts_t,opts_w;
1001: PetscInt ii,jj,ci,cj;
1003: opts_E0 = opts_E1 = 1.0;
1004: opts_nu0 = opts_nu1 = 0.333;
1005: PetscOptionsGetScalar(PETSC_NULL,"-sponge_E0",&opts_E0,&flg);
1006: PetscOptionsGetScalar(PETSC_NULL,"-sponge_nu0",&opts_nu0,&flg);
1007: PetscOptionsGetScalar(PETSC_NULL,"-sponge_E1",&opts_E1,&flg);
1008: PetscOptionsGetScalar(PETSC_NULL,"-sponge_nu1",&opts_nu1,&flg);
1010: opts_t = opts_w = 1;
1011: PetscOptionsGetInt(PETSC_NULL,"-sponge_t",&opts_t,&flg);
1012: PetscOptionsGetInt(PETSC_NULL,"-sponge_w",&opts_w,&flg);
1014: ii = (i)/(opts_t+opts_w+opts_t);
1015: jj = (j)/(opts_t+opts_w+opts_t);
1017: ci = i - ii*(opts_t+opts_w+opts_t);
1018: cj = j - jj*(opts_t+opts_w+opts_t);
1020: for (p = 0; p < GAUSS_POINTS; p++) {
1021: element_props[j][i].E[p] = opts_E0;
1022: element_props[j][i].nu[p] = opts_nu0;
1023: }
1024: if ( (ci >= opts_t) && (ci < opts_t+opts_w) ) {
1025: if ( (cj >= opts_t) && (cj < opts_t+opts_w) ) {
1026: for (p = 0; p < GAUSS_POINTS; p++) {
1027: element_props[j][i].E[p] = opts_E1;
1028: element_props[j][i].nu[p] = opts_nu1;
1029: }
1030: }
1031: }
1033: } else {
1034: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1035: }
1036: }
1037: }
1038: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
1040: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1042: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1043: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1044: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
1046: PetscOptionsGetBool(PETSC_NULL,"-no_view",&no_view,PETSC_NULL);
1047: if (!no_view) {
1048: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for elasticity eqn.","properties");
1049: DMDACoordViewGnuplot2d(elas_da,"mesh");
1050: }
1052: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1053: DMCreateMatrix(elas_da,MATAIJ,&A);
1054: DMDAGetCoordinates(elas_da,&vel_coords);
1055: MatNullSpaceCreateRigidBody(vel_coords,&matnull);
1056: MatSetNearNullSpace(A,matnull);
1057: MatNullSpaceDestroy(&matnull);
1058: MatGetVecs(A,&f,&X);
1060: /* assemble A11 */
1061: MatZeroEntries(A);
1062: VecZeroEntries(f);
1064: AssembleA_Elasticity(A,elas_da,da_prop,properties);
1065: /* build force vector */
1066: AssembleF_Elasticity(f,elas_da,da_prop,properties);
1069: KSPCreate(PETSC_COMM_WORLD,&ksp_E);
1070: KSPSetOptionsPrefix(ksp_E,"elas_"); /* elasticity */
1072: PetscOptionsGetBool(PETSC_NULL,"-use_nonsymbc",&use_nonsymbc,&flg);
1073: /* solve */
1074: if (use_nonsymbc == PETSC_FALSE) {
1075: Mat AA;
1076: Vec ff,XX;
1077: IS is;
1078: VecScatter scat;
1080: DMDABCApplySymmetricCompression(elas_da,A,f,&is,&AA,&ff);
1081: VecDuplicate(ff,&XX);
1083: KSPSetOperators(ksp_E,AA,AA,SAME_NONZERO_PATTERN);
1084: KSPSetFromOptions(ksp_E);
1086: KSPSolve(ksp_E,ff,XX);
1088: /* push XX back into X */
1089: DMDABCApplyCompression(elas_da,PETSC_NULL,X);
1091: VecScatterCreate(XX,PETSC_NULL,X,is,&scat);
1092: VecScatterBegin(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
1093: VecScatterEnd(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
1094: VecScatterDestroy(&scat);
1096: MatDestroy(&AA);
1097: VecDestroy(&ff);
1098: VecDestroy(&XX);
1099: ISDestroy(&is);
1100: } else {
1101: DMDABCApplyCompression(elas_da,A,f);
1103: KSPSetOperators(ksp_E,A,A,SAME_NONZERO_PATTERN);
1104: KSPSetFromOptions(ksp_E);
1106: KSPSolve(ksp_E,f,X);
1107: }
1109: if (!no_view) {DMDAViewGnuplot2d(elas_da,X,"Displacement solution for elasticity eqn.","X");}
1110: KSPDestroy(&ksp_E);
1112: VecDestroy(&X);
1113: VecDestroy(&f);
1114: MatDestroy(&A);
1116: DMDestroy(&elas_da);
1117: DMDestroy(&da_prop);
1119: VecDestroy(&properties);
1120: VecDestroy(&l_properties);
1122: return(0);
1123: }
1127: int main(int argc,char **args)
1128: {
1130: PetscInt mx,my;
1132: PetscInitialize(&argc,&args,(char *)0,help);
1134: mx = my = 10;
1135: PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,PETSC_NULL);
1136: PetscOptionsGetInt(PETSC_NULL,"-my",&my,PETSC_NULL);
1138: solve_elasticity_2d(mx,my);
1140: PetscFinalize();
1141: return 0;
1142: }
1144: /* -------------------------- helpers for boundary conditions -------------------------------- */
1148: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1149: {
1150: DM cda;
1151: Vec coords;
1152: PetscInt si,sj,nx,ny,i,j;
1153: PetscInt M,N;
1154: DMDACoor2d **_coords;
1155: PetscInt *g_idx;
1156: PetscInt *bc_global_ids;
1157: PetscScalar *bc_vals;
1158: PetscInt nbcs;
1159: PetscInt n_dofs;
1163: /* enforce bc's */
1164: DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);
1166: DMDAGetCoordinateDA(da,&cda);
1167: DMDAGetGhostedCoordinates(da,&coords);
1168: DMDAVecGetArray(cda,coords,&_coords);
1169: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1170: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1172: /* /// */
1174: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1175: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1177: /* init the entries to -1 so VecSetValues will ignore them */
1178: for (i = 0; i < ny*n_dofs; i++) {
1179: bc_global_ids[i] = -1;
1180: }
1182: i = nx-1;
1183: for (j = 0; j < ny; j++) {
1184: PetscInt local_id;
1185: PETSC_UNUSED PetscScalar coordx,coordy;
1187: local_id = i+j*nx;
1189: bc_global_ids[j] = g_idx[ n_dofs*local_id+d_idx ];
1191: coordx = _coords[j+sj][i+si].x;
1192: coordy = _coords[j+sj][i+si].y;
1194: bc_vals[j] = bc_val;
1195: }
1196: nbcs = 0;
1197: if ((si+nx) == (M)) {
1198: nbcs = ny;
1199: }
1201: if (b) {
1202: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1203: VecAssemblyBegin(b);
1204: VecAssemblyEnd(b);
1205: }
1206: if (A) {
1207: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1208: }
1210: PetscFree(bc_vals);
1211: PetscFree(bc_global_ids);
1213: DMDAVecRestoreArray(cda,coords,&_coords);
1214: return(0);
1215: }
1219: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1220: {
1221: DM cda;
1222: Vec coords;
1223: PetscInt si,sj,nx,ny,i,j;
1224: PetscInt M,N;
1225: DMDACoor2d **_coords;
1226: PetscInt *g_idx;
1227: PetscInt *bc_global_ids;
1228: PetscScalar *bc_vals;
1229: PetscInt nbcs;
1230: PetscInt n_dofs;
1234: /* enforce bc's */
1235: DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);
1237: DMDAGetCoordinateDA(da,&cda);
1238: DMDAGetGhostedCoordinates(da,&coords);
1239: DMDAVecGetArray(cda,coords,&_coords);
1240: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1241: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1243: /* /// */
1245: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1246: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1248: /* init the entries to -1 so VecSetValues will ignore them */
1249: for (i = 0; i < ny*n_dofs; i++) {
1250: bc_global_ids[i] = -1;
1251: }
1253: i = 0;
1254: for (j = 0; j < ny; j++) {
1255: PetscInt local_id;
1256: PETSC_UNUSED PetscScalar coordx,coordy;
1258: local_id = i+j*nx;
1260: bc_global_ids[j] = g_idx[ n_dofs*local_id+d_idx ];
1262: coordx = _coords[j+sj][i+si].x;
1263: coordy = _coords[j+sj][i+si].y;
1265: bc_vals[j] = bc_val;
1266: }
1267: nbcs = 0;
1268: if (si == 0) {
1269: nbcs = ny;
1270: }
1272: if (b) {
1273: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1274: VecAssemblyBegin(b);
1275: VecAssemblyEnd(b);
1276: }
1277: if (A) {
1278: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1279: }
1281: PetscFree(bc_vals);
1282: PetscFree(bc_global_ids);
1284: DMDAVecRestoreArray(cda,coords,&_coords);
1285: return(0);
1286: }
1290: static PetscErrorCode DMDABCApplyCompression(DM elas_da,Mat A,Vec f)
1291: {
1295: BCApply_EAST(elas_da,0,-1.0,A,f);
1296: BCApply_EAST(elas_da,1, 0.0,A,f);
1297: BCApply_WEST(elas_da,0,1.0,A,f);
1298: BCApply_WEST(elas_da,1,0.0,A,f);
1299: return(0);
1300: }
1304: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff)
1305: {
1307: PetscInt start,end,m;
1308: PetscInt *unconstrained;
1309: PetscInt cnt,i;
1310: Vec x;
1311: PetscScalar *_x;
1312: IS is;
1313: VecScatter scat;
1316: /* push bc's into f and A */
1317: VecDuplicate(f,&x);
1318: BCApply_EAST(elas_da,0,-1.0,A,x);
1319: BCApply_EAST(elas_da,1, 0.0,A,x);
1320: BCApply_WEST(elas_da,0,1.0,A,x);
1321: BCApply_WEST(elas_da,1,0.0,A,x);
1323: /* define which dofs are not constrained */
1324: VecGetLocalSize(x,&m);
1325: PetscMalloc(sizeof(PetscInt)*m,&unconstrained);
1326: VecGetOwnershipRange(x,&start,&end);
1327: VecGetArray(x,&_x);
1328: cnt = 0;
1329: for (i = 0; i < m; i++ ) {
1330: PetscReal val;
1332: val = PetscRealPart(_x[i]);
1333: if( fabs(val) < 0.1 ) {
1334: unconstrained[cnt] = start + i;
1335: cnt++;
1336: }
1337: }
1338: VecRestoreArray(x,&_x);
1340: ISCreateGeneral(PETSC_COMM_WORLD,cnt,unconstrained,PETSC_COPY_VALUES,&is);
1341: PetscFree(unconstrained);
1343: /* define correction for dirichlet in the rhs */
1344: MatMult(A,x,f);
1345: VecScale(f,-1.0);
1347: /* get new matrix */
1348: MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,AA);
1349: /* get new vector */
1350: MatGetVecs(*AA,PETSC_NULL,ff);
1352: VecScatterCreate(f,is,*ff,PETSC_NULL,&scat);
1353: VecScatterBegin(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1354: VecScatterEnd(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1355: VecScatterDestroy(&scat);
1357: *dofs = is;
1358: VecDestroy(&x);
1359: return(0);
1360: }