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: }