Actual source code: ex56.c

  1: static char help[] =
  2: "ex56: 3D, 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: 
 13: int main(int argc,char **args)
 14: {
 15:   Mat            Amat;
 17:   PetscInt       m,nn,M,its,Istart,Iend,i,j,k,ii,jj,kk,ic,ne=4,NP,Ni0,Nj0,Nk0,Ni1,Nj1,Nk1,ipx,ipy,ipz,id;
 18:   PetscReal      x,y,z,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[24][24],DD2[24][24];
 26: #if defined(PETSC_USE_LOG)
 27:   PetscLogStage  stage[2];
 28: #endif
 29:   PetscScalar DD1[24][24];
 30:   const PCType type;

 32:   PetscInitialize(&argc,&args,(char *)0,help);
 33:   wcomm = PETSC_COMM_WORLD;
 34:   MPI_Comm_rank( wcomm, &mype );
 35:   MPI_Comm_size( wcomm, &npe );
 36: 
 37:   PetscOptionsGetInt(PETSC_NULL,"-ne",&ne,PETSC_NULL);
 38:   h = 1./ne; nn = ne+1;
 39:   /* ne*ne; number of global elements */
 40:   PetscOptionsGetReal(PETSC_NULL,"-alpha",&soft_alpha,PETSC_NULL);
 41:   M = 3*nn*nn*nn; /* global number of equations */
 42:   m = nn*nn*nn/npe;
 43:   if(mype==npe-1) m = nn*nn*nn - (npe-1)*m;
 44:   m *= 3; /* number of equations local*/

 46:   /* Setup solver */
 47:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 48:   KSPSetType( ksp, KSPCG );
 49:   KSPGetPC( ksp, &pc );
 50:   PCSetType( pc, PCGAMG );  /* default */
 51:   KSPSetFromOptions( ksp );
 52:   PCGetType( pc, &type );

 54:   /* create stiffness matrix */
 55:   if( strcmp(type, PCPROMETHEUS) == 0 ){
 56:     /* prometheus needs BAIJ */
 57:     MatCreateMPIBAIJ(wcomm,3,m,m,M,M,81,PETSC_NULL,57,PETSC_NULL,&Amat);
 58:   }
 59:   else {
 60:     MatCreateMPIAIJ(wcomm,m,m,M,M,81,PETSC_NULL,57,PETSC_NULL,&Amat);
 61:   }
 62:   MatGetOwnershipRange(Amat,&Istart,&Iend);
 63:   MatSetBlockSize(Amat,3);
 64:   m = Iend - Istart;
 65:   /* Generate vectors */
 66:   VecCreate(wcomm,&xx);
 67:   VecSetSizes(xx,m,M);
 68:   VecSetFromOptions(xx);
 69:   VecDuplicate(xx,&bb);
 70:   VecSet(bb,.0);
 71:   /* generate element matrices */
 72:   {
 73:     FILE *file;
 74:     char fname[] = "elem_3d_elast_v_25.txt";
 75:     file = fopen(fname, "r");
 76:     if (file == 0) {
 77:       PetscPrintf(PETSC_COMM_WORLD,"\t%s failed to open input file '%s'\n",__FUNCT__,fname);
 78:       for(i=0;i<24;i++){
 79:         for(j=0;j<24;j++){
 80:           if(i==j)DD1[i][j] = 1.0;
 81:           else DD1[i][j] = -.25;
 82:         }
 83:       }
 84:     }
 85:     else {
 86:       for(i=0;i<24;i++){
 87:         for(j=0;j<24;j++){
 88:           fscanf(file, "%le", &DD1[i][j]);
 89:         }
 90:       }
 91:     }
 92:     /* BC version of element */
 93:     for(i=0;i<24;i++)
 94:       for(j=0;j<24;j++)
 95:         if(i<12 || j < 12)
 96:           if(i==j) DD2[i][j] = .1*DD1[i][j];
 97:           else DD2[i][j] = 0.0;
 98:         else DD2[i][j] = DD1[i][j];
 99:   }

101:   NP = (PetscInt)(pow((double)npe,1./3.) + .5);
102:   if(npe!=NP*NP*NP)SETERRQ1(wcomm,PETSC_ERR_ARG_WRONG, "npe=%d: npe^{1/3} must be integer",npe);
103:   if(nn!=NP*(nn/NP))SETERRQ1(wcomm,PETSC_ERR_ARG_WRONG, "-ne %d: (ne+1)%(npe^{1/3}) must equal zero",ne);
104:   ipx = mype%NP;
105:   ipy = (mype%(NP*NP))/NP;
106:   ipz = mype/(NP*NP);
107:   Ni0 = ipx*(nn/NP);
108:   Nj0 = ipy*(nn/NP);
109:   Nk0 = ipz*(nn/NP);
110:   Ni1 = Ni0 + (nn/NP);
111:   Nj1 = Nj0 + (nn/NP);
112:   Nk1 = Nk0 + (nn/NP);
113: 
114:   {
115:     PetscReal *coords;
116:     const PetscInt NN = nn/NP, id0 = ipz*nn*nn*NN + ipy*nn*NN*NN + ipx*NN*NN*NN;


119:     PetscMalloc( (m+1)*sizeof(PetscReal), &coords );
120:     coords[m] = -99.0;

122:     /* forms the element stiffness for the Laplacian and coordinates */
123:     for(i=Ni0,ic=0,ii=0;i<Ni1;i++,ii++){
124:       for(j=Nj0,jj=0;j<Nj1;j++,jj++){
125:         for(k=Nk0,kk=0;k<Nk1;k++,kk++,ic++){
126:           /* coords */
127:           x = coords[3*ic] = h*(PetscReal)i;
128:           y = coords[3*ic+1] = h*(PetscReal)j;
129:           z = coords[3*ic+2] = h*(PetscReal)k;
130:           /* matrix */
131:           id = id0 + ii + NN*jj + NN*NN*kk;

133:           if( i<ne && j<ne && k<ne) {
134:             /* radius */
135:             PetscReal radius = PetscSqrtScalar((x-.5+h/2)*(x-.5+h/2)+(y-.5+h/2)*(y-.5+h/2)+
136:                                                (z-.5+h/2)*(z-.5+h/2));
137:             PetscReal alpha = 1.0;
138:             PetscInt jx,ix,idx[8] = { id, id+1, id+NN+1, id+NN,
139:                                       id        + NN*NN, id+1    + NN*NN,
140:                                       id+NN+1 + NN*NN, id+NN + NN*NN };
141:             /* correct indices */
142:             if(i==Ni1-1 && Ni1!=nn){
143:               idx[1] += NN*(NN*NN-1);
144:               idx[2] += NN*(NN*NN-1);
145:               idx[5] += NN*(NN*NN-1);
146:               idx[6] += NN*(NN*NN-1);
147:             }
148:             if(j==Nj1-1 && Nj1!=nn) {
149:               idx[2] += NN*NN*(nn-1);
150:               idx[3] += NN*NN*(nn-1);
151:               idx[6] += NN*NN*(nn-1);
152:               idx[7] += NN*NN*(nn-1);
153:             }
154:             if(k==Nk1-1 && Nk1!=nn) {
155:               idx[4] += NN*(nn*nn-NN*NN);
156:               idx[5] += NN*(nn*nn-NN*NN);
157:               idx[6] += NN*(nn*nn-NN*NN);
158:               idx[7] += NN*(nn*nn-NN*NN);
159:             }
160: 
161:             if( radius < 0.25 ){
162:               alpha = soft_alpha;
163:             }
164:             for(ix=0;ix<24;ix++)for(jx=0;jx<24;jx++) DD[ix][jx] = alpha*DD1[ix][jx];
165:             if( k>0 ) {
166:               MatSetValuesBlocked(Amat,8,idx,8,idx,(const PetscScalar*)DD,ADD_VALUES);
167: 
168:             }
169:             else {
170:               /* a BC */
171:               for(ix=0;ix<24;ix++)for(jx=0;jx<24;jx++) DD[ix][jx] = alpha*DD2[ix][jx];
172:               MatSetValuesBlocked(Amat,8,idx,8,idx,(const PetscScalar*)DD,ADD_VALUES);
173: 
174:             }
175:           }
176:           if( k>0 ) {
177:             PetscScalar v = h*h;
178:             PetscInt jx = 3*id; /* load in x direction */
179:             VecSetValues(bb,1,&jx,&v,INSERT_VALUES);
180:           }
181:         }
182:       }
183:     }
184:     MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
185:     MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
186:     VecAssemblyBegin(bb);
187:     VecAssemblyEnd(bb);

189:     KSPSetOperators( ksp, Amat, Amat, SAME_NONZERO_PATTERN );
190:     PCSetCoordinates( pc, 3, coords );
191:     PetscFree( coords );
192:   }
193: 
194:   if( !PETSC_TRUE ) {
195:     PetscViewer viewer;
196:     PetscViewerASCIIOpen(wcomm, "Amat.m", &viewer);
197:     PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB);
198:     MatView(Amat,viewer);
199:     PetscViewerDestroy( &viewer );
200:   }

