Actual source code: ex49.c
petsc-3.8.4 2018-03-24
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 : indicates the structure of the coefficients to use. \n"
9: "\
10: -c_str 0 => Setup for an 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 => Setup for a 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 => Setup for a 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 => Setup for a 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: }
196: /* 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) {
208: *mxl = m;
209: if ((sx+m) == M) *mxl = m-1; /* last proc */
210: }
211: if (myl) {
212: *myl = n;
213: if ((sy+n) == N) *myl = n-1; /* last proc */
214: }
215: if (mzl) {
216: *mzl = p;
217: if ((sz+p) == P) *mzl = p-1; /* last proc */
218: }
219: return(0);
220: }
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 (sz) {
240: *sz = sk;
241: if (sk != 0) *sz = sk+1;
242: }
244: DMDAGetLocalElementSize(da,mx,my,mz);
245: return(0);
246: }
248: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
249: {
251: PetscMPIInt rank;
252: PetscInt proc_I,proc_J;
253: PetscInt cpu_x,cpu_y;
254: PetscInt local_mx,local_my;
255: Vec vlx,vly;
256: PetscInt *LX,*LY,i;
257: PetscScalar *_a;
258: Vec V_SEQ;
259: VecScatter ctx;
262: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
264: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
266: proc_J = rank/cpu_x;
267: proc_I = rank-cpu_x*proc_J;
269: PetscMalloc1(cpu_x,&LX);
270: PetscMalloc1(cpu_y,&LY);
272: DMDAGetLocalElementSize(da,&local_mx,&local_my,NULL);
273: VecCreate(PETSC_COMM_WORLD,&vlx);
274: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
275: VecSetFromOptions(vlx);
277: VecCreate(PETSC_COMM_WORLD,&vly);
278: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
279: VecSetFromOptions(vly);
281: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
282: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
283: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
284: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
288: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
289: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
290: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
291: VecGetArray(V_SEQ,&_a);
292: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
293: VecRestoreArray(V_SEQ,&_a);
294: VecScatterDestroy(&ctx);
295: VecDestroy(&V_SEQ);
297: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
298: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
299: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
300: VecGetArray(V_SEQ,&_a);
301: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
302: VecRestoreArray(V_SEQ,&_a);
303: VecScatterDestroy(&ctx);
304: VecDestroy(&V_SEQ);
306: *_lx = LX;
307: *_ly = LY;
309: VecDestroy(&vlx);
310: VecDestroy(&vly);
311: return(0);
312: }
314: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
315: {
316: DM cda;
317: Vec coords;
318: DMDACoor2d **_coords;
319: PetscInt si,sj,nx,ny,i,j;
320: FILE *fp;
321: char fname[PETSC_MAX_PATH_LEN];
322: PetscMPIInt rank;
326: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
327: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
328: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
329: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
330: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
332: DMGetCoordinateDM(da,&cda);
333: DMGetCoordinatesLocal(da,&coords);
334: DMDAVecGetArray(cda,coords,&_coords);
335: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
336: for (j = sj; j < sj+ny-1; j++) {
337: for (i = si; i < si+nx-1; i++) {
338: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
339: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i].x),PetscRealPart(_coords[j+1][i].y));
340: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i+1].x),PetscRealPart(_coords[j+1][i+1].y));
341: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i+1].x),PetscRealPart(_coords[j][i+1].y));
342: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
343: }
344: }
345: DMDAVecRestoreArray(cda,coords,&_coords);
347: PetscFClose(PETSC_COMM_SELF,fp);
348: return(0);
349: }
351: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
352: {
353: DM cda;
354: Vec coords,local_fields;
355: DMDACoor2d **_coords;
356: FILE *fp;
357: char fname[PETSC_MAX_PATH_LEN];
358: const char *field_name;
359: PetscMPIInt rank;
360: PetscInt si,sj,nx,ny,i,j;
361: PetscInt n_dofs,d;
362: PetscScalar *_fields;
366: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
367: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
368: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
369: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
371: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
372: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
373: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
374: for (d = 0; d < n_dofs; d++) {
375: DMDAGetFieldName(da,d,&field_name);
376: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
377: }
378: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
381: DMGetCoordinateDM(da,&cda);
382: DMGetCoordinatesLocal(da,&coords);
383: DMDAVecGetArray(cda,coords,&_coords);
384: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
386: DMCreateLocalVector(da,&local_fields);
387: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
388: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
389: VecGetArray(local_fields,&_fields);
391: for (j = sj; j < sj+ny; j++) {
392: for (i = si; i < si+nx; i++) {
393: PetscScalar coord_x,coord_y;
394: PetscScalar field_d;
396: coord_x = _coords[j][i].x;
397: coord_y = _coords[j][i].y;
399: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
400: for (d = 0; d < n_dofs; d++) {
401: field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
402: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",PetscRealPart(field_d));
403: }
404: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
405: }
406: }
407: VecRestoreArray(local_fields,&_fields);
408: VecDestroy(&local_fields);
410: DMDAVecRestoreArray(cda,coords,&_coords);
412: PetscFClose(PETSC_COMM_SELF,fp);
413: return(0);
414: }
416: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
417: {
418: DM cda;
419: Vec local_fields;
420: FILE *fp;
421: char fname[PETSC_MAX_PATH_LEN];
422: const char *field_name;
423: PetscMPIInt rank;
424: PetscInt si,sj,nx,ny,i,j,p;
425: PetscInt n_dofs,d;
426: GaussPointCoefficients **_coefficients;
427: PetscErrorCode ierr;
430: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
431: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
432: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
433: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
435: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
436: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
437: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
438: for (d = 0; d < n_dofs; d++) {
439: DMDAGetFieldName(da,d,&field_name);
440: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
441: }
442: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
445: DMGetCoordinateDM(da,&cda);
446: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
448: DMCreateLocalVector(da,&local_fields);
449: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
450: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
451: DMDAVecGetArray(da,local_fields,&_coefficients);
453: for (j = sj; j < sj+ny; j++) {
454: for (i = si; i < si+nx; i++) {
455: PetscScalar coord_x,coord_y;
457: for (p = 0; p < GAUSS_POINTS; p++) {
458: coord_x = _coefficients[j][i].gp_coords[2*p];
459: coord_y = _coefficients[j][i].gp_coords[2*p+1];
461: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
463: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e %1.6e",
464: PetscRealPart(_coefficients[j][i].E[p]),PetscRealPart(_coefficients[j][i].nu[p]),
465: PetscRealPart(_coefficients[j][i].fx[p]),PetscRealPart(_coefficients[j][i].fy[p]));
466: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
467: }
468: }
469: }
470: DMDAVecRestoreArray(da,local_fields,&_coefficients);
471: VecDestroy(&local_fields);
473: PetscFClose(PETSC_COMM_SELF,fp);
474: return(0);
475: }
477: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar E[],PetscScalar nu[])
478: {
479: PetscInt ngp;
480: PetscScalar gp_xi[GAUSS_POINTS][2];
481: PetscScalar gp_weight[GAUSS_POINTS];
482: PetscInt p,i,j,k,l;
483: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
484: PetscScalar J_p;
485: PetscScalar B[3][U_DOFS*NODES_PER_EL];
486: PetscScalar prop_E,prop_nu,factor,constit_D[3][3];
488: /* define quadrature rule */
489: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
491: /* evaluate integral */
492: for (p = 0; p < ngp; p++) {
493: ConstructQ12D_GNi(gp_xi[p],GNi_p);
494: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
496: for (i = 0; i < NODES_PER_EL; i++) {
497: PetscScalar d_dx_i = GNx_p[0][i];
498: PetscScalar d_dy_i = GNx_p[1][i];
500: B[0][2*i] = d_dx_i; B[0][2*i+1] = 0.0;
501: B[1][2*i] = 0.0; B[1][2*i+1] = d_dy_i;
502: B[2][2*i] = d_dy_i; B[2][2*i+1] = d_dx_i;
503: }
505: /* form D for the quadrature point */
506: prop_E = E[p];
507: prop_nu = nu[p];
508: factor = prop_E / ((1.0+prop_nu)*(1.0-2.0*prop_nu));
509: constit_D[0][0] = 1.0-prop_nu; constit_D[0][1] = prop_nu; constit_D[0][2] = 0.0;
510: constit_D[1][0] = prop_nu; constit_D[1][1] = 1.0-prop_nu; constit_D[1][2] = 0.0;
511: constit_D[2][0] = 0.0; constit_D[2][1] = 0.0; constit_D[2][2] = 0.5*(1.0-2.0*prop_nu);
512: for (i = 0; i < 3; i++) {
513: for (j = 0; j < 3; j++) {
514: constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
515: }
516: }
518: /* form Bt tildeD B */
519: /*
520: Ke_ij = Bt_ik . D_kl . B_lj
521: = B_ki . D_kl . B_lj
522: */
523: for (i = 0; i < 8; i++) {
524: for (j = 0; j < 8; j++) {
525: for (k = 0; k < 3; k++) {
526: for (l = 0; l < 3; l++) {
527: Ke[8*i+j] = Ke[8*i+j] + B[k][i] * constit_D[k][l] * B[l][j];
528: }
529: }
530: }
531: }
533: } /* end quadrature */
534: }
536: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
537: {
538: PetscInt ngp;
539: PetscScalar gp_xi[GAUSS_POINTS][2];
540: PetscScalar gp_weight[GAUSS_POINTS];
541: PetscInt p,i;
542: PetscScalar Ni_p[NODES_PER_EL];
543: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
544: PetscScalar J_p,fac;
546: /* define quadrature rule */
547: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
549: /* evaluate integral */
550: for (p = 0; p < ngp; p++) {
551: ConstructQ12D_Ni(gp_xi[p],Ni_p);
552: ConstructQ12D_GNi(gp_xi[p],GNi_p);
553: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
554: fac = gp_weight[p]*J_p;
556: for (i = 0; i < NODES_PER_EL; i++) {
557: Fe[NSD*i] += fac*Ni_p[i]*fx[p];
558: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
559: }
560: }
561: }
563: /*
564: i,j are the element indices
565: The unknown is a vector quantity.
566: The s[].c is used to indicate the degree of freedom.
567: */
568: static PetscErrorCode DMDAGetElementEqnums_u(MatStencil s_u[],PetscInt i,PetscInt j)
569: {
571: /* displacement */
572: /* node 0 */
573: s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0; /* Ux0 */
574: s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1; /* Uy0 */
576: /* node 1 */
577: s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0; /* Ux1 */
578: s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1; /* Uy1 */
580: /* node 2 */
581: s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0; /* Ux2 */
582: s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1; /* Uy2 */
584: /* node 3 */
585: s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0; /* Ux3 */
586: s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1; /* Uy3 */
587: return(0);
588: }
590: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
591: {
593: /* get coords for the element */
594: el_coords[NSD*0+0] = _coords[ej][ei].x; el_coords[NSD*0+1] = _coords[ej][ei].y;
595: el_coords[NSD*1+0] = _coords[ej+1][ei].x; el_coords[NSD*1+1] = _coords[ej+1][ei].y;
596: el_coords[NSD*2+0] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
597: el_coords[NSD*3+0] = _coords[ej][ei+1].x; el_coords[NSD*3+1] = _coords[ej][ei+1].y;
598: return(0);
599: }
601: static PetscErrorCode AssembleA_Elasticity(Mat A,DM elas_da,DM properties_da,Vec properties)
602: {
603: DM cda;
604: Vec coords;
605: DMDACoor2d **_coords;
606: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
607: PetscInt sex,sey,mx,my;
608: PetscInt ei,ej;
609: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
610: PetscScalar el_coords[NODES_PER_EL*NSD];
611: Vec local_properties;
612: GaussPointCoefficients **props;
613: PetscScalar *prop_E,*prop_nu;
614: PetscErrorCode ierr;
617: /* setup for coords */
618: DMGetCoordinateDM(elas_da,&cda);
619: DMGetCoordinatesLocal(elas_da,&coords);
620: DMDAVecGetArray(cda,coords,&_coords);
622: /* setup for coefficients */
623: DMCreateLocalVector(properties_da,&local_properties);
624: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
625: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
626: DMDAVecGetArray(properties_da,local_properties,&props);
628: DMDAGetElementCorners(elas_da,&sex,&sey,0,&mx,&my,0);
629: for (ej = sey; ej < sey+my; ej++) {
630: for (ei = sex; ei < sex+mx; ei++) {
631: /* get coords for the element */
632: GetElementCoords(_coords,ei,ej,el_coords);
634: /* get coefficients for the element */
635: prop_E = props[ej][ei].E;
636: prop_nu = props[ej][ei].nu;
638: /* initialise element stiffness matrix */
639: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
641: /* form element stiffness matrix */
642: FormStressOperatorQ1(Ae,el_coords,prop_E,prop_nu);
644: /* insert element matrix into global matrix */
645: DMDAGetElementEqnums_u(u_eqn,ei,ej);
646: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
647: }
648: }
649: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
650: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
652: DMDAVecRestoreArray(cda,coords,&_coords);
654: DMDAVecRestoreArray(properties_da,local_properties,&props);
655: VecDestroy(&local_properties);
656: return(0);
657: }
660: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F,MatStencil u_eqn[],PetscScalar Fe_u[])
661: {
662: PetscInt n;
665: for (n = 0; n < 4; n++) {
666: 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];
667: 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];
668: }
669: return(0);
670: }
672: static PetscErrorCode AssembleF_Elasticity(Vec F,DM elas_da,DM properties_da,Vec properties)
673: {
674: DM cda;
675: Vec coords;
676: DMDACoor2d **_coords;
677: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
678: PetscInt sex,sey,mx,my;
679: PetscInt ei,ej;
680: PetscScalar Fe[NODES_PER_EL*U_DOFS];
681: PetscScalar el_coords[NODES_PER_EL*NSD];
682: Vec local_properties;
683: GaussPointCoefficients **props;
684: PetscScalar *prop_fx,*prop_fy;
685: Vec local_F;
686: ElasticityDOF **ff;
687: PetscErrorCode ierr;
690: /* setup for coords */
691: DMGetCoordinateDM(elas_da,&cda);
692: DMGetCoordinatesLocal(elas_da,&coords);
693: DMDAVecGetArray(cda,coords,&_coords);
695: /* setup for coefficients */
696: DMGetLocalVector(properties_da,&local_properties);
697: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
698: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
699: DMDAVecGetArray(properties_da,local_properties,&props);
701: /* get acces to the vector */
702: DMGetLocalVector(elas_da,&local_F);
703: VecZeroEntries(local_F);
704: DMDAVecGetArray(elas_da,local_F,&ff);
707: DMDAGetElementCorners(elas_da,&sex,&sey,0,&mx,&my,0);
708: for (ej = sey; ej < sey+my; ej++) {
709: for (ei = sex; ei < sex+mx; ei++) {
710: /* get coords for the element */
711: GetElementCoords(_coords,ei,ej,el_coords);
713: /* get coefficients for the element */
714: prop_fx = props[ej][ei].fx;
715: prop_fy = props[ej][ei].fy;
717: /* initialise element stiffness matrix */
718: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
720: /* form element stiffness matrix */
721: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
723: /* insert element matrix into global matrix */
724: DMDAGetElementEqnums_u(u_eqn,ei,ej);
726: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,Fe);
727: }
728: }
730: DMDAVecRestoreArray(elas_da,local_F,&ff);
731: DMLocalToGlobalBegin(elas_da,local_F,ADD_VALUES,F);
732: DMLocalToGlobalEnd(elas_da,local_F,ADD_VALUES,F);
733: DMRestoreLocalVector(elas_da,&local_F);
735: DMDAVecRestoreArray(cda,coords,&_coords);
737: DMDAVecRestoreArray(properties_da,local_properties,&props);
738: DMRestoreLocalVector(properties_da,&local_properties);
739: return(0);
740: }
742: static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
743: {
744: DM elas_da,da_prop;
745: PetscInt u_dof,dof,stencil_width;
746: Mat A;
747: PetscInt mxl,myl;
748: DM prop_cda,vel_cda;
749: Vec prop_coords,vel_coords;
750: PetscInt si,sj,nx,ny,i,j,p;
751: Vec f,X;
752: PetscInt prop_dof,prop_stencil_width;
753: Vec properties,l_properties;
754: MatNullSpace matnull;
755: PetscReal dx,dy;
756: PetscInt M,N;
757: DMDACoor2d **_prop_coords,**_vel_coords;
758: GaussPointCoefficients **element_props;
759: KSP ksp_E;
760: PetscInt coefficient_structure = 0;
761: PetscInt cpu_x,cpu_y,*lx = NULL,*ly = NULL;
762: PetscBool use_gp_coords = PETSC_FALSE;
763: PetscBool use_nonsymbc = PETSC_FALSE;
764: PetscBool no_view = PETSC_FALSE;
765: PetscBool flg;
766: PetscErrorCode ierr;
769: /* Generate the da for velocity and pressure */
770: /*
771: We use Q1 elements for the temperature.
772: FEM has a 9-point stencil (BOX) or connectivity pattern
773: Num nodes in each direction is mx+1, my+1
774: */
775: u_dof = U_DOFS; /* Vx, Vy - velocities */
776: dof = u_dof;
777: stencil_width = 1;
778: 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);
779: DMSetFromOptions(elas_da);
780: DMSetUp(elas_da);
782: DMDASetFieldName(elas_da,0,"Ux");
783: DMDASetFieldName(elas_da,1,"Uy");
785: /* unit box [0,1] x [0,1] */
786: DMDASetUniformCoordinates(elas_da,0.0,1.0,0.0,1.0,0.0,1.0);
789: /* Generate element properties, we will assume all material properties are constant over the element */
790: /* local number of elements */
791: DMDAGetLocalElementSize(elas_da,&mxl,&myl,NULL);
793: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
794: DMDAGetInfo(elas_da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
795: DMDAGetElementOwnershipRanges2d(elas_da,&lx,&ly);
797: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
798: prop_stencil_width = 0;
799: 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);
800: DMSetFromOptions(da_prop);
801: DMSetUp(da_prop);
803: PetscFree(lx);
804: PetscFree(ly);
806: /* define centroid positions */
807: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
808: dx = 1.0/((PetscReal)(M));
809: dy = 1.0/((PetscReal)(N));
811: 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);
813: /* define coefficients */
814: PetscOptionsGetInt(NULL,NULL,"-c_str",&coefficient_structure,NULL);
816: DMCreateGlobalVector(da_prop,&properties);
817: DMCreateLocalVector(da_prop,&l_properties);
818: DMDAVecGetArray(da_prop,l_properties,&element_props);
820: DMGetCoordinateDM(da_prop,&prop_cda);
821: DMGetCoordinatesLocal(da_prop,&prop_coords);
822: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
824: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
826: DMGetCoordinateDM(elas_da,&vel_cda);
827: DMGetCoordinatesLocal(elas_da,&vel_coords);
828: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
831: /* interpolate the coordinates */
832: for (j = sj; j < sj+ny; j++) {
833: for (i = si; i < si+nx; i++) {
834: PetscInt ngp;
835: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
836: PetscScalar el_coords[8];
838: GetElementCoords(_vel_coords,i,j,el_coords);
839: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
841: for (p = 0; p < GAUSS_POINTS; p++) {
842: PetscScalar gp_x,gp_y;
843: PetscInt n;
844: PetscScalar xi_p[2],Ni_p[4];
846: xi_p[0] = gp_xi[p][0];
847: xi_p[1] = gp_xi[p][1];
848: ConstructQ12D_Ni(xi_p,Ni_p);
850: gp_x = 0.0;
851: gp_y = 0.0;
852: for (n = 0; n < NODES_PER_EL; n++) {
853: gp_x = gp_x+Ni_p[n]*el_coords[2*n];
854: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
855: }
856: element_props[j][i].gp_coords[2*p] = gp_x;
857: element_props[j][i].gp_coords[2*p+1] = gp_y;
858: }
859: }
860: }
862: /* define the coefficients */
863: PetscOptionsGetBool(NULL,NULL,"-use_gp_coords",&use_gp_coords,&flg);
865: for (j = sj; j < sj+ny; j++) {
866: for (i = si; i < si+nx; i++) {
867: PetscScalar centroid_x = _prop_coords[j][i].x; /* centroids of cell */
868: PetscScalar centroid_y = _prop_coords[j][i].y;
869: PETSC_UNUSED PetscScalar coord_x,coord_y;
872: if (coefficient_structure == 0) { /* isotropic */
873: PetscScalar opts_E,opts_nu;
875: opts_E = 1.0;
876: opts_nu = 0.33;
877: PetscOptionsGetScalar(NULL,NULL,"-iso_E",&opts_E,&flg);
878: PetscOptionsGetScalar(NULL,NULL,"-iso_nu",&opts_nu,&flg);
880: for (p = 0; p < GAUSS_POINTS; p++) {
881: element_props[j][i].E[p] = opts_E;
882: element_props[j][i].nu[p] = opts_nu;
884: element_props[j][i].fx[p] = 0.0;
885: element_props[j][i].fy[p] = 0.0;
886: }
887: } else if (coefficient_structure == 1) { /* step */
888: PetscScalar opts_E0,opts_nu0,opts_xc;
889: PetscScalar opts_E1,opts_nu1;
891: opts_E0 = opts_E1 = 1.0;
892: opts_nu0 = opts_nu1 = 0.333;
893: opts_xc = 0.5;
894: PetscOptionsGetScalar(NULL,NULL,"-step_E0",&opts_E0,&flg);
895: PetscOptionsGetScalar(NULL,NULL,"-step_nu0",&opts_nu0,&flg);
896: PetscOptionsGetScalar(NULL,NULL,"-step_E1",&opts_E1,&flg);
897: PetscOptionsGetScalar(NULL,NULL,"-step_nu1",&opts_nu1,&flg);
898: PetscOptionsGetScalar(NULL,NULL,"-step_xc",&opts_xc,&flg);
900: for (p = 0; p < GAUSS_POINTS; p++) {
901: coord_x = centroid_x;
902: coord_y = centroid_y;
903: if (use_gp_coords) {
904: coord_x = element_props[j][i].gp_coords[2*p];
905: coord_y = element_props[j][i].gp_coords[2*p+1];
906: }
908: element_props[j][i].E[p] = opts_E0;
909: element_props[j][i].nu[p] = opts_nu0;
910: if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
911: element_props[j][i].E[p] = opts_E1;
912: element_props[j][i].nu[p] = opts_nu1;
913: }
915: element_props[j][i].fx[p] = 0.0;
916: element_props[j][i].fy[p] = 0.0;
917: }
918: } else if (coefficient_structure == 2) { /* brick */
919: PetscReal values_E[10];
920: PetscReal values_nu[10];
921: PetscInt nbricks,maxnbricks;
922: PetscInt index,span;
923: PetscInt jj;
925: flg = PETSC_FALSE;
926: maxnbricks = 10;
927: PetscOptionsGetRealArray(NULL,NULL, "-brick_E",values_E,&maxnbricks,&flg);
928: nbricks = maxnbricks;
929: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of E values for each brick");
931: flg = PETSC_FALSE;
932: maxnbricks = 10;
933: PetscOptionsGetRealArray(NULL,NULL, "-brick_nu",values_nu,&maxnbricks,&flg);
934: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of nu values for each brick");
935: if (maxnbricks != nbricks) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply equal numbers of values for E and nu");
937: span = 1;
938: PetscOptionsGetInt(NULL,NULL,"-brick_span",&span,&flg);
940: /* cycle through the indices so that no two material properties are repeated in lines of x or y */
941: jj = (j/span)%nbricks;
942: index = (jj+i/span)%nbricks;
943: /*printf("j=%d: index = %d \n", j,index); */
945: for (p = 0; p < GAUSS_POINTS; p++) {
946: element_props[j][i].E[p] = values_E[index];
947: element_props[j][i].nu[p] = values_nu[index];
948: }
949: } else if (coefficient_structure == 3) { /* sponge */
950: PetscScalar opts_E0,opts_nu0;
951: PetscScalar opts_E1,opts_nu1;
952: PetscInt opts_t,opts_w;
953: PetscInt ii,jj,ci,cj;
955: opts_E0 = opts_E1 = 1.0;
956: opts_nu0 = opts_nu1 = 0.333;
957: PetscOptionsGetScalar(NULL,NULL,"-sponge_E0",&opts_E0,&flg);
958: PetscOptionsGetScalar(NULL,NULL,"-sponge_nu0",&opts_nu0,&flg);
959: PetscOptionsGetScalar(NULL,NULL,"-sponge_E1",&opts_E1,&flg);
960: PetscOptionsGetScalar(NULL,NULL,"-sponge_nu1",&opts_nu1,&flg);
962: opts_t = opts_w = 1;
963: PetscOptionsGetInt(NULL,NULL,"-sponge_t",&opts_t,&flg);
964: PetscOptionsGetInt(NULL,NULL,"-sponge_w",&opts_w,&flg);
966: ii = (i)/(opts_t+opts_w+opts_t);
967: jj = (j)/(opts_t+opts_w+opts_t);
969: ci = i - ii*(opts_t+opts_w+opts_t);
970: cj = j - jj*(opts_t+opts_w+opts_t);
972: for (p = 0; p < GAUSS_POINTS; p++) {
973: element_props[j][i].E[p] = opts_E0;
974: element_props[j][i].nu[p] = opts_nu0;
975: }
976: if ((ci >= opts_t) && (ci < opts_t+opts_w)) {
977: if ((cj >= opts_t) && (cj < opts_t+opts_w)) {
978: for (p = 0; p < GAUSS_POINTS; p++) {
979: element_props[j][i].E[p] = opts_E1;
980: element_props[j][i].nu[p] = opts_nu1;
981: }
982: }
983: }
985: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
986: }
987: }
988: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
990: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
992: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
993: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
994: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
996: PetscOptionsGetBool(NULL,NULL,"-no_view",&no_view,NULL);
997: if (!no_view) {
998: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for elasticity eqn.","properties");
999: DMDACoordViewGnuplot2d(elas_da,"mesh");
1000: }
1002: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1003: DMSetMatType(elas_da,MATAIJ);
1004: DMCreateMatrix(elas_da,&A);
1005: MatSetBlockSize(A,2);
1006: DMGetCoordinates(elas_da,&vel_coords);
1007: MatNullSpaceCreateRigidBody(vel_coords,&matnull);
1008: MatSetNearNullSpace(A,matnull);
1009: MatNullSpaceDestroy(&matnull);
1010: MatCreateVecs(A,&f,&X);
1012: /* assemble A11 */
1013: MatZeroEntries(A);
1014: VecZeroEntries(f);
1016: AssembleA_Elasticity(A,elas_da,da_prop,properties);
1017: /* build force vector */
1018: AssembleF_Elasticity(f,elas_da,da_prop,properties);
1021: KSPCreate(PETSC_COMM_WORLD,&ksp_E);
1022: KSPSetOptionsPrefix(ksp_E,"elas_"); /* elasticity */
1024: PetscOptionsGetBool(NULL,NULL,"-use_nonsymbc",&use_nonsymbc,&flg);
1025: /* solve */
1026: if (!use_nonsymbc) {
1027: Mat AA;
1028: Vec ff,XX;
1029: IS is;
1030: VecScatter scat;
1032: DMDABCApplySymmetricCompression(elas_da,A,f,&is,&AA,&ff);
1033: VecDuplicate(ff,&XX);
1035: KSPSetOperators(ksp_E,AA,AA);
1036: KSPSetFromOptions(ksp_E);
1038: KSPSolve(ksp_E,ff,XX);
1040: /* push XX back into X */
1041: DMDABCApplyCompression(elas_da,NULL,X);
1043: VecScatterCreate(XX,NULL,X,is,&scat);
1044: VecScatterBegin(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
1045: VecScatterEnd(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
1046: VecScatterDestroy(&scat);
1048: MatDestroy(&AA);
1049: VecDestroy(&ff);
1050: VecDestroy(&XX);
1051: ISDestroy(&is);
1052: } else {
1053: DMDABCApplyCompression(elas_da,A,f);
1055: KSPSetOperators(ksp_E,A,A);
1056: KSPSetFromOptions(ksp_E);
1058: KSPSolve(ksp_E,f,X);
1059: }
1061: if (!no_view) {DMDAViewGnuplot2d(elas_da,X,"Displacement solution for elasticity eqn.","X");}
1062: KSPDestroy(&ksp_E);
1064: VecDestroy(&X);
1065: VecDestroy(&f);
1066: MatDestroy(&A);
1068: DMDestroy(&elas_da);
1069: DMDestroy(&da_prop);
1071: VecDestroy(&properties);
1072: VecDestroy(&l_properties);
1073: return(0);
1074: }
1076: int main(int argc,char **args)
1077: {
1079: PetscInt mx,my;
1081: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
1082: mx = my = 10;
1083: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1084: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1085: solve_elasticity_2d(mx,my);
1086: PetscFinalize();
1087: return ierr;
1088: }
1090: /* -------------------------- helpers for boundary conditions -------------------------------- */
1092: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1093: {
1094: DM cda;
1095: Vec coords;
1096: PetscInt si,sj,nx,ny,i,j;
1097: PetscInt M,N;
1098: DMDACoor2d **_coords;
1099: const PetscInt *g_idx;
1100: PetscInt *bc_global_ids;
1101: PetscScalar *bc_vals;
1102: PetscInt nbcs;
1103: PetscInt n_dofs;
1104: PetscErrorCode ierr;
1105: ISLocalToGlobalMapping ltogm;
1108: /* enforce bc's */
1109: DMGetLocalToGlobalMapping(da,<ogm);
1110: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1112: DMGetCoordinateDM(da,&cda);
1113: DMGetCoordinatesLocal(da,&coords);
1114: DMDAVecGetArray(cda,coords,&_coords);
1115: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1116: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1118: /* --- */
1120: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1121: PetscMalloc1(ny*n_dofs,&bc_vals);
1123: /* init the entries to -1 so VecSetValues will ignore them */
1124: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1126: i = nx-1;
1127: for (j = 0; j < ny; j++) {
1128: PetscInt local_id;
1129: PETSC_UNUSED PetscScalar coordx,coordy;
1131: local_id = i+j*nx;
1133: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1135: coordx = _coords[j+sj][i+si].x;
1136: coordy = _coords[j+sj][i+si].y;
1138: bc_vals[j] = bc_val;
1139: }
1140: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1141: nbcs = 0;
1142: if ((si+nx) == (M)) nbcs = ny;
1144: if (b) {
1145: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1146: VecAssemblyBegin(b);
1147: VecAssemblyEnd(b);
1148: }
1149: if (A) {
1150: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1151: }
1153: PetscFree(bc_vals);
1154: PetscFree(bc_global_ids);
1156: DMDAVecRestoreArray(cda,coords,&_coords);
1157: return(0);
1158: }
1160: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1161: {
1162: DM cda;
1163: Vec coords;
1164: PetscInt si,sj,nx,ny,i,j;
1165: PetscInt M,N;
1166: DMDACoor2d **_coords;
1167: const PetscInt *g_idx;
1168: PetscInt *bc_global_ids;
1169: PetscScalar *bc_vals;
1170: PetscInt nbcs;
1171: PetscInt n_dofs;
1172: PetscErrorCode ierr;
1173: ISLocalToGlobalMapping ltogm;
1176: /* enforce bc's */
1177: DMGetLocalToGlobalMapping(da,<ogm);
1178: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1180: DMGetCoordinateDM(da,&cda);
1181: DMGetCoordinatesLocal(da,&coords);
1182: DMDAVecGetArray(cda,coords,&_coords);
1183: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1184: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1186: /* --- */
1188: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1189: PetscMalloc1(ny*n_dofs,&bc_vals);
1191: /* init the entries to -1 so VecSetValues will ignore them */
1192: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1194: i = 0;
1195: for (j = 0; j < ny; j++) {
1196: PetscInt local_id;
1197: PETSC_UNUSED PetscScalar coordx,coordy;
1199: local_id = i+j*nx;
1201: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1203: coordx = _coords[j+sj][i+si].x;
1204: coordy = _coords[j+sj][i+si].y;
1206: bc_vals[j] = bc_val;
1207: }
1208: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1209: nbcs = 0;
1210: if (si == 0) nbcs = ny;
1212: if (b) {
1213: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1214: VecAssemblyBegin(b);
1215: VecAssemblyEnd(b);
1216: }
1217: if (A) {
1218: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1219: }
1221: PetscFree(bc_vals);
1222: PetscFree(bc_global_ids);
1224: DMDAVecRestoreArray(cda,coords,&_coords);
1225: return(0);
1226: }
1228: static PetscErrorCode DMDABCApplyCompression(DM elas_da,Mat A,Vec f)
1229: {
1233: BCApply_EAST(elas_da,0,-1.0,A,f);
1234: BCApply_EAST(elas_da,1, 0.0,A,f);
1235: BCApply_WEST(elas_da,0,1.0,A,f);
1236: BCApply_WEST(elas_da,1,0.0,A,f);
1237: return(0);
1238: }
1240: static PetscErrorCode Orthogonalize(PetscInt n,Vec *vecs)
1241: {
1242: PetscInt i,j;
1243: PetscScalar dot;
1247: for (i=0; i<n; i++) {
1248: VecNormalize(vecs[i],NULL);
1249: for (j=i+1; j<n; j++) {
1250: VecDot(vecs[i],vecs[j],&dot);
1251: VecAXPY(vecs[j],-dot,vecs[i]);
1252: }
1253: }
1254: return(0);
1255: }
1257: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff)
1258: {
1260: PetscInt start,end,m;
1261: PetscInt *unconstrained;
1262: PetscInt cnt,i;
1263: Vec x;
1264: PetscScalar *_x;
1265: IS is;
1266: VecScatter scat;
1269: /* push bc's into f and A */
1270: VecDuplicate(f,&x);
1271: BCApply_EAST(elas_da,0,-1.0,A,x);
1272: BCApply_EAST(elas_da,1, 0.0,A,x);
1273: BCApply_WEST(elas_da,0,1.0,A,x);
1274: BCApply_WEST(elas_da,1,0.0,A,x);
1276: /* define which dofs are not constrained */
1277: VecGetLocalSize(x,&m);
1278: PetscMalloc1(m,&unconstrained);
1279: VecGetOwnershipRange(x,&start,&end);
1280: VecGetArray(x,&_x);
1281: cnt = 0;
1282: for (i = 0; i < m; i+=2) {
1283: PetscReal val1,val2;
1285: val1 = PetscRealPart(_x[i]);
1286: val2 = PetscRealPart(_x[i+1]);
1287: if (PetscAbs(val1) < 0.1 && PetscAbs(val2) < 0.1) {
1288: unconstrained[cnt] = start + i;
1289: cnt++;
1290: unconstrained[cnt] = start + i + 1;
1291: cnt++;
1292: }
1293: }
1294: VecRestoreArray(x,&_x);
1296: ISCreateGeneral(PETSC_COMM_WORLD,cnt,unconstrained,PETSC_COPY_VALUES,&is);
1297: PetscFree(unconstrained);
1298: ISSetBlockSize(is,2);
1300: /* define correction for dirichlet in the rhs */
1301: MatMult(A,x,f);
1302: VecScale(f,-1.0);
1304: /* get new matrix */
1305: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,AA);
1306: /* get new vector */
1307: MatCreateVecs(*AA,NULL,ff);
1309: VecScatterCreate(f,is,*ff,NULL,&scat);
1310: VecScatterBegin(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1311: VecScatterEnd(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1313: { /* Constrain near-null space */
1314: PetscInt nvecs;
1315: const Vec *vecs;
1316: Vec *uvecs;
1317: PetscBool has_const;
1318: MatNullSpace mnull,unull;
1320: MatGetNearNullSpace(A,&mnull);
1321: MatNullSpaceGetVecs(mnull,&has_const,&nvecs,&vecs);
1322: VecDuplicateVecs(*ff,nvecs,&uvecs);
1323: for (i=0; i<nvecs; i++) {
1324: VecScatterBegin(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1325: VecScatterEnd(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1326: }
1327: Orthogonalize(nvecs,uvecs);
1328: MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,nvecs,uvecs,&unull);
1329: MatSetNearNullSpace(*AA,unull);
1330: MatNullSpaceDestroy(&unull);
1331: VecDestroyVecs(nvecs,&uvecs);
1332: }
1334: VecScatterDestroy(&scat);
1336: *dofs = is;
1337: VecDestroy(&x);
1338: return(0);
1339: }