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