202:   /* solve */
203: #if defined(PETSC_USE_LOG)
204:   PetscLogStageRegister("Setup", &stage[0]);
205:   PetscLogStageRegister("Solve", &stage[1]);
206:   PetscLogStagePush(stage[0]);
207: #endif
208:   KSPSetUp( ksp );
209: #if defined(PETSC_USE_LOG)
210:   PetscLogStagePop();
211: #endif

213:   VecSet(xx,.0);

215: #if defined(PETSC_USE_LOG)
216:   PetscLogStagePush(stage[1]);
217: #endif
218:   KSPSolve( ksp, bb, xx );
219: #if defined(PETSC_USE_LOG)
220:   PetscLogStagePop();
221: #endif

223:   KSPGetIterationNumber(ksp,&its);

225:   if( !PETSC_TRUE ) {
226:     PetscReal norm,norm2;
227:     PetscViewer viewer;
228:     Vec res;

230:     VecNorm( bb, NORM_2, &norm2 );

232:     VecDuplicate( xx, &res );
233:     MatMult( Amat, xx, res );
234:     VecAXPY( bb, -1.0, res );
235:     VecDestroy( &res );
236:     VecNorm( bb, NORM_2, &norm );
237:     PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e\n",0,__FUNCT__,norm/norm2,norm2);
238:     PetscViewerASCIIOpen(wcomm, "residual.m", &viewer);
239:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
240:     VecView(bb,viewer);
241:     PetscViewerDestroy( &viewer );


244:     /* PetscViewerASCIIOpen(wcomm, "rhs.m", &viewer);   */
245:     /* PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB ); */
246:     /* CHKERRQ( ierr ); */
247:     /* VecView( bb,viewer );            */
248:     /* PetscViewerDestroy( &viewer );   */

250:     /* PetscViewerASCIIOpen(wcomm, "solution.m", &viewer);   */
251:     /* PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB ); */
252:     /*  */
253:     /* VecView( xx, viewer );  */
254:     /* PetscViewerDestroy( &viewer );  */
255:   }

257:   /* Free work space */
258:   KSPDestroy(&ksp);
259:   VecDestroy(&xx);
260:   VecDestroy(&bb);
261:   MatDestroy(&Amat);

263:   PetscFinalize();
264:   return 0;
265: }