Actual source code: ex10.c
2: static char help[] = "This example calculates the stiffness matrix for a brick in three\n\
3: dimensions using 20 node serendipity elements and the equations of linear\n\
4: elasticity. This also demonstrates use of block\n\
5: diagonal data structure. Input arguments are:\n\
6: -m : problem size\n\n";
8: #include petscksp.h
10: /* This code is not intended as an efficient implementation, it is only
11: here to produce an interesting sparse matrix quickly.
13: PLEASE DO NOT BASE ANY OF YOUR CODES ON CODE LIKE THIS, THERE ARE MUCH
14: BETTER WAYS TO DO THIS. */
24: int main(int argc,char **args)
25: {
26: Mat mat;
28: PetscInt i,its,m = 3,rdim,cdim,rstart,rend;
29: PetscMPIInt rank,size;
30: PetscScalar v,neg1 = -1.0;
31: Vec u,x,b;
32: KSP ksp;
33: PetscReal norm;
35: PetscInitialize(&argc,&args,(char *)0,help);
36: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
37: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
38: MPI_Comm_size(PETSC_COMM_WORLD,&size);
40: /* Form matrix */
41: GetElasticityMatrix(m,&mat);
43: /* Generate vectors */
44: MatGetSize(mat,&rdim,&cdim);
45: MatGetOwnershipRange(mat,&rstart,&rend);
46: VecCreate(PETSC_COMM_WORLD,&u);
47: VecSetSizes(u,PETSC_DECIDE,rdim);
48: VecSetFromOptions(u);
49: VecDuplicate(u,&b);
50: VecDuplicate(b,&x);
51: for (i=rstart; i<rend; i++) {
52: v = (PetscScalar)(i-rstart + 100*rank);
53: VecSetValues(u,1,&i,&v,INSERT_VALUES);
54: }
55: VecAssemblyBegin(u);
56: VecAssemblyEnd(u);
57:
58: /* Compute right-hand-side */
59: MatMult(mat,u,b);
60:
61: /* Solve linear system */
62: KSPCreate(PETSC_COMM_WORLD,&ksp);
63: KSPSetOperators(ksp,mat,mat,SAME_NONZERO_PATTERN);
64: KSPGMRESSetRestart(ksp,2*m);
65: KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
66: KSPSetType(ksp,KSPCG);
67: KSPSetFromOptions(ksp);
68: KSPSolve(ksp,b,x);
69: KSPGetIterationNumber(ksp,&its);
70: /* Check error */
71: VecAXPY(x,neg1,u);
72: VecNorm(x,NORM_2,&norm);
74: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Number of iterations %D\n",norm,its);
76: /* Free work space */
77: KSPDestroy(ksp);
78: VecDestroy(u);
79: VecDestroy(x);
80: VecDestroy(b);
81: MatDestroy(mat);
83: PetscFinalize();
84: return 0;
85: }
86: /* -------------------------------------------------------------------- */
89: /*
90: GetElasticityMatrix - Forms 3D linear elasticity matrix.
91: */
92: PetscErrorCode GetElasticityMatrix(PetscInt m,Mat *newmat)
93: {
94: PetscInt i,j,k,i1,i2,j_1,j2,k1,k2,h1,h2,shiftx,shifty,shiftz;
95: PetscInt ict,nz,base,r1,r2,N,*rowkeep,nstart;
97: IS iskeep;
98: PetscReal **K,norm;
99: Mat mat,submat = 0,*submatb;
100: const MatType type = MATSEQBAIJ;
102: m /= 2; /* This is done just to be consistent with the old example */
103: N = 3*(2*m+1)*(2*m+1)*(2*m+1);
104: PetscPrintf(PETSC_COMM_SELF,"m = %D, N=%D\n",m,N);
105: MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,80,PETSC_NULL,&mat);
107: /* Form stiffness for element */
108: PetscMalloc(81*sizeof(PetscReal *),&K);
109: for (i=0; i<81; i++) {
110: PetscMalloc(81*sizeof(PetscReal),&K[i]);
111: }
112: Elastic20Stiff(K);
114: /* Loop over elements and add contribution to stiffness */
115: shiftx = 3; shifty = 3*(2*m+1); shiftz = 3*(2*m+1)*(2*m+1);
116: for (k=0; k<m; k++) {
117: for (j=0; j<m; j++) {
118: for (i=0; i<m; i++) {
119: h1 = 0;
120: base = 2*k*shiftz + 2*j*shifty + 2*i*shiftx;
121: for (k1=0; k1<3; k1++) {
122: for (j_1=0; j_1<3; j_1++) {
123: for (i1=0; i1<3; i1++) {
124: h2 = 0;
125: r1 = base + i1*shiftx + j_1*shifty + k1*shiftz;
126: for (k2=0; k2<3; k2++) {
127: for (j2=0; j2<3; j2++) {
128: for (i2=0; i2<3; i2++) {
129: r2 = base + i2*shiftx + j2*shifty + k2*shiftz;
130: AddElement(mat,r1,r2,K,h1,h2);
131: h2 += 3;
132: }
133: }
134: }
135: h1 += 3;
136: }
137: }
138: }
139: }
140: }
141: }
143: for (i=0; i<81; i++) {
144: PetscFree(K[i]);
145: }
146: PetscFree(K);
148: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
149: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
151: /* Exclude any superfluous rows and columns */
152: nstart = 3*(2*m+1)*(2*m+1);
153: ict = 0;
154: PetscMalloc((N-nstart)*sizeof(PetscInt),&rowkeep);
155: for (i=nstart; i<N; i++) {
156: MatGetRow(mat,i,&nz,0,0);
157: if (nz) rowkeep[ict++] = i;
158: MatRestoreRow(mat,i,&nz,0,0);
159: }
160: ISCreateGeneral(PETSC_COMM_SELF,ict,rowkeep,&iskeep);
161: MatGetSubMatrices(mat,1,&iskeep,&iskeep,MAT_INITIAL_MATRIX,&submatb);
162: submat = *submatb;
163: PetscFree(submatb);
164: PetscFree(rowkeep);
165: ISDestroy(iskeep);
166: MatDestroy(mat);
168: /* Convert storage formats -- just to demonstrate conversion to various
169: formats (in particular, block diagonal storage). This is NOT the
170: recommended means to solve such a problem. */
171: MatConvert(submat,type,MAT_INITIAL_MATRIX,newmat);
172: MatDestroy(submat);
174: PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO);
175: MatView(*newmat,PETSC_VIEWER_STDOUT_WORLD);
176: MatNorm(*newmat,NORM_1,&norm);
177: PetscPrintf(PETSC_COMM_WORLD,"matrix 1 norm = %G\n",norm);
179: return 0;
180: }
181: /* -------------------------------------------------------------------- */
184: PetscErrorCode AddElement(Mat mat,PetscInt r1,PetscInt r2,PetscReal **K,PetscInt h1,PetscInt h2)
185: {
186: PetscScalar val;
187: PetscInt l1,l2,row,col;
190: for (l1=0; l1<3; l1++) {
191: for (l2=0; l2<3; l2++) {
192: /*
193: NOTE you should never do this! Inserting values 1 at a time is
194: just too expensive!
195: */
196: if (K[h1+l1][h2+l2] != 0.0) {
197: row = r1+l1; col = r2+l2; val = K[h1+l1][h2+l2];
198: MatSetValues(mat,1,&row,1,&col,&val,ADD_VALUES);
199: row = r2+l2; col = r1+l1;
200: MatSetValues(mat,1,&row,1,&col,&val,ADD_VALUES);
201: }
202: }
203: }
204: return 0;
205: }
206: /* -------------------------------------------------------------------- */
207: PetscReal N[20][64]; /* Interpolation function. */
208: PetscReal part_N[3][20][64]; /* Partials of interpolation function. */
209: PetscReal rst[3][64]; /* Location of integration pts in (r,s,t) */
210: PetscReal weight[64]; /* Gaussian quadrature weights. */
211: PetscReal xyz[20][3]; /* (x,y,z) coordinates of nodes */
212: PetscReal E,nu; /* Physcial constants. */
213: PetscInt n_int,N_int; /* N_int = n_int^3, number of int. pts. */
214: /* Ordering of the vertices, (r,s,t) coordinates, of the canonical cell. */
215: PetscReal r2[20] = {-1.0,0.0,1.0,-1.0,1.0,-1.0,0.0,1.0,
216: -1.0,1.0,-1.0,1.0,
217: -1.0,0.0,1.0,-1.0,1.0,-1.0,0.0,1.0};
218: PetscReal s2[20] = {-1.0,-1.0, -1.0,0.0,0.0,1.0, 1.0, 1.0,
219: -1.0,-1.0,1.0,1.0,
220: -1.0,-1.0, -1.0,0.0,0.0,1.0, 1.0, 1.0};
221: PetscReal t2[20] = {-1.0,-1.0,-1.0,-1.0,-1.0,-1.0,-1.0,-1.0,
222: 0.0,0.0,0.0,0.0,
223: 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
224: PetscInt rmap[20] = {0,1,2,3,5,6,7,8,9,11,15,17,18,19,20,21,23,24,25,26};
225: /* -------------------------------------------------------------------- */
228: /*
229: Elastic20Stiff - Forms 20 node elastic stiffness for element.
230: */
231: PetscErrorCode Elastic20Stiff(PetscReal **Ke)
232: {
233: PetscReal K[60][60],x,y,z,dx,dy,dz,m,v;
234: PetscInt i,j,k,l,Ii,J;
236: paulsetup20();
238: x = -1.0; y = -1.0; z = -1.0; dx = 2.0; dy = 2.0; dz = 2.0;
239: xyz[0][0] = x; xyz[0][1] = y; xyz[0][2] = z;
240: xyz[1][0] = x + dx; xyz[1][1] = y; xyz[1][2] = z;
241: xyz[2][0] = x + 2.*dx; xyz[2][1] = y; xyz[2][2] = z;
242: xyz[3][0] = x; xyz[3][1] = y + dy; xyz[3][2] = z;
243: xyz[4][0] = x + 2.*dx; xyz[4][1] = y + dy; xyz[4][2] = z;
244: xyz[5][0] = x; xyz[5][1] = y + 2.*dy; xyz[5][2] = z;
245: xyz[6][0] = x + dx; xyz[6][1] = y + 2.*dy; xyz[6][2] = z;
246: xyz[7][0] = x + 2.*dx; xyz[7][1] = y + 2.*dy; xyz[7][2] = z;
247: xyz[8][0] = x; xyz[8][1] = y; xyz[8][2] = z + dz;
248: xyz[9][0] = x + 2.*dx; xyz[9][1] = y; xyz[9][2] = z + dz;
249: xyz[10][0] = x; xyz[10][1] = y + 2.*dy; xyz[10][2] = z + dz;
250: xyz[11][0] = x + 2.*dx; xyz[11][1] = y + 2.*dy; xyz[11][2] = z + dz;
251: xyz[12][0] = x; xyz[12][1] = y; xyz[12][2] = z + 2.*dz;
252: xyz[13][0] = x + dx; xyz[13][1] = y; xyz[13][2] = z + 2.*dz;
253: xyz[14][0] = x + 2.*dx; xyz[14][1] = y; xyz[14][2] = z + 2.*dz;
254: xyz[15][0] = x; xyz[15][1] = y + dy; xyz[15][2] = z + 2.*dz;
255: xyz[16][0] = x + 2.*dx; xyz[16][1] = y + dy; xyz[16][2] = z + 2.*dz;
256: xyz[17][0] = x; xyz[17][1] = y + 2.*dy; xyz[17][2] = z + 2.*dz;
257: xyz[18][0] = x + dx; xyz[18][1] = y + 2.*dy; xyz[18][2] = z + 2.*dz;
258: xyz[19][0] = x + 2.*dx; xyz[19][1] = y + 2.*dy; xyz[19][2] = z + 2.*dz;
259: paulintegrate20(K);
261: /* copy the stiffness from K into format used by Ke */
262: for (i=0; i<81; i++) {
263: for (j=0; j<81; j++) {
264: Ke[i][j] = 0.0;
265: }
266: }
267: Ii = 0;
268: m = 0.0;
269: for (i=0; i<20; i++) {
270: J = 0;
271: for (j=0; j<20; j++) {
272: for (k=0; k<3; k++) {
273: for (l=0; l<3; l++) {
274: Ke[3*rmap[i]+k][3*rmap[j]+l] = v = K[Ii+k][J+l];
275: m = PetscMax(m,PetscAbsReal(v));
276: }
277: }
278: J += 3;
279: }
280: Ii += 3;
281: }
282: /* zero out the extremely small values */
283: m = (1.e-8)*m;
284: for (i=0; i<81; i++) {
285: for (j=0; j<81; j++) {
286: if (PetscAbsReal(Ke[i][j]) < m) Ke[i][j] = 0.0;
287: }
288: }
289: /* force the matrix to be exactly symmetric */
290: for (i=0; i<81; i++) {
291: for (j=0; j<i; j++) {
292: Ke[i][j] = (Ke[i][j] + Ke[j][i])/2.0;
293: }
294: }
295: return 0;
296: }
297: /* -------------------------------------------------------------------- */
300: /*
301: paulsetup20 - Sets up data structure for forming local elastic stiffness.
302: */
303: PetscErrorCode paulsetup20(void)
304: {
305: PetscInt i,j,k,cnt;
306: PetscReal x[4],w[4];
307: PetscReal c;
309: n_int = 3;
310: nu = 0.3;
311: E = 1.0;
313: /* Assign integration points and weights for
314: Gaussian quadrature formulae. */
315: if(n_int == 2) {
316: x[0] = (-0.577350269189626);
317: x[1] = (0.577350269189626);
318: w[0] = 1.0000000;
319: w[1] = 1.0000000;
320: }
321: else if(n_int == 3) {
322: x[0] = (-0.774596669241483);
323: x[1] = 0.0000000;
324: x[2] = 0.774596669241483;
325: w[0] = 0.555555555555555;
326: w[1] = 0.888888888888888;
327: w[2] = 0.555555555555555;
328: }
329: else if(n_int == 4) {
330: x[0] = (-0.861136311594053);
331: x[1] = (-0.339981043584856);
332: x[2] = 0.339981043584856;
333: x[3] = 0.861136311594053;
334: w[0] = 0.347854845137454;
335: w[1] = 0.652145154862546;
336: w[2] = 0.652145154862546;
337: w[3] = 0.347854845137454;
338: }
339: else {
340: SETERRQ(1,"Unknown value for n_int");
341: }
343: /* rst[][i] contains the location of the i-th integration point
344: in the canonical (r,s,t) coordinate system. weight[i] contains
345: the Gaussian weighting factor. */
347: cnt = 0;
348: for (i=0; i<n_int;i++) {
349: for (j=0; j<n_int;j++) {
350: for (k=0; k<n_int;k++) {
351: rst[0][cnt]=x[i];
352: rst[1][cnt]=x[j];
353: rst[2][cnt]=x[k];
354: weight[cnt] = w[i]*w[j]*w[k];
355: ++cnt;
356: }
357: }
358: }
359: N_int = cnt;
361: /* N[][j] is the interpolation vector, N[][j] .* xyz[] */
362: /* yields the (x,y,z) locations of the integration point. */
363: /* part_N[][][j] is the partials of the N function */
364: /* w.r.t. (r,s,t). */
366: c = 1.0/8.0;
367: for (j=0; j<N_int; j++) {
368: for (i=0; i<20; i++) {
369: if (i==0 || i==2 || i==5 || i==7 || i==12 || i==14 || i== 17 || i==19){
370: N[i][j] = c*(1.0 + r2[i]*rst[0][j])*
371: (1.0 + s2[i]*rst[1][j])*(1.0 + t2[i]*rst[2][j])*
372: (-2.0 + r2[i]*rst[0][j] + s2[i]*rst[1][j] + t2[i]*rst[2][j]);
373: part_N[0][i][j] = c*r2[i]*(1 + s2[i]*rst[1][j])*(1 + t2[i]*rst[2][j])*
374: (-1.0 + 2.0*r2[i]*rst[0][j] + s2[i]*rst[1][j] +
375: t2[i]*rst[2][j]);
376: part_N[1][i][j] = c*s2[i]*(1 + r2[i]*rst[0][j])*(1 + t2[i]*rst[2][j])*
377: (-1.0 + r2[i]*rst[0][j] + 2.0*s2[i]*rst[1][j] +
378: t2[i]*rst[2][j]);
379: part_N[2][i][j] = c*t2[i]*(1 + r2[i]*rst[0][j])*(1 + s2[i]*rst[1][j])*
380: (-1.0 + r2[i]*rst[0][j] + s2[i]*rst[1][j] +
381: 2.0*t2[i]*rst[2][j]);
382: }
383: else if (i==1 || i==6 || i==13 || i==18) {
384: N[i][j] = .25*(1.0 - rst[0][j]*rst[0][j])*
385: (1.0 + s2[i]*rst[1][j])*(1.0 + t2[i]*rst[2][j]);
386: part_N[0][i][j] = -.5*rst[0][j]*(1 + s2[i]*rst[1][j])*
387: (1 + t2[i]*rst[2][j]);
388: part_N[1][i][j] = .25*s2[i]*(1 + t2[i]*rst[2][j])*
389: (1.0 - rst[0][j]*rst[0][j]);
390: part_N[2][i][j] = .25*t2[i]*(1.0 - rst[0][j]*rst[0][j])*
391: (1 + s2[i]*rst[1][j]);
392: }
393: else if (i==3 || i==4 || i==15 || i==16) {
394: N[i][j] = .25*(1.0 - rst[1][j]*rst[1][j])*
395: (1.0 + r2[i]*rst[0][j])*(1.0 + t2[i]*rst[2][j]);
396: part_N[0][i][j] = .25*r2[i]*(1 + t2[i]*rst[2][j])*
397: (1.0 - rst[1][j]*rst[1][j]);
398: part_N[1][i][j] = -.5*rst[1][j]*(1 + r2[i]*rst[0][j])*
399: (1 + t2[i]*rst[2][j]);
400: part_N[2][i][j] = .25*t2[i]*(1.0 - rst[1][j]*rst[1][j])*
401: (1 + r2[i]*rst[0][j]);
402: }
403: else if (i==8 || i==9 || i==10 || i==11) {
404: N[i][j] = .25*(1.0 - rst[2][j]*rst[2][j])*
405: (1.0 + r2[i]*rst[0][j])*(1.0 + s2[i]*rst[1][j]);
406: part_N[0][i][j] = .25*r2[i]*(1 + s2[i]*rst[1][j])*
407: (1.0 - rst[2][j]*rst[2][j]);
408: part_N[1][i][j] = .25*s2[i]*(1.0 - rst[2][j]*rst[2][j])*
409: (1 + r2[i]*rst[0][j]);
410: part_N[2][i][j] = -.5*rst[2][j]*(1 + r2[i]*rst[0][j])*
411: (1 + s2[i]*rst[1][j]);
412: }
413: }
414: }
415: return 0;
416: }
417: /* -------------------------------------------------------------------- */
420: /*
421: paulintegrate20 - Does actual numerical integration on 20 node element.
422: */
423: PetscErrorCode paulintegrate20(PetscReal K[60][60])
424: {
425: PetscReal det_jac,jac[3][3],inv_jac[3][3];
426: PetscReal B[6][60],B_temp[6][60],C[6][6];
427: PetscReal temp;
428: PetscInt i,j,k,step;
430: /* Zero out K, since we will accumulate the result here */
431: for (i=0; i<60; i++) {
432: for (j=0; j<60; j++) {
433: K[i][j] = 0.0;
434: }
435: }
437: /* Loop over integration points ... */
438: for (step=0; step<N_int; step++) {
440: /* Compute the Jacobian, its determinant, and inverse. */
441: for (i=0; i<3; i++) {
442: for (j=0; j<3; j++) {
443: jac[i][j] = 0;
444: for (k=0; k<20; k++) {
445: jac[i][j] += part_N[i][k][step]*xyz[k][j];
446: }
447: }
448: }
449: det_jac = jac[0][0]*(jac[1][1]*jac[2][2]-jac[1][2]*jac[2][1])
450: + jac[0][1]*(jac[1][2]*jac[2][0]-jac[1][0]*jac[2][2])
451: + jac[0][2]*(jac[1][0]*jac[2][1]-jac[1][1]*jac[2][0]);
452: inv_jac[0][0] = (jac[1][1]*jac[2][2]-jac[1][2]*jac[2][1])/det_jac;
453: inv_jac[0][1] = (jac[0][2]*jac[2][1]-jac[0][1]*jac[2][2])/det_jac;
454: inv_jac[0][2] = (jac[0][1]*jac[1][2]-jac[1][1]*jac[0][2])/det_jac;
455: inv_jac[1][0] = (jac[1][2]*jac[2][0]-jac[1][0]*jac[2][2])/det_jac;
456: inv_jac[1][1] = (jac[0][0]*jac[2][2]-jac[2][0]*jac[0][2])/det_jac;
457: inv_jac[1][2] = (jac[0][2]*jac[1][0]-jac[0][0]*jac[1][2])/det_jac;
458: inv_jac[2][0] = (jac[1][0]*jac[2][1]-jac[1][1]*jac[2][0])/det_jac;
459: inv_jac[2][1] = (jac[0][1]*jac[2][0]-jac[0][0]*jac[2][1])/det_jac;
460: inv_jac[2][2] = (jac[0][0]*jac[1][1]-jac[1][0]*jac[0][1])/det_jac;
462: /* Compute the B matrix. */
463: for (i=0; i<3; i++) {
464: for (j=0; j<20; j++) {
465: B_temp[i][j] = 0.0;
466: for (k=0; k<3; k++) {
467: B_temp[i][j] += inv_jac[i][k]*part_N[k][j][step];
468: }
469: }
470: }
471: for (i=0; i<6; i++) {
472: for (j=0; j<60; j++) {
473: B[i][j] = 0.0;
474: }
475: }
477: /* Put values in correct places in B. */
478: for (k=0; k<20; k++) {
479: B[0][3*k] = B_temp[0][k];
480: B[1][3*k+1] = B_temp[1][k];
481: B[2][3*k+2] = B_temp[2][k];
482: B[3][3*k] = B_temp[1][k];
483: B[3][3*k+1] = B_temp[0][k];
484: B[4][3*k+1] = B_temp[2][k];
485: B[4][3*k+2] = B_temp[1][k];
486: B[5][3*k] = B_temp[2][k];
487: B[5][3*k+2] = B_temp[0][k];
488: }
489:
490: /* Construct the C matrix, uses the constants "nu" and "E". */
491: for (i=0; i<6; i++) {
492: for (j=0; j<6; j++) {
493: C[i][j] = 0.0;
494: }
495: }
496: temp = (1.0 + nu)*(1.0 - 2.0*nu);
497: temp = E/temp;
498: C[0][0] = temp*(1.0 - nu);
499: C[1][1] = C[0][0];
500: C[2][2] = C[0][0];
501: C[3][3] = temp*(0.5 - nu);
502: C[4][4] = C[3][3];
503: C[5][5] = C[3][3];
504: C[0][1] = temp*nu;
505: C[0][2] = C[0][1];
506: C[1][0] = C[0][1];
507: C[1][2] = C[0][1];
508: C[2][0] = C[0][1];
509: C[2][1] = C[0][1];
510:
511: for (i=0; i<6; i++) {
512: for (j=0; j<60; j++) {
513: B_temp[i][j] = 0.0;
514: for (k=0; k<6; k++) {
515: B_temp[i][j] += C[i][k]*B[k][j];
516: }
517: B_temp[i][j] *= det_jac;
518: }
519: }
521: /* Accumulate B'*C*B*det(J)*weight, as a function of (r,s,t), in K. */
522: for (i=0; i<60; i++) {
523: for (j=0; j<60; j++) {
524: temp = 0.0;
525: for (k=0; k<6; k++) {
526: temp += B[k][i]*B_temp[k][j];
527: }
528: K[i][j] += temp*weight[step];
529: }
530: }
531: } /* end of loop over integration points */
532: return 0;
533: }