Actual source code: ex49.c
1: static char help[] = " Solves the compressible plane strain elasticity equations in 2d on the unit domain using Q1 finite elements. \n\
2: Material properties E (Youngs modulus) and nu (Poisson ratio) may vary as a function of space. \n\
3: The model utilises boundary conditions which produce compression in the x direction. \n\
4: Options: \n"
5: "\
6: -mx : number of elements in x-direction \n\
7: -my : number of elements in y-direction \n\
8: -c_str : structure of the coefficients to use. \n"
9: "\
10: -c_str 0 => isotropic material with constant coefficients. \n\
11: Parameters: \n\
12: -iso_E : Youngs modulus \n\
13: -iso_nu : Poisson ratio \n\
14: -c_str 1 => step function in the material properties in x. \n\
15: Parameters: \n\
16: -step_E0 : Youngs modulus to the left of the step \n\
17: -step_nu0 : Poisson ratio to the left of the step \n\
18: -step_E1 : Youngs modulus to the right of the step \n\
19: -step_n1 : Poisson ratio to the right of the step \n\
20: -step_xc : x coordinate of the step \n"
21: "\
22: -c_str 2 => checkerboard material with alternating properties. \n\
23: Repeats the following pattern throughout the domain. For example with 4 materials specified, we would heve \n\
24: -------------------------\n\
25: | D | A | B | C |\n\
26: ------|-----|-----|------\n\
27: | C | D | A | B |\n\
28: ------|-----|-----|------\n\
29: | B | C | D | A |\n\
30: ------|-----|-----|------\n\
31: | A | B | C | D |\n\
32: -------------------------\n\
33: \n\
34: Parameters: \n\
35: -brick_E : a comma separated list of Young's modulii \n\
36: -brick_nu : a comma separated list of Poisson ratios \n\
37: -brick_span : the number of elements in x and y each brick will span \n\
38: -c_str 3 => sponge-like material with alternating properties. \n\
39: Repeats the following pattern throughout the domain \n"
40: "\
41: -----------------------------\n\
42: | [background] |\n\
43: | E0,nu0 |\n\
44: | ----------------- |\n\
45: | | [inclusion] | |\n\
46: | | E1,nu1 | |\n\
47: | | | |\n\
48: | | <---- w ----> | |\n\
49: | | | |\n\
50: | | | |\n\
51: | ----------------- |\n\
52: | |\n\
53: | |\n\
54: -----------------------------\n\
55: <-------- t + w + t ------->\n\
56: \n\
57: Parameters: \n\
58: -sponge_E0 : Youngs modulus of the surrounding material \n\
59: -sponge_E1 : Youngs modulus of the inclusion \n\
60: -sponge_nu0 : Poisson ratio of the surrounding material \n\
61: -sponge_nu1 : Poisson ratio of the inclusion \n\
62: -sponge_t : the number of elements defining the border around each inclusion \n\
63: -sponge_w : the number of elements in x and y each inclusion will span\n\
64: -use_gp_coords : Evaluate the Youngs modulus, Poisson ratio and the body force at the global coordinates of the quadrature points.\n\
65: By default, E, nu and the body force are evaulated at the element center and applied as a constant over the entire element.\n\
66: -use_nonsymbc : Option to use non-symmetric boundary condition imposition. This choice will use less memory.";
68: /* Contributed by Dave May */
70: #include <petscksp.h>
71: #include <petscdm.h>
72: #include <petscdmda.h>
74: static PetscErrorCode DMDABCApplyCompression(DM,Mat,Vec);
75: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff);
77: #define NSD 2 /* number of spatial dimensions */
78: #define NODES_PER_EL 4 /* nodes per element */
79: #define U_DOFS 2 /* degrees of freedom per displacement node */
80: #define GAUSS_POINTS 4
82: /* cell based evaluation */
83: typedef struct {
84: PetscScalar E,nu,fx,fy;
85: } Coefficients;
87: /* Gauss point based evaluation 8+4+4+4 = 20 */
88: typedef struct {
89: PetscScalar gp_coords[2*GAUSS_POINTS];
90: PetscScalar E[GAUSS_POINTS];
91: PetscScalar nu[GAUSS_POINTS];
92: PetscScalar fx[GAUSS_POINTS];
93: PetscScalar fy[GAUSS_POINTS];
94: } GaussPointCoefficients;
96: typedef struct {
97: PetscScalar ux_dof;
98: PetscScalar uy_dof;
99: } ElasticityDOF;
101: /*
103: D = E/((1+nu)(1-2nu)) * [ 1-nu nu 0 ]
104: [ nu 1-nu 0 ]
105: [ 0 0 0.5*(1-2nu) ]
107: B = [ d_dx 0 ]
108: [ 0 d_dy ]
109: [ d_dy d_dx ]
111: */
113: /* FEM routines */
114: /*
115: Element: Local basis function ordering
116: 1-----2
117: | |
118: | |
119: 0-----3
120: */
121: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
122: {
123: PetscScalar xi = _xi[0];
124: PetscScalar eta = _xi[1];
126: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
127: Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
128: Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
129: Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
130: }
132: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
133: {
134: PetscScalar xi = _xi[0];
135: PetscScalar eta = _xi[1];
137: GNi[0][0] = -0.25*(1.0-eta);
138: GNi[0][1] = -0.25*(1.0+eta);
139: GNi[0][2] = 0.25*(1.0+eta);
140: GNi[0][3] = 0.25*(1.0-eta);
142: GNi[1][0] = -0.25*(1.0-xi);
143: GNi[1][1] = 0.25*(1.0-xi);
144: GNi[1][2] = 0.25*(1.0+xi);
145: GNi[1][3] = -0.25*(1.0+xi);
146: }
148: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
149: {
150: PetscScalar J00,J01,J10,J11,J;
151: PetscScalar iJ00,iJ01,iJ10,iJ11;
152: PetscInt i;
154: J00 = J01 = J10 = J11 = 0.0;
155: for (i = 0; i < NODES_PER_EL; i++) {
156: PetscScalar cx = coords[2*i+0];
157: PetscScalar cy = coords[2*i+1];
159: J00 = J00+GNi[0][i]*cx; /* J_xx = dx/dxi */
160: J01 = J01+GNi[0][i]*cy; /* J_xy = dy/dxi */
161: J10 = J10+GNi[1][i]*cx; /* J_yx = dx/deta */
162: J11 = J11+GNi[1][i]*cy; /* J_yy = dy/deta */
163: }
164: J = (J00*J11)-(J01*J10);
166: iJ00 = J11/J;
167: iJ01 = -J01/J;
168: iJ10 = -J10/J;
169: iJ11 = J00/J;
171: for (i = 0; i < NODES_PER_EL; i++) {
172: GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
173: GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
174: }
176: if (det_J) *det_J = J;
177: }
179: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
180: {
181: *ngp = 4;
182: gp_xi[0][0] = -0.57735026919;gp_xi[0][1] = -0.57735026919;
183: gp_xi[1][0] = -0.57735026919;gp_xi[1][1] = 0.57735026919;
184: gp_xi[2][0] = 0.57735026919;gp_xi[2][1] = 0.57735026919;
185: gp_xi[3][0] = 0.57735026919;gp_xi[3][1] = -0.57735026919;
186: gp_weight[0] = 1.0;
187: gp_weight[1] = 1.0;
188: gp_weight[2] = 1.0;
189: gp_weight[3] = 1.0;
190: }
192: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
193: {
194: PetscMPIInt rank;
195: PetscInt proc_I,proc_J;
196: PetscInt cpu_x,cpu_y;
197: PetscInt local_mx,local_my;
198: Vec vlx,vly;
199: PetscInt *LX,*LY,i;
200: PetscScalar *_a;
201: Vec V_SEQ;
202: VecScatter ctx;
205: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
207: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
209: proc_J = rank/cpu_x;
210: proc_I = rank-cpu_x*proc_J;
212: PetscMalloc1(cpu_x,&LX);
213: PetscMalloc1(cpu_y,&LY);
215: DMDAGetElementsSizes(da,&local_mx,&local_my,NULL);
216: VecCreate(PETSC_COMM_WORLD,&vlx);
217: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
218: VecSetFromOptions(vlx);
220: VecCreate(PETSC_COMM_WORLD,&vly);
221: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
222: VecSetFromOptions(vly);
224: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
225: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
226: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
227: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
229: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
230: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
231: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
232: VecGetArray(V_SEQ,&_a);
233: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
234: VecRestoreArray(V_SEQ,&_a);
235: VecScatterDestroy(&ctx);
236: VecDestroy(&V_SEQ);
238: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
239: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
240: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
241: VecGetArray(V_SEQ,&_a);
242: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
243: VecRestoreArray(V_SEQ,&_a);
244: VecScatterDestroy(&ctx);
245: VecDestroy(&V_SEQ);
247: *_lx = LX;
248: *_ly = LY;
250: VecDestroy(&vlx);
251: VecDestroy(&vly);
252: return 0;
253: }
255: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
256: {
257: DM cda;
258: Vec coords;
259: DMDACoor2d **_coords;
260: PetscInt si,sj,nx,ny,i,j;
261: FILE *fp;
262: char fname[PETSC_MAX_PATH_LEN];
263: PetscMPIInt rank;
266: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
267: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
268: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
270: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
272: DMGetCoordinateDM(da,&cda);
273: DMGetCoordinatesLocal(da,&coords);
274: DMDAVecGetArray(cda,coords,&_coords);
275: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
276: for (j = sj; j < sj+ny-1; j++) {
277: for (i = si; i < si+nx-1; i++) {
278: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
279: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i].x),(double)PetscRealPart(_coords[j+1][i].y));
280: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i+1].x),(double)PetscRealPart(_coords[j+1][i+1].y));
281: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i+1].x),(double)PetscRealPart(_coords[j][i+1].y));
282: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
283: }
284: }
285: DMDAVecRestoreArray(cda,coords,&_coords);
287: PetscFClose(PETSC_COMM_SELF,fp);
288: return 0;
289: }
291: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
292: {
293: DM cda;
294: Vec coords,local_fields;
295: DMDACoor2d **_coords;
296: FILE *fp;
297: char fname[PETSC_MAX_PATH_LEN];
298: const char *field_name;
299: PetscMPIInt rank;
300: PetscInt si,sj,nx,ny,i,j;
301: PetscInt n_dofs,d;
302: PetscScalar *_fields;
305: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
306: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
307: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
310: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
311: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
312: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
313: for (d = 0; d < n_dofs; d++) {
314: DMDAGetFieldName(da,d,&field_name);
315: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
316: }
317: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
319: DMGetCoordinateDM(da,&cda);
320: DMGetCoordinatesLocal(da,&coords);
321: DMDAVecGetArray(cda,coords,&_coords);
322: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
324: DMCreateLocalVector(da,&local_fields);
325: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
326: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
327: VecGetArray(local_fields,&_fields);
329: for (j = sj; j < sj+ny; j++) {
330: for (i = si; i < si+nx; i++) {
331: PetscScalar coord_x,coord_y;
332: PetscScalar field_d;
334: coord_x = _coords[j][i].x;
335: coord_y = _coords[j][i].y;
337: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
338: for (d = 0; d < n_dofs; d++) {
339: field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
340: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
341: }
342: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
343: }
344: }
345: VecRestoreArray(local_fields,&_fields);
346: VecDestroy(&local_fields);
348: DMDAVecRestoreArray(cda,coords,&_coords);
350: PetscFClose(PETSC_COMM_SELF,fp);
351: return 0;
352: }
354: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
355: {
356: DM cda;
357: Vec local_fields;
358: FILE *fp;
359: char fname[PETSC_MAX_PATH_LEN];
360: const char *field_name;
361: PetscMPIInt rank;
362: PetscInt si,sj,nx,ny,i,j,p;
363: PetscInt n_dofs,d;
364: GaussPointCoefficients **_coefficients;
367: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
368: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
369: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
372: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
373: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
374: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
375: for (d = 0; d < n_dofs; d++) {
376: DMDAGetFieldName(da,d,&field_name);
377: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
378: }
379: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
381: DMGetCoordinateDM(da,&cda);
382: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
384: DMCreateLocalVector(da,&local_fields);
385: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
386: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
387: DMDAVecGetArray(da,local_fields,&_coefficients);
389: for (j = sj; j < sj+ny; j++) {
390: for (i = si; i < si+nx; i++) {
391: PetscScalar coord_x,coord_y;
393: for (p = 0; p < GAUSS_POINTS; p++) {
394: coord_x = _coefficients[j][i].gp_coords[2*p];
395: coord_y = _coefficients[j][i].gp_coords[2*p+1];
397: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
399: PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e %1.6e\n",
400: (double)PetscRealPart(_coefficients[j][i].E[p]),(double)PetscRealPart(_coefficients[j][i].nu[p]),
401: (double)PetscRealPart(_coefficients[j][i].fx[p]),(double)PetscRealPart(_coefficients[j][i].fy[p])));
402: }
403: }
404: }
405: DMDAVecRestoreArray(da,local_fields,&_coefficients);
406: VecDestroy(&local_fields);
408: PetscFClose(PETSC_COMM_SELF,fp);
409: return 0;
410: }
412: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar E[],PetscScalar nu[])
413: {
414: PetscInt ngp;
415: PetscScalar gp_xi[GAUSS_POINTS][2];
416: PetscScalar gp_weight[GAUSS_POINTS];
417: PetscInt p,i,j,k,l;
418: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
419: PetscScalar J_p;
420: PetscScalar B[3][U_DOFS*NODES_PER_EL];
421: PetscScalar prop_E,prop_nu,factor,constit_D[3][3];
423: /* define quadrature rule */
424: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
426: /* evaluate integral */
427: for (p = 0; p < ngp; p++) {
428: ConstructQ12D_GNi(gp_xi[p],GNi_p);
429: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
431: for (i = 0; i < NODES_PER_EL; i++) {
432: PetscScalar d_dx_i = GNx_p[0][i];
433: PetscScalar d_dy_i = GNx_p[1][i];
435: B[0][2*i] = d_dx_i; B[0][2*i+1] = 0.0;
436: B[1][2*i] = 0.0; B[1][2*i+1] = d_dy_i;
437: B[2][2*i] = d_dy_i; B[2][2*i+1] = d_dx_i;
438: }
440: /* form D for the quadrature point */
441: prop_E = E[p];
442: prop_nu = nu[p];
443: factor = prop_E / ((1.0+prop_nu)*(1.0-2.0*prop_nu));
444: constit_D[0][0] = 1.0-prop_nu; constit_D[0][1] = prop_nu; constit_D[0][2] = 0.0;
445: constit_D[1][0] = prop_nu; constit_D[1][1] = 1.0-prop_nu; constit_D[1][2] = 0.0;
446: constit_D[2][0] = 0.0; constit_D[2][1] = 0.0; constit_D[2][2] = 0.5*(1.0-2.0*prop_nu);
447: for (i = 0; i < 3; i++) {
448: for (j = 0; j < 3; j++) {
449: constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
450: }
451: }
453: /* form Bt tildeD B */
454: /*
455: Ke_ij = Bt_ik . D_kl . B_lj
456: = B_ki . D_kl . B_lj
457: */
458: for (i = 0; i < 8; i++) {
459: for (j = 0; j < 8; j++) {
460: for (k = 0; k < 3; k++) {
461: for (l = 0; l < 3; l++) {
462: Ke[8*i+j] = Ke[8*i+j] + B[k][i] * constit_D[k][l] * B[l][j];
463: }
464: }
465: }
466: }
468: } /* end quadrature */
469: }
471: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
472: {
473: PetscInt ngp;
474: PetscScalar gp_xi[GAUSS_POINTS][2];
475: PetscScalar gp_weight[GAUSS_POINTS];
476: PetscInt p,i;
477: PetscScalar Ni_p[NODES_PER_EL];
478: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
479: PetscScalar J_p,fac;
481: /* define quadrature rule */
482: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
484: /* evaluate integral */
485: for (p = 0; p < ngp; p++) {
486: ConstructQ12D_Ni(gp_xi[p],Ni_p);
487: ConstructQ12D_GNi(gp_xi[p],GNi_p);
488: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
489: fac = gp_weight[p]*J_p;
491: for (i = 0; i < NODES_PER_EL; i++) {
492: Fe[NSD*i] += fac*Ni_p[i]*fx[p];
493: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
494: }
495: }
496: }
498: /*
499: i,j are the element indices
500: The unknown is a vector quantity.
501: The s[].c is used to indicate the degree of freedom.
502: */
503: static PetscErrorCode DMDAGetElementEqnums_u(MatStencil s_u[],PetscInt i,PetscInt j)
504: {
506: /* displacement */
507: /* node 0 */
508: s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0; /* Ux0 */
509: s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1; /* Uy0 */
511: /* node 1 */
512: s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0; /* Ux1 */
513: s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1; /* Uy1 */
515: /* node 2 */
516: s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0; /* Ux2 */
517: s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1; /* Uy2 */
519: /* node 3 */
520: s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0; /* Ux3 */
521: s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1; /* Uy3 */
522: return 0;
523: }
525: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
526: {
528: /* get coords for the element */
529: el_coords[NSD*0+0] = _coords[ej][ei].x; el_coords[NSD*0+1] = _coords[ej][ei].y;
530: el_coords[NSD*1+0] = _coords[ej+1][ei].x; el_coords[NSD*1+1] = _coords[ej+1][ei].y;
531: el_coords[NSD*2+0] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
532: el_coords[NSD*3+0] = _coords[ej][ei+1].x; el_coords[NSD*3+1] = _coords[ej][ei+1].y;
533: return 0;
534: }
536: static PetscErrorCode AssembleA_Elasticity(Mat A,DM elas_da,DM properties_da,Vec properties)
537: {
538: DM cda;
539: Vec coords;
540: DMDACoor2d **_coords;
541: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
542: PetscInt sex,sey,mx,my;
543: PetscInt ei,ej;
544: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
545: PetscScalar el_coords[NODES_PER_EL*NSD];
546: Vec local_properties;
547: GaussPointCoefficients **props;
548: PetscScalar *prop_E,*prop_nu;
551: /* setup for coords */
552: DMGetCoordinateDM(elas_da,&cda);
553: DMGetCoordinatesLocal(elas_da,&coords);
554: DMDAVecGetArray(cda,coords,&_coords);
556: /* setup for coefficients */
557: DMCreateLocalVector(properties_da,&local_properties);
558: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
559: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
560: DMDAVecGetArray(properties_da,local_properties,&props);
562: DMDAGetElementsCorners(elas_da,&sex,&sey,0);
563: DMDAGetElementsSizes(elas_da,&mx,&my,0);
564: for (ej = sey; ej < sey+my; ej++) {
565: for (ei = sex; ei < sex+mx; ei++) {
566: /* get coords for the element */
567: GetElementCoords(_coords,ei,ej,el_coords);
569: /* get coefficients for the element */
570: prop_E = props[ej][ei].E;
571: prop_nu = props[ej][ei].nu;
573: /* initialise element stiffness matrix */
574: PetscMemzero(Ae,sizeof(Ae));
576: /* form element stiffness matrix */
577: FormStressOperatorQ1(Ae,el_coords,prop_E,prop_nu);
579: /* insert element matrix into global matrix */
580: DMDAGetElementEqnums_u(u_eqn,ei,ej);
581: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
582: }
583: }
584: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
585: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
587: DMDAVecRestoreArray(cda,coords,&_coords);
589: DMDAVecRestoreArray(properties_da,local_properties,&props);
590: VecDestroy(&local_properties);
591: return 0;
592: }
594: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F,MatStencil u_eqn[],PetscScalar Fe_u[])
595: {
596: PetscInt n;
599: for (n = 0; n < 4; n++) {
600: 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];
601: 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];
602: }
603: return 0;
604: }
606: static PetscErrorCode AssembleF_Elasticity(Vec F,DM elas_da,DM properties_da,Vec properties)
607: {
608: DM cda;
609: Vec coords;
610: DMDACoor2d **_coords;
611: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
612: PetscInt sex,sey,mx,my;
613: PetscInt ei,ej;
614: PetscScalar Fe[NODES_PER_EL*U_DOFS];
615: PetscScalar el_coords[NODES_PER_EL*NSD];
616: Vec local_properties;
617: GaussPointCoefficients **props;
618: PetscScalar *prop_fx,*prop_fy;
619: Vec local_F;
620: ElasticityDOF **ff;
623: /* setup for coords */
624: DMGetCoordinateDM(elas_da,&cda);
625: DMGetCoordinatesLocal(elas_da,&coords);
626: DMDAVecGetArray(cda,coords,&_coords);
628: /* setup for coefficients */
629: DMGetLocalVector(properties_da,&local_properties);
630: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
631: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
632: DMDAVecGetArray(properties_da,local_properties,&props);
634: /* get access to the vector */
635: DMGetLocalVector(elas_da,&local_F);
636: VecZeroEntries(local_F);
637: DMDAVecGetArray(elas_da,local_F,&ff);
639: DMDAGetElementsCorners(elas_da,&sex,&sey,0);
640: DMDAGetElementsSizes(elas_da,&mx,&my,0);
641: for (ej = sey; ej < sey+my; ej++) {
642: for (ei = sex; ei < sex+mx; ei++) {
643: /* get coords for the element */
644: GetElementCoords(_coords,ei,ej,el_coords);
646: /* get coefficients for the element */
647: prop_fx = props[ej][ei].fx;
648: prop_fy = props[ej][ei].fy;
650: /* initialise element stiffness matrix */
651: PetscMemzero(Fe,sizeof(Fe));
653: /* form element stiffness matrix */
654: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
656: /* insert element matrix into global matrix */
657: DMDAGetElementEqnums_u(u_eqn,ei,ej);
659: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,Fe);
660: }
661: }
663: DMDAVecRestoreArray(elas_da,local_F,&ff);
664: DMLocalToGlobalBegin(elas_da,local_F,ADD_VALUES,F);
665: DMLocalToGlobalEnd(elas_da,local_F,ADD_VALUES,F);
666: DMRestoreLocalVector(elas_da,&local_F);
668: DMDAVecRestoreArray(cda,coords,&_coords);
670: DMDAVecRestoreArray(properties_da,local_properties,&props);
671: DMRestoreLocalVector(properties_da,&local_properties);
672: return 0;
673: }
675: static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
676: {
677: DM elas_da,da_prop;
678: PetscInt u_dof,dof,stencil_width;
679: Mat A;
680: PetscInt mxl,myl;
681: DM prop_cda,vel_cda;
682: Vec prop_coords,vel_coords;
683: PetscInt si,sj,nx,ny,i,j,p;
684: Vec f,X;
685: PetscInt prop_dof,prop_stencil_width;
686: Vec properties,l_properties;
687: MatNullSpace matnull;
688: PetscReal dx,dy;
689: PetscInt M,N;
690: DMDACoor2d **_prop_coords,**_vel_coords;
691: GaussPointCoefficients **element_props;
692: KSP ksp_E;
693: PetscInt coefficient_structure = 0;
694: PetscInt cpu_x,cpu_y,*lx = NULL,*ly = NULL;
695: PetscBool use_gp_coords = PETSC_FALSE;
696: PetscBool use_nonsymbc = PETSC_FALSE;
697: PetscBool no_view = PETSC_FALSE;
698: PetscBool flg;
701: /* Generate the da for velocity and pressure */
702: /*
703: We use Q1 elements for the temperature.
704: FEM has a 9-point stencil (BOX) or connectivity pattern
705: Num nodes in each direction is mx+1, my+1
706: */
707: u_dof = U_DOFS; /* Vx, Vy - velocities */
708: dof = u_dof;
709: stencil_width = 1;
710: 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);
712: DMSetMatType(elas_da,MATAIJ);
713: DMSetFromOptions(elas_da);
714: DMSetUp(elas_da);
716: DMDASetFieldName(elas_da,0,"Ux");
717: DMDASetFieldName(elas_da,1,"Uy");
719: /* unit box [0,1] x [0,1] */
720: DMDASetUniformCoordinates(elas_da,0.0,1.0,0.0,1.0,0.0,1.0);
722: /* Generate element properties, we will assume all material properties are constant over the element */
723: /* local number of elements */
724: DMDAGetElementsSizes(elas_da,&mxl,&myl,NULL);
726: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
727: DMDAGetInfo(elas_da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
728: DMDAGetElementOwnershipRanges2d(elas_da,&lx,&ly);
730: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
731: prop_stencil_width = 0;
732: 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);
733: DMSetFromOptions(da_prop);
734: DMSetUp(da_prop);
736: PetscFree(lx);
737: PetscFree(ly);
739: /* define centroid positions */
740: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
741: dx = 1.0/((PetscReal)(M));
742: dy = 1.0/((PetscReal)(N));
744: 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);
746: /* define coefficients */
747: PetscOptionsGetInt(NULL,NULL,"-c_str",&coefficient_structure,NULL);
749: DMCreateGlobalVector(da_prop,&properties);
750: DMCreateLocalVector(da_prop,&l_properties);
751: DMDAVecGetArray(da_prop,l_properties,&element_props);
753: DMGetCoordinateDM(da_prop,&prop_cda);
754: DMGetCoordinatesLocal(da_prop,&prop_coords);
755: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
757: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
759: DMGetCoordinateDM(elas_da,&vel_cda);
760: DMGetCoordinatesLocal(elas_da,&vel_coords);
761: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
763: /* interpolate the coordinates */
764: for (j = sj; j < sj+ny; j++) {
765: for (i = si; i < si+nx; i++) {
766: PetscInt ngp;
767: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
768: PetscScalar el_coords[8];
770: GetElementCoords(_vel_coords,i,j,el_coords);
771: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
773: for (p = 0; p < GAUSS_POINTS; p++) {
774: PetscScalar gp_x,gp_y;
775: PetscInt n;
776: PetscScalar xi_p[2],Ni_p[4];
778: xi_p[0] = gp_xi[p][0];
779: xi_p[1] = gp_xi[p][1];
780: ConstructQ12D_Ni(xi_p,Ni_p);
782: gp_x = 0.0;
783: gp_y = 0.0;
784: for (n = 0; n < NODES_PER_EL; n++) {
785: gp_x = gp_x+Ni_p[n]*el_coords[2*n];
786: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
787: }
788: element_props[j][i].gp_coords[2*p] = gp_x;
789: element_props[j][i].gp_coords[2*p+1] = gp_y;
790: }
791: }
792: }
794: /* define the coefficients */
795: PetscOptionsGetBool(NULL,NULL,"-use_gp_coords",&use_gp_coords,&flg);
797: for (j = sj; j < sj+ny; j++) {
798: for (i = si; i < si+nx; i++) {
799: PetscScalar centroid_x = _prop_coords[j][i].x; /* centroids of cell */
800: PetscScalar centroid_y = _prop_coords[j][i].y;
801: PETSC_UNUSED PetscScalar coord_x,coord_y;
803: if (coefficient_structure == 0) { /* isotropic */
804: PetscScalar opts_E,opts_nu;
806: opts_E = 1.0;
807: opts_nu = 0.33;
808: PetscOptionsGetScalar(NULL,NULL,"-iso_E",&opts_E,&flg);
809: PetscOptionsGetScalar(NULL,NULL,"-iso_nu",&opts_nu,&flg);
811: for (p = 0; p < GAUSS_POINTS; p++) {
812: element_props[j][i].E[p] = opts_E;
813: element_props[j][i].nu[p] = opts_nu;
815: element_props[j][i].fx[p] = 0.0;
816: element_props[j][i].fy[p] = 0.0;
817: }
818: } else if (coefficient_structure == 1) { /* step */
819: PetscScalar opts_E0,opts_nu0,opts_xc;
820: PetscScalar opts_E1,opts_nu1;
822: opts_E0 = opts_E1 = 1.0;
823: opts_nu0 = opts_nu1 = 0.333;
824: opts_xc = 0.5;
825: PetscOptionsGetScalar(NULL,NULL,"-step_E0",&opts_E0,&flg);
826: PetscOptionsGetScalar(NULL,NULL,"-step_nu0",&opts_nu0,&flg);
827: PetscOptionsGetScalar(NULL,NULL,"-step_E1",&opts_E1,&flg);
828: PetscOptionsGetScalar(NULL,NULL,"-step_nu1",&opts_nu1,&flg);
829: PetscOptionsGetScalar(NULL,NULL,"-step_xc",&opts_xc,&flg);
831: for (p = 0; p < GAUSS_POINTS; p++) {
832: coord_x = centroid_x;
833: coord_y = centroid_y;
834: if (use_gp_coords) {
835: coord_x = element_props[j][i].gp_coords[2*p];
836: coord_y = element_props[j][i].gp_coords[2*p+1];
837: }
839: element_props[j][i].E[p] = opts_E0;
840: element_props[j][i].nu[p] = opts_nu0;
841: if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
842: element_props[j][i].E[p] = opts_E1;
843: element_props[j][i].nu[p] = opts_nu1;
844: }
846: element_props[j][i].fx[p] = 0.0;
847: element_props[j][i].fy[p] = 0.0;
848: }
849: } else if (coefficient_structure == 2) { /* brick */
850: PetscReal values_E[10];
851: PetscReal values_nu[10];
852: PetscInt nbricks,maxnbricks;
853: PetscInt index,span;
854: PetscInt jj;
856: flg = PETSC_FALSE;
857: maxnbricks = 10;
858: PetscOptionsGetRealArray(NULL,NULL, "-brick_E",values_E,&maxnbricks,&flg);
859: nbricks = maxnbricks;
862: flg = PETSC_FALSE;
863: maxnbricks = 10;
864: PetscOptionsGetRealArray(NULL,NULL, "-brick_nu",values_nu,&maxnbricks,&flg);
868: span = 1;
869: PetscOptionsGetInt(NULL,NULL,"-brick_span",&span,&flg);
871: /* cycle through the indices so that no two material properties are repeated in lines of x or y */
872: jj = (j/span)%nbricks;
873: index = (jj+i/span)%nbricks;
874: /*printf("j=%d: index = %d \n", j,index); */
876: for (p = 0; p < GAUSS_POINTS; p++) {
877: element_props[j][i].E[p] = values_E[index];
878: element_props[j][i].nu[p] = values_nu[index];
879: }
880: } else if (coefficient_structure == 3) { /* sponge */
881: PetscScalar opts_E0,opts_nu0;
882: PetscScalar opts_E1,opts_nu1;
883: PetscInt opts_t,opts_w;
884: PetscInt ii,jj,ci,cj;
886: opts_E0 = opts_E1 = 1.0;
887: opts_nu0 = opts_nu1 = 0.333;
888: PetscOptionsGetScalar(NULL,NULL,"-sponge_E0",&opts_E0,&flg);
889: PetscOptionsGetScalar(NULL,NULL,"-sponge_nu0",&opts_nu0,&flg);
890: PetscOptionsGetScalar(NULL,NULL,"-sponge_E1",&opts_E1,&flg);
891: PetscOptionsGetScalar(NULL,NULL,"-sponge_nu1",&opts_nu1,&flg);
893: opts_t = opts_w = 1;
894: PetscOptionsGetInt(NULL,NULL,"-sponge_t",&opts_t,&flg);
895: PetscOptionsGetInt(NULL,NULL,"-sponge_w",&opts_w,&flg);
897: ii = (i)/(opts_t+opts_w+opts_t);
898: jj = (j)/(opts_t+opts_w+opts_t);
900: ci = i - ii*(opts_t+opts_w+opts_t);
901: cj = j - jj*(opts_t+opts_w+opts_t);
903: for (p = 0; p < GAUSS_POINTS; p++) {
904: element_props[j][i].E[p] = opts_E0;
905: element_props[j][i].nu[p] = opts_nu0;
906: }
907: if ((ci >= opts_t) && (ci < opts_t+opts_w)) {
908: if ((cj >= opts_t) && (cj < opts_t+opts_w)) {
909: for (p = 0; p < GAUSS_POINTS; p++) {
910: element_props[j][i].E[p] = opts_E1;
911: element_props[j][i].nu[p] = opts_nu1;
912: }
913: }
914: }
915: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
916: }
917: }
918: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
920: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
922: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
923: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
924: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
926: PetscOptionsGetBool(NULL,NULL,"-no_view",&no_view,NULL);
927: if (!no_view) {
928: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for elasticity eqn.","properties");
929: DMDACoordViewGnuplot2d(elas_da,"mesh");
930: }
932: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
933: DMCreateMatrix(elas_da,&A);
934: DMGetCoordinates(elas_da,&vel_coords);
935: MatNullSpaceCreateRigidBody(vel_coords,&matnull);
936: MatSetNearNullSpace(A,matnull);
937: MatNullSpaceDestroy(&matnull);
938: MatCreateVecs(A,&f,&X);
940: /* assemble A11 */
941: MatZeroEntries(A);
942: VecZeroEntries(f);
944: AssembleA_Elasticity(A,elas_da,da_prop,properties);
945: /* build force vector */
946: AssembleF_Elasticity(f,elas_da,da_prop,properties);
948: KSPCreate(PETSC_COMM_WORLD,&ksp_E);
949: KSPSetOptionsPrefix(ksp_E,"elas_"); /* elasticity */
951: PetscOptionsGetBool(NULL,NULL,"-use_nonsymbc",&use_nonsymbc,&flg);
952: /* solve */
953: if (!use_nonsymbc) {
954: Mat AA;
955: Vec ff,XX;
956: IS is;
957: VecScatter scat;
959: DMDABCApplySymmetricCompression(elas_da,A,f,&is,&AA,&ff);
960: VecDuplicate(ff,&XX);
962: KSPSetOperators(ksp_E,AA,AA);
963: KSPSetFromOptions(ksp_E);
965: KSPSolve(ksp_E,ff,XX);
967: /* push XX back into X */
968: DMDABCApplyCompression(elas_da,NULL,X);
970: VecScatterCreate(XX,NULL,X,is,&scat);
971: VecScatterBegin(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
972: VecScatterEnd(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
973: VecScatterDestroy(&scat);
975: MatDestroy(&AA);
976: VecDestroy(&ff);
977: VecDestroy(&XX);
978: ISDestroy(&is);
979: } else {
980: DMDABCApplyCompression(elas_da,A,f);
982: KSPSetOperators(ksp_E,A,A);
983: KSPSetFromOptions(ksp_E);
985: KSPSolve(ksp_E,f,X);
986: }
988: if (!no_view) DMDAViewGnuplot2d(elas_da,X,"Displacement solution for elasticity eqn.","X");
989: KSPDestroy(&ksp_E);
991: VecDestroy(&X);
992: VecDestroy(&f);
993: MatDestroy(&A);
995: DMDestroy(&elas_da);
996: DMDestroy(&da_prop);
998: VecDestroy(&properties);
999: VecDestroy(&l_properties);
1000: return 0;
1001: }
1003: int main(int argc,char **args)
1004: {
1005: PetscInt mx,my;
1007: PetscInitialize(&argc,&args,(char*)0,help);
1008: mx = my = 10;
1009: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1010: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1011: solve_elasticity_2d(mx,my);
1012: PetscFinalize();
1013: return 0;
1014: }
1016: /* -------------------------- helpers for boundary conditions -------------------------------- */
1018: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1019: {
1020: DM cda;
1021: Vec coords;
1022: PetscInt si,sj,nx,ny,i,j;
1023: PetscInt M,N;
1024: DMDACoor2d **_coords;
1025: const PetscInt *g_idx;
1026: PetscInt *bc_global_ids;
1027: PetscScalar *bc_vals;
1028: PetscInt nbcs;
1029: PetscInt n_dofs;
1030: ISLocalToGlobalMapping ltogm;
1033: /* enforce bc's */
1034: DMGetLocalToGlobalMapping(da,<ogm);
1035: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1037: DMGetCoordinateDM(da,&cda);
1038: DMGetCoordinatesLocal(da,&coords);
1039: DMDAVecGetArray(cda,coords,&_coords);
1040: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1041: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1043: /* --- */
1045: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1046: PetscMalloc1(ny*n_dofs,&bc_vals);
1048: /* init the entries to -1 so VecSetValues will ignore them */
1049: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1051: i = nx-1;
1052: for (j = 0; j < ny; j++) {
1053: PetscInt local_id;
1054: PETSC_UNUSED PetscScalar coordx,coordy;
1056: local_id = i+j*nx;
1058: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1060: coordx = _coords[j+sj][i+si].x;
1061: coordy = _coords[j+sj][i+si].y;
1063: bc_vals[j] = bc_val;
1064: }
1065: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1066: nbcs = 0;
1067: if ((si+nx) == (M)) nbcs = ny;
1069: if (b) {
1070: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1071: VecAssemblyBegin(b);
1072: VecAssemblyEnd(b);
1073: }
1074: if (A) {
1075: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1076: }
1078: PetscFree(bc_vals);
1079: PetscFree(bc_global_ids);
1081: DMDAVecRestoreArray(cda,coords,&_coords);
1082: return 0;
1083: }
1085: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1086: {
1087: DM cda;
1088: Vec coords;
1089: PetscInt si,sj,nx,ny,i,j;
1090: PetscInt M,N;
1091: DMDACoor2d **_coords;
1092: const PetscInt *g_idx;
1093: PetscInt *bc_global_ids;
1094: PetscScalar *bc_vals;
1095: PetscInt nbcs;
1096: PetscInt n_dofs;
1097: ISLocalToGlobalMapping ltogm;
1100: /* enforce bc's */
1101: DMGetLocalToGlobalMapping(da,<ogm);
1102: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1104: DMGetCoordinateDM(da,&cda);
1105: DMGetCoordinatesLocal(da,&coords);
1106: DMDAVecGetArray(cda,coords,&_coords);
1107: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1108: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1110: /* --- */
1112: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1113: PetscMalloc1(ny*n_dofs,&bc_vals);
1115: /* init the entries to -1 so VecSetValues will ignore them */
1116: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1118: i = 0;
1119: for (j = 0; j < ny; j++) {
1120: PetscInt local_id;
1121: PETSC_UNUSED PetscScalar coordx,coordy;
1123: local_id = i+j*nx;
1125: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1127: coordx = _coords[j+sj][i+si].x;
1128: coordy = _coords[j+sj][i+si].y;
1130: bc_vals[j] = bc_val;
1131: }
1132: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1133: nbcs = 0;
1134: if (si == 0) nbcs = ny;
1136: if (b) {
1137: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1138: VecAssemblyBegin(b);
1139: VecAssemblyEnd(b);
1140: }
1141: if (A) {
1142: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1143: }
1145: PetscFree(bc_vals);
1146: PetscFree(bc_global_ids);
1148: DMDAVecRestoreArray(cda,coords,&_coords);
1149: return 0;
1150: }
1152: static PetscErrorCode DMDABCApplyCompression(DM elas_da,Mat A,Vec f)
1153: {
1155: BCApply_EAST(elas_da,0,-1.0,A,f);
1156: BCApply_EAST(elas_da,1, 0.0,A,f);
1157: BCApply_WEST(elas_da,0,1.0,A,f);
1158: BCApply_WEST(elas_da,1,0.0,A,f);
1159: return 0;
1160: }
1162: static PetscErrorCode Orthogonalize(PetscInt n,Vec *vecs)
1163: {
1164: PetscInt i,j;
1165: PetscScalar dot;
1167: for (i=0; i<n; i++) {
1168: VecNormalize(vecs[i],NULL);
1169: for (j=i+1; j<n; j++) {
1170: VecDot(vecs[i],vecs[j],&dot);
1171: VecAXPY(vecs[j],-dot,vecs[i]);
1172: }
1173: }
1174: return 0;
1175: }
1177: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff)
1178: {
1179: PetscInt start,end,m;
1180: PetscInt *unconstrained;
1181: PetscInt cnt,i;
1182: Vec x;
1183: PetscScalar *_x;
1184: IS is;
1185: VecScatter scat;
1188: /* push bc's into f and A */
1189: VecDuplicate(f,&x);
1190: BCApply_EAST(elas_da,0,-1.0,A,x);
1191: BCApply_EAST(elas_da,1, 0.0,A,x);
1192: BCApply_WEST(elas_da,0,1.0,A,x);
1193: BCApply_WEST(elas_da,1,0.0,A,x);
1195: /* define which dofs are not constrained */
1196: VecGetLocalSize(x,&m);
1197: PetscMalloc1(m,&unconstrained);
1198: VecGetOwnershipRange(x,&start,&end);
1199: VecGetArray(x,&_x);
1200: cnt = 0;
1201: for (i = 0; i < m; i+=2) {
1202: PetscReal val1,val2;
1204: val1 = PetscRealPart(_x[i]);
1205: val2 = PetscRealPart(_x[i+1]);
1206: if (PetscAbs(val1) < 0.1 && PetscAbs(val2) < 0.1) {
1207: unconstrained[cnt] = start + i;
1208: cnt++;
1209: unconstrained[cnt] = start + i + 1;
1210: cnt++;
1211: }
1212: }
1213: VecRestoreArray(x,&_x);
1215: ISCreateGeneral(PETSC_COMM_WORLD,cnt,unconstrained,PETSC_COPY_VALUES,&is);
1216: PetscFree(unconstrained);
1217: ISSetBlockSize(is,2);
1219: /* define correction for dirichlet in the rhs */
1220: MatMult(A,x,f);
1221: VecScale(f,-1.0);
1223: /* get new matrix */
1224: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,AA);
1225: /* get new vector */
1226: MatCreateVecs(*AA,NULL,ff);
1228: VecScatterCreate(f,is,*ff,NULL,&scat);
1229: VecScatterBegin(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1230: VecScatterEnd(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1232: { /* Constrain near-null space */
1233: PetscInt nvecs;
1234: const Vec *vecs;
1235: Vec *uvecs;
1236: PetscBool has_const;
1237: MatNullSpace mnull,unull;
1239: MatGetNearNullSpace(A,&mnull);
1240: MatNullSpaceGetVecs(mnull,&has_const,&nvecs,&vecs);
1241: VecDuplicateVecs(*ff,nvecs,&uvecs);
1242: for (i=0; i<nvecs; i++) {
1243: VecScatterBegin(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1244: VecScatterEnd(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1245: }
1246: Orthogonalize(nvecs,uvecs);
1247: MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,nvecs,uvecs,&unull);
1248: MatSetNearNullSpace(*AA,unull);
1249: MatNullSpaceDestroy(&unull);
1250: VecDestroyVecs(nvecs,&uvecs);
1251: }
1253: VecScatterDestroy(&scat);
1255: *dofs = is;
1256: VecDestroy(&x);
1257: return 0;
1258: }
1260: /*TEST
1262: build:
1263: requires: !complex !single
1265: test:
1266: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_rtol 5e-3 -elas_ksp_view
1267: output_file: output/ex49_1.out
1269: test:
1270: suffix: 2
1271: nsize: 4
1272: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type gcr -elas_pc_type asm -elas_sub_pc_type lu -elas_ksp_rtol 5e-3
1274: test:
1275: suffix: 3
1276: nsize: 4
1277: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,10,1000,100 -brick_nu 0.4,0.2,0.3,0.1 -brick_span 3 -elas_pc_type asm -elas_sub_pc_type lu -elas_ksp_rtol 5e-3
1279: test:
1280: suffix: 4
1281: nsize: 4
1282: args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type unpreconditioned -mx 40 -my 40 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_mg_levels_ksp_type chebyshev -elas_pc_type ml -elas_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -elas_mg_levels_pc_type pbjacobi -elas_mg_levels_ksp_max_it 3 -use_nonsymbc -elas_pc_ml_nullspace user
1283: requires: ml
1285: test:
1286: suffix: 5
1287: nsize: 3
1288: args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_pc_type gamg -elas_mg_levels_ksp_type chebyshev -elas_mg_levels_ksp_max_it 1 -elas_mg_levels_ksp_chebyshev_esteig 0.2,1.1 -elas_mg_levels_pc_type jacobi
1290: test:
1291: suffix: 6
1292: nsize: 4
1293: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipegcr -elas_pc_type asm -elas_sub_pc_type lu
1295: test:
1296: suffix: 7
1297: nsize: 4
1298: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipegcr -elas_pc_type asm -elas_sub_pc_type ksp -elas_sub_ksp_ksp_type cg -elas_sub_ksp_ksp_max_it 15
1300: test:
1301: suffix: 8
1302: nsize: 4
1303: args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipefgmres -elas_pc_type asm -elas_sub_pc_type ksp -elas_sub_ksp_ksp_type cg -elas_sub_ksp_ksp_max_it 15
1305: test:
1306: suffix: hypre_nullspace
1307: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
1308: args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_pc_type hypre -elas_pc_hypre_boomeramg_nodal_coarsen 6 -elas_pc_hypre_boomeramg_vec_interp_variant 3 -elas_pc_hypre_boomeramg_interp_type ext+i -elas_ksp_view
1310: test:
1311: nsize: 4
1312: suffix: bddc
1313: args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic
1315: test:
1316: nsize: 4
1317: suffix: bddc_unsym
1318: args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic -use_nonsymbc -elas_pc_bddc_symmetric 0
1320: test:
1321: nsize: 4
1322: suffix: bddc_unsym_deluxe
1323: args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic -use_nonsymbc -elas_pc_bddc_symmetric 0 -elas_pc_bddc_use_deluxe_scaling -elas_sub_schurs_symmetric 0
1325: test:
1326: nsize: 4
1327: suffix: fetidp_unsym_deluxe
1328: args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type fetidp -elas_fetidp_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_fetidp_bddc_pc_bddc_monolithic -use_nonsymbc -elas_fetidp_bddc_pc_bddc_use_deluxe_scaling -elas_fetidp_bddc_sub_schurs_symmetric 0 -elas_fetidp_bddc_pc_bddc_deluxe_singlemat
1330: test:
1331: nsize: 4
1332: suffix: bddc_layerjump
1333: args: -mx 40 -my 40 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_ksp_norm_type natural
1335: test:
1336: nsize: 4
1337: suffix: bddc_subdomainjump
1338: args: -mx 40 -my 40 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,1000 -brick_nu 0.4,0.2 -brick_span 20 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_pc_is_use_stiffness_scaling -elas_ksp_norm_type natural
1340: test:
1341: nsize: 9
1342: suffix: bddc_subdomainjump_deluxe
1343: args: -mx 30 -my 30 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,1000 -brick_nu 0.4,0.2 -brick_span 10 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_pc_bddc_use_deluxe_scaling -elas_ksp_norm_type natural -elas_pc_bddc_schur_layers 1
1344: TEST*/