Actual source code: ex55.c
1: static char help[] =
2: "ex55: 2D, bi-linear quadrilateral (Q1), displacement finite element formulation\n\
3: of plain strain linear elasticity, that uses the GAMG PC. E=1.0, nu=0.25.\n\
4: Unit square domain with Dirichelet boundary condition on the y=0 side only.\n\
5: Load of 1.0 in x direction on all nodes (not a true uniform load).\n\
6: -ne <size> : number of (square) quadrilateral elements in each dimention\n\
7: -alpha <v> : scaling of material coeficient in embedded circle\n\n";
9: #include <petscksp.h>
13: int main(int argc,char **args)
14: {
15: Mat Amat,Pmat;
17: PetscInt i,m,M,its,Istart,Iend,j,Ii,ix,ne=4;
18: PetscReal x,y,h;
19: Vec xx,bb;
20: KSP ksp;
21: PetscReal soft_alpha = 1.e-3;
22: MPI_Comm wcomm;
23: PetscMPIInt npe,mype;
24: PC pc;
25: PetscScalar DD[8][8],DD2[8][8];
26: #if defined(PETSC_USE_LOG)
27: PetscLogStage stage[2];
28: #endif
29: PetscScalar DD1[8][8] = { {5.333333333333333E-01, 2.0000E-01, -3.333333333333333E-01, 0.0000E+00, -2.666666666666667E-01, -2.0000E-01, 6.666666666666667E-02, 0.0000E-00 },
30: {2.0000E-01, 5.333333333333333E-01, 0.0000E-00, 6.666666666666667E-02, -2.0000E-01, -2.666666666666667E-01, 0.0000E-00, -3.333333333333333E-01 },
31: {-3.333333333333333E-01, 0.0000E-00, 5.333333333333333E-01, -2.0000E-01, 6.666666666666667E-02, 0.0000E-00, -2.666666666666667E-01, 2.0000E-01 },
32: {0.0000E+00, 6.666666666666667E-02, -2.0000E-01, 5.333333333333333E-01, 0.0000E-00, -3.333333333333333E-01, 2.0000E-01, -2.666666666666667E-01 },
33: {-2.666666666666667E-01, -2.0000E-01, 6.666666666666667E-02, 0.0000E-00, 5.333333333333333E-01, 2.0000E-01, -3.333333333333333E-01, 0.0000E+00 },
34: {-2.0000E-01, -2.666666666666667E-01, 0.0000E-00, -3.333333333333333E-01, 2.0000E-01, 5.333333333333333E-01, 0.0000E-00, 6.666666666666667E-02 },
35: {6.666666666666667E-02, 0.0000E-00, -2.666666666666667E-01, 2.0000E-01, -3.333333333333333E-01, 0.0000E-00, 5.333333333333333E-01, -2.0000E-01 },
36: {0.0000E-00, -3.333333333333333E-01, 2.0000E-01, -2.666666666666667E-01, 0.0000E-00, 6.666666666666667E-02, -2.0000E-01, 5.333333333333333E-01 } };
39: PetscInitialize(&argc,&args,(char *)0,help);
40: wcomm = PETSC_COMM_WORLD;
41: MPI_Comm_rank( wcomm, &mype );
42: MPI_Comm_size( wcomm, &npe );
43: PetscOptionsGetInt(PETSC_NULL,"-ne",&ne,PETSC_NULL);
44: h = 1./ne;
45: /* ne*ne; number of global elements */
46: PetscOptionsGetReal(PETSC_NULL,"-alpha",&soft_alpha,PETSC_NULL);
47: M = 2*(ne+1)*(ne+1); /* global number of equations */
48: m = (ne+1)*(ne+1)/npe;
49: if(mype==npe-1) m = (ne+1)*(ne+1) - (npe-1)*m;
50: m *= 2;
51: /* create stiffness matrix */
52: MatCreateMPIAIJ(wcomm,m,m,M,M,18,PETSC_NULL,6,PETSC_NULL,&Amat);
53: MatCreateMPIAIJ(wcomm,m,m,M,M,18,PETSC_NULL,6,PETSC_NULL,&Pmat);
54: MatGetOwnershipRange(Amat,&Istart,&Iend);
55: MatSetBlockSize(Amat,2);
56: MatSetBlockSize(Pmat,2);
57: m = Iend - Istart;
58: /* Generate vectors */
59: VecCreate(wcomm,&xx);
60: VecSetSizes(xx,m,M);
61: VecSetFromOptions(xx);
62: VecDuplicate(xx,&bb);
63: VecSet(bb,.0);
64: /* generate element matrices */
65: {
66: FILE *file;
67: char fname[] = "elem_2d_pln_strn_v_25.txt";
68: file = fopen(fname, "r");
69: if (file == 0) {
70: DD[0][0] = 0.53333333333333321 ;
71: DD[0][1] = 0.20000000000000001 ;
72: DD[0][2] = -0.33333333333333331 ;
73: DD[0][3] = 0.0000000000000000 ;
74: DD[0][4] = -0.26666666666666666 ;
75: DD[0][5] = -0.20000000000000001 ;
76: DD[0][6] = 6.66666666666666796E-002;
77: DD[0][7] = 6.93889390390722838E-018;
78: DD[1][0] = 0.20000000000000001 ;
79: DD[1][1] = 0.53333333333333333 ;
80: DD[1][2] = 7.80625564189563192E-018;
81: DD[1][3] = 6.66666666666666935E-002;
82: DD[1][4] = -0.20000000000000001 ;
83: DD[1][5] = -0.26666666666666666 ;
84: DD[1][6] = -3.46944695195361419E-018;
85: DD[1][7] = -0.33333333333333331 ;
86: DD[2][0] = -0.33333333333333331 ;
87: DD[2][1] = 1.12757025938492461E-017;
88: DD[2][2] = 0.53333333333333333 ;
89: DD[2][3] = -0.20000000000000001 ;
90: DD[2][4] = 6.66666666666666935E-002;
91: DD[2][5] = -6.93889390390722838E-018;
92: DD[2][6] = -0.26666666666666666 ;
93: DD[2][7] = 0.19999999999999998 ;
94: DD[3][0] = 0.0000000000000000 ;
95: DD[3][1] = 6.66666666666666935E-002;
96: DD[3][2] = -0.20000000000000001 ;
97: DD[3][3] = 0.53333333333333333 ;
98: DD[3][4] = 4.33680868994201774E-018;
99: DD[3][5] = -0.33333333333333331 ;
100: DD[3][6] = 0.20000000000000001 ;
101: DD[3][7] = -0.26666666666666666 ;
102: DD[4][0] = -0.26666666666666666 ;
103: DD[4][1] = -0.20000000000000001 ;
104: DD[4][2] = 6.66666666666666935E-002;
105: DD[4][3] = 8.67361737988403547E-019;
106: DD[4][4] = 0.53333333333333333 ;
107: DD[4][5] = 0.19999999999999998 ;
108: DD[4][6] = -0.33333333333333331 ;
109: DD[4][7] = -3.46944695195361419E-018;
110: DD[5][0] = -0.20000000000000001 ;
111: DD[5][1] = -0.26666666666666666 ;
112: DD[5][2] = -1.04083408558608426E-017;
113: DD[5][3] = -0.33333333333333331 ;
114: DD[5][4] = 0.19999999999999998 ;
115: DD[5][5] = 0.53333333333333333 ;
116: DD[5][6] = 6.93889390390722838E-018;
117: DD[5][7] = 6.66666666666666519E-002;
118: DD[6][0] = 6.66666666666666796E-002;
119: DD[6][1] = -6.93889390390722838E-018;
120: DD[6][2] = -0.26666666666666666 ;
121: DD[6][3] = 0.19999999999999998 ;
122: DD[6][4] = -0.33333333333333331 ;
123: DD[6][5] = 6.93889390390722838E-018;
124: DD[6][6] = 0.53333333333333321 ;
125: DD[6][7] = -0.20000000000000001 ;
126: DD[7][0] = 6.93889390390722838E-018;
127: DD[7][1] = -0.33333333333333331 ;
128: DD[7][2] = 0.19999999999999998 ;
129: DD[7][3] = -0.26666666666666666 ;
130: DD[7][4] = 0.0000000000000000 ;
131: DD[7][5] = 6.66666666666666519E-002;
132: DD[7][6] = -0.20000000000000001 ;
133: DD[7][7] = 0.53333333333333321 ;
134: }
135: else {
136: for(i=0;i<8;i++)
137: for(j=0;j<8;j++)
138: fscanf(file, "%le", &DD1[i][j]);
139: }
140: /* BC version of element */
141: for(i=0;i<8;i++)
142: for(j=0;j<8;j++)
143: if(i<4 || j < 4)
144: if(i==j) DD2[i][j] = .1*DD1[i][j];
145: else DD2[i][j] = 0.0;
146: else DD2[i][j] = DD1[i][j];
147: }
148: {
149: PetscReal coords[2*m];
150: /* forms the element stiffness for the Laplacian and coordinates */
151: for (Ii = Istart/2, ix = 0; Ii < Iend/2; Ii++, ix++ ) {
152: j = Ii/(ne+1); i = Ii%(ne+1);
153: /* coords */
154: x = h*(Ii % (ne+1)); y = h*(Ii/(ne+1));
155: coords[2*ix] = x; coords[2*ix+1] = y;
156: if( i<ne && j<ne ) {
157: PetscInt jj,ii,idx[4] = {Ii, Ii+1, Ii + (ne+1) + 1, Ii + (ne+1)};
158: /* radius */
159: PetscReal radius = PetscSqrtScalar( (x-.5+h/2)*(x-.5+h/2) + (y-.5+h/2)*(y-.5+h/2) );
160: PetscReal alpha = 1.0;
161: if( radius < 0.25 ){
162: alpha = soft_alpha;
163: }
164: for(ii=0;ii<8;ii++)for(jj=0;jj<8;jj++) DD[ii][jj] = alpha*DD1[ii][jj];
165: MatSetValuesBlocked(Pmat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
166: if( j>0 ) {
167: MatSetValuesBlocked(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
168: }
169: else {
170: /* a BC */
171: for(ii=0;ii<8;ii++)for(jj=0;jj<8;jj++) DD[ii][jj] = alpha*DD2[ii][jj];
172: MatSetValuesBlocked(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
173: }
174: }
175: if( j>0 ) {
176: PetscScalar v = h*h;
177: PetscInt jj = 2*Ii; /* load in x direction */
178: VecSetValues(bb,1,&jj,&v,INSERT_VALUES);
179: }
180: }
181: MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
182: MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
183: MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
184: MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
185: VecAssemblyBegin(bb);
186: VecAssemblyEnd(bb);
188: /* Setup solver */
189: KSPCreate(PETSC_COMM_WORLD,&ksp);
190: KSPSetOperators( ksp, Amat, Amat, SAME_NONZERO_PATTERN );
191: KSPSetType( ksp, KSPCG );
192: KSPGetPC( ksp, &pc );
193: PCSetType( pc, PCGAMG );
194: PCSetCoordinates( pc, 2, coords );
195: KSPSetFromOptions( ksp );
196: }
198: if( !PETSC_TRUE ) {
199: PetscViewer viewer;
200: PetscViewerASCIIOpen(wcomm, "Amat.m", &viewer);
201: PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB);
202: MatView(Amat,viewer);
203: PetscViewerDestroy( &viewer );
204: }
206: /* solve */
207: #if defined(PETSC_USE_LOG)
208: PetscLogStageRegister("Setup", &stage[0]);
209: PetscLogStageRegister("Solve", &stage[1]);
210: PetscLogStagePush(stage[0]);
211: #endif
212: KSPSetUp( ksp );
213: #if defined(PETSC_USE_LOG)
214: PetscLogStagePop();
215: #endif
217: VecSet(xx,.0);
219: #if defined(PETSC_USE_LOG)
220: PetscLogStagePush(stage[1]);
221: #endif
222: KSPSolve( ksp, bb, xx );
223: #if defined(PETSC_USE_LOG)
224: PetscLogStagePop();
225: #endif
227: KSPGetIterationNumber(ksp,&its);
229: if( !PETSC_TRUE ) {
230: PetscReal norm,norm2;
231: PetscViewer viewer;
232: Vec res;
233: MPI_Comm wcomm = ((PetscObject)bb)->comm;
234:
235: VecNorm( bb, NORM_2, &norm2 );
237: VecDuplicate( xx, &res );
238: MatMult( Amat, xx, res );
239: VecAXPY( bb, -1.0, res );
240: VecDestroy( &res );
241: VecNorm( bb, NORM_2, &norm );
242: PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e\n",0,__FUNCT__,norm/norm2,norm2);
243: PetscViewerASCIIOpen(wcomm, "residual.m", &viewer);
244: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
245: VecView(bb,viewer);
246: PetscViewerDestroy( &viewer );
249: /* PetscViewerASCIIOpen(wcomm, "rhs.m", &viewer); */
250: /* PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB ); */
251: /* CHKERRQ( ierr ); */
252: /* VecView( bb,viewer ); */
253: /* PetscViewerDestroy( &viewer ); */
255: /* PetscViewerASCIIOpen(wcomm, "solution.m", &viewer); */
256: /* PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB ); */
257: /* */
258: /* VecView( xx, viewer ); */
259: /* PetscViewerDestroy( &viewer ); */
260: }
262: /* Free work space */
263: KSPDestroy(&ksp);
264: VecDestroy(&xx);
265: VecDestroy(&bb);
266: MatDestroy(&Amat);
267: MatDestroy(&Pmat);
269: PetscFinalize();
270: return 0;
271: }