Actual source code: ex55.c
petsc-3.3-p7 2013-05-11
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>
10: #include <assert.h>
14: int main(int argc,char **args)
15: {
16: Mat Amat,Pmat;
18: PetscInt i,m,M,its,Istart,Iend,j,Ii,ix,ne=4;
19: PetscReal x,y,h;
20: Vec xx,bb;
21: KSP ksp;
22: PetscReal soft_alpha = 1.e-3;
23: MPI_Comm wcomm;
24: PetscBool use_coords = PETSC_FALSE;
25: PetscMPIInt npe,mype;
26: PC pc;
27: PetscScalar DD[8][8],DD2[8][8];
28: #if defined(PETSC_USE_LOG)
29: PetscLogStage stage[2];
30: #endif
31: 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 },
32: {2.0000E-01, 5.333333333333333E-01, 0.0000E-00, 6.666666666666667E-02, -2.0000E-01, -2.666666666666667E-01, 0.0000E-00, -3.333333333333333E-01 },
33: {-3.333333333333333E-01, 0.0000E-00, 5.333333333333333E-01, -2.0000E-01, 6.666666666666667E-02, 0.0000E-00, -2.666666666666667E-01, 2.0000E-01 },
34: {0.0000E+00, 6.666666666666667E-02, -2.0000E-01, 5.333333333333333E-01, 0.0000E-00, -3.333333333333333E-01, 2.0000E-01, -2.666666666666667E-01 },
35: {-2.666666666666667E-01, -2.0000E-01, 6.666666666666667E-02, 0.0000E-00, 5.333333333333333E-01, 2.0000E-01, -3.333333333333333E-01, 0.0000E+00 },
36: {-2.0000E-01, -2.666666666666667E-01, 0.0000E-00, -3.333333333333333E-01, 2.0000E-01, 5.333333333333333E-01, 0.0000E-00, 6.666666666666667E-02 },
37: {6.666666666666667E-02, 0.0000E-00, -2.666666666666667E-01, 2.0000E-01, -3.333333333333333E-01, 0.0000E-00, 5.333333333333333E-01, -2.0000E-01 },
38: {0.0000E-00, -3.333333333333333E-01, 2.0000E-01, -2.666666666666667E-01, 0.0000E-00, 6.666666666666667E-02, -2.0000E-01, 5.333333333333333E-01 } };
41: PetscInitialize(&argc,&args,(char *)0,help);
42: wcomm = PETSC_COMM_WORLD;
43: MPI_Comm_rank( wcomm, &mype );
44: MPI_Comm_size( wcomm, &npe );
45: PetscOptionsGetInt(PETSC_NULL,"-ne",&ne,PETSC_NULL);
46: h = 1./ne;
47: /* ne*ne; number of global elements */
48: PetscOptionsGetReal(PETSC_NULL,"-alpha",&soft_alpha,PETSC_NULL);
49: PetscOptionsGetBool(PETSC_NULL,"-use_coordinates",&use_coords,PETSC_NULL);
50: M = 2*(ne+1)*(ne+1); /* global number of equations */
51: m = (ne+1)*(ne+1)/npe;
52: if(mype==npe-1) m = (ne+1)*(ne+1) - (npe-1)*m;
53: m *= 2;
54: /* create stiffness matrix */
55: MatCreate(wcomm,&Amat);
56: MatSetSizes(Amat,m,m,M,M);
57: MatSetBlockSize(Amat,2);
58: MatSetType(Amat,MATAIJ);
59: MatSeqAIJSetPreallocation(Amat,18,PETSC_NULL);
60: MatMPIAIJSetPreallocation(Amat,18,PETSC_NULL,12,PETSC_NULL);
62: MatCreate(wcomm,&Pmat);
63: MatSetSizes(Pmat,m,m,M,M);
64: MatSetBlockSize(Pmat,2);
65: MatSetType(Pmat,MATAIJ);
66: MatSeqAIJSetPreallocation(Pmat,18,PETSC_NULL);
67: MatMPIAIJSetPreallocation(Pmat,18,PETSC_NULL,12,PETSC_NULL);
69: MatGetOwnershipRange(Amat,&Istart,&Iend);
70: assert(m == Iend - Istart);
71: /* Generate vectors */
72: VecCreate(wcomm,&xx);
73: VecSetSizes(xx,m,M);
74: VecSetFromOptions(xx);
75: VecDuplicate(xx,&bb);
76: VecSet(bb,.0);
77: /* generate element matrices */
78: {
79: FILE *file;
80: char fname[] = "data/elem_2d_pln_strn_v_25.txt";
81: file = fopen(fname, "r");
82: if (file == 0) {
83: DD[0][0] = 0.53333333333333321 ;
84: DD[0][1] = 0.20000000000000001 ;
85: DD[0][2] = -0.33333333333333331 ;
86: DD[0][3] = 0.0000000000000000 ;
87: DD[0][4] = -0.26666666666666666 ;
88: DD[0][5] = -0.20000000000000001 ;
89: DD[0][6] = 6.66666666666666796E-002;
90: DD[0][7] = 6.93889390390722838E-018;
91: DD[1][0] = 0.20000000000000001 ;
92: DD[1][1] = 0.53333333333333333 ;
93: DD[1][2] = 7.80625564189563192E-018;
94: DD[1][3] = 6.66666666666666935E-002;
95: DD[1][4] = -0.20000000000000001 ;
96: DD[1][5] = -0.26666666666666666 ;
97: DD[1][6] = -3.46944695195361419E-018;
98: DD[1][7] = -0.33333333333333331 ;
99: DD[2][0] = -0.33333333333333331 ;
100: DD[2][1] = 1.12757025938492461E-017;
101: DD[2][2] = 0.53333333333333333 ;
102: DD[2][3] = -0.20000000000000001 ;
103: DD[2][4] = 6.66666666666666935E-002;
104: DD[2][5] = -6.93889390390722838E-018;
105: DD[2][6] = -0.26666666666666666 ;
106: DD[2][7] = 0.19999999999999998 ;
107: DD[3][0] = 0.0000000000000000 ;
108: DD[3][1] = 6.66666666666666935E-002;
109: DD[3][2] = -0.20000000000000001 ;
110: DD[3][3] = 0.53333333333333333 ;
111: DD[3][4] = 4.33680868994201774E-018;
112: DD[3][5] = -0.33333333333333331 ;
113: DD[3][6] = 0.20000000000000001 ;
114: DD[3][7] = -0.26666666666666666 ;
115: DD[4][0] = -0.26666666666666666 ;
116: DD[4][1] = -0.20000000000000001 ;
117: DD[4][2] = 6.66666666666666935E-002;
118: DD[4][3] = 8.67361737988403547E-019;
119: DD[4][4] = 0.53333333333333333 ;
120: DD[4][5] = 0.19999999999999998 ;
121: DD[4][6] = -0.33333333333333331 ;
122: DD[4][7] = -3.46944695195361419E-018;
123: DD[5][0] = -0.20000000000000001 ;
124: DD[5][1] = -0.26666666666666666 ;
125: DD[5][2] = -1.04083408558608426E-017;
126: DD[5][3] = -0.33333333333333331 ;
127: DD[5][4] = 0.19999999999999998 ;
128: DD[5][5] = 0.53333333333333333 ;
129: DD[5][6] = 6.93889390390722838E-018;
130: DD[5][7] = 6.66666666666666519E-002;
131: DD[6][0] = 6.66666666666666796E-002;
132: DD[6][1] = -6.93889390390722838E-018;
133: DD[6][2] = -0.26666666666666666 ;
134: DD[6][3] = 0.19999999999999998 ;
135: DD[6][4] = -0.33333333333333331 ;
136: DD[6][5] = 6.93889390390722838E-018;
137: DD[6][6] = 0.53333333333333321 ;
138: DD[6][7] = -0.20000000000000001 ;
139: DD[7][0] = 6.93889390390722838E-018;
140: DD[7][1] = -0.33333333333333331 ;
141: DD[7][2] = 0.19999999999999998 ;
142: DD[7][3] = -0.26666666666666666 ;
143: DD[7][4] = 0.0000000000000000 ;
144: DD[7][5] = 6.66666666666666519E-002;
145: DD[7][6] = -0.20000000000000001 ;
146: DD[7][7] = 0.53333333333333321 ;
147: }
148: else {
149: for(i=0;i<8;i++) {
150: for(j=0;j<8;j++) {
151: double val;
152: fscanf(file, "%le", &val);
153: DD1[i][j] = val;
154: }
155: }
156: }
157: /* BC version of element */
158: for(i=0;i<8;i++)
159: for(j=0;j<8;j++)
160: if(i<4 || j < 4)
161: if(i==j) DD2[i][j] = .1*DD1[i][j];
162: else DD2[i][j] = 0.0;
163: else DD2[i][j] = DD1[i][j];
164: }
165: {
166: PetscReal *coords;
167: PetscMalloc(m*sizeof(PetscReal),&coords);
168: /* forms the element stiffness for the Laplacian and coordinates */
169: for (Ii = Istart/2, ix = 0; Ii < Iend/2; Ii++, ix++ ) {
170: j = Ii/(ne+1); i = Ii%(ne+1);
171: /* coords */
172: x = h*(Ii % (ne+1)); y = h*(Ii/(ne+1));
173: coords[2*ix] = x; coords[2*ix+1] = y;
174: if( i<ne && j<ne ) {
175: PetscInt jj,ii,idx[4] = {Ii, Ii+1, Ii + (ne+1) + 1, Ii + (ne+1)};
176: /* radius */
177: PetscReal radius = PetscSqrtScalar( (x-.5+h/2)*(x-.5+h/2) + (y-.5+h/2)*(y-.5+h/2) );
178: PetscReal alpha = 1.0;
179: if( radius < 0.25 ){
180: alpha = soft_alpha;
181: }
182: for(ii=0;ii<8;ii++)for(jj=0;jj<8;jj++) DD[ii][jj] = alpha*DD1[ii][jj];
183: MatSetValuesBlocked(Pmat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
184: if( j>0 ) {
185: MatSetValuesBlocked(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
186: }
187: else {
188: /* a BC */
189: for(ii=0;ii<8;ii++)for(jj=0;jj<8;jj++) DD[ii][jj] = alpha*DD2[ii][jj];
190: MatSetValuesBlocked(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
191: }
192: }
193: if( j>0 ) {
194: PetscScalar v = h*h;
195: PetscInt jj = 2*Ii; /* load in x direction */
196: VecSetValues(bb,1,&jj,&v,INSERT_VALUES);
197: }
198: }
199: MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
200: MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
201: MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
202: MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
203: VecAssemblyBegin(bb);
204: VecAssemblyEnd(bb);
206: /* Setup solver */
207: KSPCreate(PETSC_COMM_WORLD,&ksp);
208: KSPSetType( ksp, KSPCG );
209: KSPGetPC( ksp, &pc );
210: PCSetType( pc, PCGAMG );
211: KSPSetFromOptions( ksp );
212:
213: /* finish KSP/PC setup */
214: KSPSetOperators( ksp, Amat, Amat, SAME_NONZERO_PATTERN );
215: if( use_coords ) {
216: PCSetCoordinates( pc, 2, m/2, coords );
217: }
218: PetscFree(coords);
219: }
221: if( !PETSC_TRUE ) {
222: PetscViewer viewer;
223: PetscViewerASCIIOpen(wcomm, "Amat.m", &viewer);
224: PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB);
225: MatView(Amat,viewer);
226: PetscViewerDestroy( &viewer );
227: }
229: /* solve */
230: #if defined(PETSC_USE_LOG)
231: PetscLogStageRegister("Setup", &stage[0]);
232: PetscLogStageRegister("Solve", &stage[1]);
233: PetscLogStagePush(stage[0]);
234: #endif
235: KSPSetUp( ksp );
236: #if defined(PETSC_USE_LOG)
237: PetscLogStagePop();
238: #endif
240: VecSet(xx,.0);
242: #if defined(PETSC_USE_LOG)
243: PetscLogStagePush(stage[1]);
244: #endif
245: KSPSolve( ksp, bb, xx );
246: #if defined(PETSC_USE_LOG)
247: PetscLogStagePop();
248: #endif
250: KSPGetIterationNumber(ksp,&its);
252: if( !PETSC_TRUE ) {
253: PetscReal norm,norm2;
254: PetscViewer viewer;
255: Vec res;
256: MPI_Comm wcomm = ((PetscObject)bb)->comm;
257:
258: VecNorm( bb, NORM_2, &norm2 );
260: VecDuplicate( xx, &res );
261: MatMult( Amat, xx, res );
262: VecAXPY( bb, -1.0, res );
263: VecDestroy( &res );
264: VecNorm( bb, NORM_2, &norm );
265: PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e\n",0,__FUNCT__,norm/norm2,norm2);
266: PetscViewerASCIIOpen(wcomm, "residual.m", &viewer);
267: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
268: VecView(bb,viewer);
269: PetscViewerDestroy( &viewer );
272: /* PetscViewerASCIIOpen(wcomm, "rhs.m", &viewer); */
273: /* PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB ); */
274: /* CHKERRQ( ierr ); */
275: /* VecView( bb,viewer ); */
276: /* PetscViewerDestroy( &viewer ); */
278: /* PetscViewerASCIIOpen(wcomm, "solution.m", &viewer); */
279: /* PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB ); */
280: /* */
281: /* VecView( xx, viewer ); */
282: /* PetscViewerDestroy( &viewer ); */
283: }
285: /* Free work space */
286: KSPDestroy(&ksp);
287: VecDestroy(&xx);
288: VecDestroy(&bb);
289: MatDestroy(&Amat);
290: MatDestroy(&Pmat);
292: PetscFinalize();
293: return 0;
294: }