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