Actual source code: ex55.c
petsc-3.5.4 2015-05-23
1: static char help[] = "2D, bi-linear quadrilateral (Q1), displacement finite element formulation\n\
2: of plain strain linear elasticity. E=1.0, nu=0.25.\n\
3: Unit square domain with Dirichelet boundary condition on the y=0 side only.\n\
4: Load of 1.0 in x direction on all nodes (not a true uniform load).\n\
5: -ne <size> : number of (square) quadrilateral elements in each dimension\n\
6: -alpha <v> : scaling of material coeficient in embedded circle\n\n";
8: #include <petscksp.h>
12: int main(int argc,char **args)
13: {
14: Mat Amat,Pmat;
16: PetscInt i,m,M,its,Istart,Iend,j,Ii,ix,ne=4;
17: PetscReal x,y,h;
18: Vec xx,bb;
19: KSP ksp;
20: PetscReal soft_alpha = 1.e-3;
21: MPI_Comm comm;
22: PetscBool use_coords = PETSC_FALSE;
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: comm = PETSC_COMM_WORLD;
41: MPI_Comm_rank(comm, &mype);
42: MPI_Comm_size(comm, &npe);
43: PetscOptionsGetInt(NULL,"-ne",&ne,NULL);
44: h = 1./ne;
45: /* ne*ne; number of global elements */
46: PetscOptionsGetReal(NULL,"-alpha",&soft_alpha,NULL);
47: PetscOptionsGetBool(NULL,"-use_coordinates",&use_coords,NULL);
48: M = 2*(ne+1)*(ne+1); /* global number of equations */
49: m = (ne+1)*(ne+1)/npe;
50: if (mype==npe-1) m = (ne+1)*(ne+1) - (npe-1)*m;
51: m *= 2;
52: /* create stiffness matrix */
53: MatCreate(comm,&Amat);
54: MatSetSizes(Amat,m,m,M,M);
55: MatSetBlockSize(Amat,2);
56: MatSetType(Amat,MATAIJ);
57: MatSeqAIJSetPreallocation(Amat,18,NULL);
58: MatMPIAIJSetPreallocation(Amat,18,NULL,18,NULL);
60: MatCreate(comm,&Pmat);
61: MatSetSizes(Pmat,m,m,M,M);
62: MatSetBlockSize(Pmat,2);
63: MatSetType(Pmat,MATAIJ);
64: MatSeqAIJSetPreallocation(Pmat,18,NULL);
65: MatMPIAIJSetPreallocation(Pmat,18,NULL,12,NULL);
67: MatGetOwnershipRange(Amat,&Istart,&Iend);
68: if (m != Iend - Istart) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"m %D does not equal Iend %D - Istart %D",m,Iend,Istart);
69: /* Generate vectors */
70: VecCreate(comm,&xx);
71: VecSetSizes(xx,m,M);
72: VecSetFromOptions(xx);
73: VecDuplicate(xx,&bb);
74: VecSet(bb,.0);
75: /* generate element matrices */
76: {
77: FILE *file;
78: char fname[] = "data/elem_2d_pln_strn_v_25.txt";
79: file = fopen(fname, "r");
80: if (file == 0) {
81: DD[0][0] = 0.53333333333333321;
82: DD[0][1] = 0.20000000000000001;
83: DD[0][2] = -0.33333333333333331;
84: DD[0][3] = 0.0000000000000000;
85: DD[0][4] = -0.26666666666666666;
86: DD[0][5] = -0.20000000000000001;
87: DD[0][6] = 6.66666666666666796E-002;
88: DD[0][7] = 6.93889390390722838E-018;
89: DD[1][0] = 0.20000000000000001;
90: DD[1][1] = 0.53333333333333333;
91: DD[1][2] = 7.80625564189563192E-018;
92: DD[1][3] = 6.66666666666666935E-002;
93: DD[1][4] = -0.20000000000000001;
94: DD[1][5] = -0.26666666666666666;
95: DD[1][6] = -3.46944695195361419E-018;
96: DD[1][7] = -0.33333333333333331;
97: DD[2][0] = -0.33333333333333331;
98: DD[2][1] = 1.12757025938492461E-017;
99: DD[2][2] = 0.53333333333333333;
100: DD[2][3] = -0.20000000000000001;
101: DD[2][4] = 6.66666666666666935E-002;
102: DD[2][5] = -6.93889390390722838E-018;
103: DD[2][6] = -0.26666666666666666;
104: DD[2][7] = 0.19999999999999998;
105: DD[3][0] = 0.0000000000000000;
106: DD[3][1] = 6.66666666666666935E-002;
107: DD[3][2] = -0.20000000000000001;
108: DD[3][3] = 0.53333333333333333;
109: DD[3][4] = 4.33680868994201774E-018;
110: DD[3][5] = -0.33333333333333331;
111: DD[3][6] = 0.20000000000000001;
112: DD[3][7] = -0.26666666666666666;
113: DD[4][0] = -0.26666666666666666;
114: DD[4][1] = -0.20000000000000001;
115: DD[4][2] = 6.66666666666666935E-002;
116: DD[4][3] = 8.67361737988403547E-019;
117: DD[4][4] = 0.53333333333333333;
118: DD[4][5] = 0.19999999999999998;
119: DD[4][6] = -0.33333333333333331;
120: DD[4][7] = -3.46944695195361419E-018;
121: DD[5][0] = -0.20000000000000001;
122: DD[5][1] = -0.26666666666666666;
123: DD[5][2] = -1.04083408558608426E-017;
124: DD[5][3] = -0.33333333333333331;
125: DD[5][4] = 0.19999999999999998;
126: DD[5][5] = 0.53333333333333333;
127: DD[5][6] = 6.93889390390722838E-018;
128: DD[5][7] = 6.66666666666666519E-002;
129: DD[6][0] = 6.66666666666666796E-002;
130: DD[6][1] = -6.93889390390722838E-018;
131: DD[6][2] = -0.26666666666666666;
132: DD[6][3] = 0.19999999999999998;
133: DD[6][4] = -0.33333333333333331;
134: DD[6][5] = 6.93889390390722838E-018;
135: DD[6][6] = 0.53333333333333321;
136: DD[6][7] = -0.20000000000000001;
137: DD[7][0] = 6.93889390390722838E-018;
138: DD[7][1] = -0.33333333333333331;
139: DD[7][2] = 0.19999999999999998;
140: DD[7][3] = -0.26666666666666666;
141: DD[7][4] = 0.0000000000000000;
142: DD[7][5] = 6.66666666666666519E-002;
143: DD[7][6] = -0.20000000000000001;
144: DD[7][7] = 0.53333333333333321;
145: } else {
146: for (i=0; i<8; i++) {
147: for (j=0; j<8; j++) {
148: double val;
149: fscanf(file, "%le", &val);
150: DD1[i][j] = val;
151: }
152: }
153: }
154: /* BC version of element */
155: for (i=0; i<8; i++) {
156: for (j=0; j<8; j++) {
157: if (i<4 || j < 4) {
158: if (i==j) DD2[i][j] = .1*DD1[i][j];
159: else DD2[i][j] = 0.0;
160: } else DD2[i][j] = DD1[i][j];
161: }
162: }
163: }
164: {
165: PetscReal *coords;
166: PetscMalloc1(m,&coords);
167: /* forms the element stiffness and coordinates */
168: for (Ii = Istart/2, ix = 0; Ii < Iend/2; Ii++, ix++) {
169: j = Ii/(ne+1); i = Ii%(ne+1);
170: /* coords */
171: x = h*(Ii % (ne+1)); y = h*(Ii/(ne+1));
172: coords[2*ix] = x; coords[2*ix+1] = y;
173: if (i<ne && j<ne) {
174: PetscInt jj,ii,idx[4] = {Ii, Ii+1, Ii + (ne+1) + 1, Ii + (ne+1)};
175: /* radius */
176: PetscReal radius = PetscSqrtReal((x-.5+h/2)*(x-.5+h/2) + (y-.5+h/2)*(y-.5+h/2));
177: PetscReal alpha = 1.0;
178: if (radius < 0.25) alpha = soft_alpha;
180: for (ii=0; ii<8; ii++) {
181: for (jj=0;jj<8;jj++) DD[ii][jj] = alpha*DD1[ii][jj];
182: }
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: } else {
187: /* a BC */
188: for (ii=0; ii<8; ii++) {
189: for (jj=0;jj<8;jj++) DD[ii][jj] = alpha*DD2[ii][jj];
190: }
191: MatSetValuesBlocked(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
192: }
193: }
194: if (j>0) {
195: PetscScalar v = h*h;
196: PetscInt jj = 2*Ii; /* load in x direction */
197: VecSetValues(bb,1,&jj,&v,INSERT_VALUES);
198: }
199: }
200: MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
201: MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
202: MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
203: MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
204: VecAssemblyBegin(bb);
205: VecAssemblyEnd(bb);
207: /* Setup solver */
208: KSPCreate(PETSC_COMM_WORLD,&ksp);
209: KSPSetType(ksp, KSPCG);
210: KSPGetPC(ksp, &pc);
211: PCSetType(pc, PCGAMG);
212: KSPSetFromOptions(ksp);
214: /* finish KSP/PC setup */
215: KSPSetOperators(ksp, Amat, Amat);
216: if (use_coords) {
217: PCSetCoordinates(pc, 2, m/2, coords);
218: }
219: PetscFree(coords);
220: }
222: if (!PETSC_TRUE) {
223: PetscViewer viewer;
224: PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
225: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
226: MatView(Amat,viewer);
227: PetscViewerDestroy(&viewer);
228: }
230: /* solve */
231: #if defined(PETSC_USE_LOG)
232: PetscLogStageRegister("Setup", &stage[0]);
233: PetscLogStageRegister("Solve", &stage[1]);
234: PetscLogStagePush(stage[0]);
235: #endif
236: KSPSetUp(ksp);
237: #if defined(PETSC_USE_LOG)
238: PetscLogStagePop();
239: #endif
241: VecSet(xx,.0);
243: #if defined(PETSC_USE_LOG)
244: PetscLogStagePush(stage[1]);
245: #endif
246: KSPSolve(ksp, bb, xx);
247: #if defined(PETSC_USE_LOG)
248: PetscLogStagePop();
249: #endif
251: KSPGetIterationNumber(ksp,&its);
253: if (0) {
254: PetscReal norm,norm2;
255: PetscViewer viewer;
256: Vec res;
258: PetscObjectGetComm((PetscObject)bb,&comm);
259: VecNorm(bb, NORM_2, &norm2);
261: VecDuplicate(xx, &res);
262: MatMult(Amat, xx, res);
263: VecAXPY(bb, -1.0, res);
264: VecDestroy(&res);
265: VecNorm(bb, NORM_2, &norm);
266: PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e\n",0,__FUNCT__,norm/norm2,norm2);
267: PetscViewerASCIIOpen(comm, "residual.m", &viewer);
268: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
269: VecView(bb,viewer);
270: PetscViewerDestroy(&viewer);
273: /* PetscViewerASCIIOpen(comm, "rhs.m", &viewer); */
274: /* PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB); */
275: /* */
276: /* VecView(bb,viewer); */
277: /* PetscViewerDestroy(&viewer); */
279: /* PetscViewerASCIIOpen(comm, "solution.m", &viewer); */
280: /* PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB); */
281: /* */
282: /* VecView(xx, viewer); */
283: /* PetscViewerDestroy(&viewer); */
284: }
286: /* Free work space */
287: KSPDestroy(&ksp);
288: VecDestroy(&xx);
289: VecDestroy(&bb);
290: MatDestroy(&Amat);
291: MatDestroy(&Pmat);
293: PetscFinalize();
294: return 0;
295: }