Actual source code: ex10.c


  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;
 24:   PetscInt    i, its, m = 3, rdim, cdim, rstart, rend;
 25:   PetscMPIInt rank, size;
 26:   PetscScalar v, neg1 = -1.0;
 27:   Vec         u, x, b;
 28:   KSP         ksp;
 29:   PetscReal   norm;

 32:   PetscInitialize(&argc, &args, (char *)0, help);
 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 %" PetscInt_FMT "\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 0;
 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;
 88:   IS          iskeep;
 89:   PetscReal **K, norm;
 90:   Mat         mat, submat = 0, *submatb;
 91:   MatType     type = MATSEQBAIJ;

 93:   m /= 2; /* This is done just to be consistent with the old example */
 94:   N = 3 * (2 * m + 1) * (2 * m + 1) * (2 * m + 1);
 95:   PetscPrintf(PETSC_COMM_SELF, "m = %" PetscInt_FMT ", N=%" PetscInt_FMT "\n", m, N);
 96:   MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 80, NULL, &mat);

 98:   /* Form stiffness for element */
 99:   PetscMalloc1(81, &K);
100:   for (i = 0; i < 81; i++) PetscMalloc1(81, &K[i]);
101:   Elastic20Stiff(K);

103:   /* Loop over elements and add contribution to stiffness */
104:   shiftx = 3;
105:   shifty = 3 * (2 * m + 1);
106:   shiftz = 3 * (2 * m + 1) * (2 * m + 1);
107:   for (k = 0; k < m; k++) {
108:     for (j = 0; j < m; j++) {
109:       for (i = 0; i < m; i++) {
110:         h1   = 0;
111:         base = 2 * k * shiftz + 2 * j * shifty + 2 * i * shiftx;
112:         for (k1 = 0; k1 < 3; k1++) {
113:           for (j_1 = 0; j_1 < 3; j_1++) {
114:             for (i1 = 0; i1 < 3; i1++) {
115:               h2 = 0;
116:               r1 = base + i1 * shiftx + j_1 * shifty + k1 * shiftz;
117:               for (k2 = 0; k2 < 3; k2++) {
118:                 for (j2 = 0; j2 < 3; j2++) {
119:                   for (i2 = 0; i2 < 3; i2++) {
120:                     r2 = base + i2 * shiftx + j2 * shifty + k2 * shiftz;
121:                     AddElement(mat, r1, r2, K, h1, h2);
122:                     h2 += 3;
123:                   }
124:                 }
125:               }
126:               h1 += 3;
127:             }
128:           }
129:         }
130:       }
131:     }
132:   }

134:   for (i = 0; i < 81; i++) PetscFree(K[i]);
135:   PetscFree(K);

137:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
138:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);

140:   /* Exclude any superfluous rows and columns */
141:   nstart = 3 * (2 * m + 1) * (2 * m + 1);
142:   ict    = 0;
143:   PetscMalloc1(N - nstart, &rowkeep);
144:   for (i = nstart; i < N; i++) {
145:     MatGetRow(mat, i, &nz, 0, 0);
146:     if (nz) rowkeep[ict++] = i;
147:     MatRestoreRow(mat, i, &nz, 0, 0);
148:   }
149:   ISCreateGeneral(PETSC_COMM_SELF, ict, rowkeep, PETSC_COPY_VALUES, &iskeep);
150:   MatCreateSubMatrices(mat, 1, &iskeep, &iskeep, MAT_INITIAL_MATRIX, &submatb);
151:   submat = *submatb;
152:   PetscFree(submatb);
153:   PetscFree(rowkeep);
154:   ISDestroy(&iskeep);
155:   MatDestroy(&mat);

157:   /* Convert storage formats -- just to demonstrate conversion to various
158:      formats (in particular, block diagonal storage).  This is NOT the
159:      recommended means to solve such a problem.  */
160:   MatConvert(submat, type, MAT_INITIAL_MATRIX, newmat);
161:   MatDestroy(&submat);

163:   MatNorm(*newmat, NORM_1, &norm);
164:   PetscPrintf(PETSC_COMM_WORLD, "matrix 1 norm = %g\n", (double)norm);

166:   return 0;
167: }
168: /* -------------------------------------------------------------------- */
169: PetscErrorCode AddElement(Mat mat, PetscInt r1, PetscInt r2, PetscReal **K, PetscInt h1, PetscInt h2)
170: {
171:   PetscScalar val;
172:   PetscInt    l1, l2, row, col;

174:   for (l1 = 0; l1 < 3; l1++) {
175:     for (l2 = 0; l2 < 3; l2++) {
176:       /*
177:    NOTE you should never do this! Inserting values 1 at a time is
178:    just too expensive!
179: */
180:       if (K[h1 + l1][h2 + l2] != 0.0) {
181:         row = r1 + l1;
182:         col = r2 + l2;
183:         val = K[h1 + l1][h2 + l2];
184:         MatSetValues(mat, 1, &row, 1, &col, &val, ADD_VALUES);
185:         row = r2 + l2;
186:         col = r1 + l1;
187:         MatSetValues(mat, 1, &row, 1, &col, &val, ADD_VALUES);
188:       }
189:     }
190:   }
191:   return 0;
192: }
193: /* -------------------------------------------------------------------- */
194: PetscReal N[20][64];         /* Interpolation function. */
195: PetscReal part_N[3][20][64]; /* Partials of interpolation function. */
196: PetscReal rst[3][64];        /* Location of integration pts in (r,s,t) */
197: PetscReal weight[64];        /* Gaussian quadrature weights. */
198: PetscReal xyz[20][3];        /* (x,y,z) coordinates of nodes  */
199: PetscReal E, nu;             /* Physical constants. */
200: PetscInt  n_int, N_int;      /* N_int = n_int^3, number of int. pts. */
201: /* Ordering of the vertices, (r,s,t) coordinates, of the canonical cell. */
202: PetscReal r2[20]   = {-1.0, 0.0, 1.0, -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 0.0, 1.0, -1.0, 1.0, -1.0, 0.0, 1.0};
203: PetscReal s2[20]   = {-1.0, -1.0, -1.0, 0.0, 0.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 1.0, 1.0, 1.0};
204: PetscReal t2[20]   = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
205: PetscInt  rmap[20] = {0, 1, 2, 3, 5, 6, 7, 8, 9, 11, 15, 17, 18, 19, 20, 21, 23, 24, 25, 26};
206: /* -------------------------------------------------------------------- */
207: /*
208:   Elastic20Stiff - Forms 20 node elastic stiffness for element.
209:  */
210: PetscErrorCode Elastic20Stiff(PetscReal **Ke)
211: {
212:   PetscReal K[60][60], x, y, z, dx, dy, dz;
213:   PetscInt  i, j, k, l, Ii, J;

215:   paulsetup20();

217:   x          = -1.0;
218:   y          = -1.0;
219:   z          = -1.0;
220:   dx         = 2.0;
221:   dy         = 2.0;
222:   dz         = 2.0;
223:   xyz[0][0]  = x;
224:   xyz[0][1]  = y;
225:   xyz[0][2]  = z;
226:   xyz[1][0]  = x + dx;
227:   xyz[1][1]  = y;
228:   xyz[1][2]  = z;
229:   xyz[2][0]  = x + 2. * dx;
230:   xyz[2][1]  = y;
231:   xyz[2][2]  = z;
232:   xyz[3][0]  = x;
233:   xyz[3][1]  = y + dy;
234:   xyz[3][2]  = z;
235:   xyz[4][0]  = x + 2. * dx;
236:   xyz[4][1]  = y + dy;
237:   xyz[4][2]  = z;
238:   xyz[5][0]  = x;
239:   xyz[5][1]  = y + 2. * dy;
240:   xyz[5][2]  = z;
241:   xyz[6][0]  = x + dx;
242:   xyz[6][1]  = y + 2. * dy;
243:   xyz[6][2]  = z;
244:   xyz[7][0]  = x + 2. * dx;
245:   xyz[7][1]  = y + 2. * dy;
246:   xyz[7][2]  = z;
247:   xyz[8][0]  = x;
248:   xyz[8][1]  = y;
249:   xyz[8][2]  = z + dz;
250:   xyz[9][0]  = x + 2. * dx;
251:   xyz[9][1]  = y;
252:   xyz[9][2]  = z + dz;
253:   xyz[10][0] = x;
254:   xyz[10][1] = y + 2. * dy;
255:   xyz[10][2] = z + dz;
256:   xyz[11][0] = x + 2. * dx;
257:   xyz[11][1] = y + 2. * dy;
258:   xyz[11][2] = z + dz;
259:   xyz[12][0] = x;
260:   xyz[12][1] = y;
261:   xyz[12][2] = z + 2. * dz;
262:   xyz[13][0] = x + dx;
263:   xyz[13][1] = y;
264:   xyz[13][2] = z + 2. * dz;
265:   xyz[14][0] = x + 2. * dx;
266:   xyz[14][1] = y;
267:   xyz[14][2] = z + 2. * dz;
268:   xyz[15][0] = x;
269:   xyz[15][1] = y + dy;
270:   xyz[15][2] = z + 2. * dz;
271:   xyz[16][0] = x + 2. * dx;
272:   xyz[16][1] = y + dy;
273:   xyz[16][2] = z + 2. * dz;
274:   xyz[17][0] = x;
275:   xyz[17][1] = y + 2. * dy;
276:   xyz[17][2] = z + 2. * dz;
277:   xyz[18][0] = x + dx;
278:   xyz[18][1] = y + 2. * dy;
279:   xyz[18][2] = z + 2. * dz;
280:   xyz[19][0] = x + 2. * dx;
281:   xyz[19][1] = y + 2. * dy;
282:   xyz[19][2] = z + 2. * dz;
283:   paulintegrate20(K);

285:   /* copy the stiffness from K into format used by Ke */
286:   for (i = 0; i < 81; i++) {
287:     for (j = 0; j < 81; j++) Ke[i][j] = 0.0;
288:   }
289:   Ii = 0;
290:   for (i = 0; i < 20; i++) {
291:     J = 0;
292:     for (j = 0; j < 20; j++) {
293:       for (k = 0; k < 3; k++) {
294:         for (l = 0; l < 3; l++) Ke[3 * rmap[i] + k][3 * rmap[j] + l] = K[Ii + k][J + l];
295:       }
296:       J += 3;
297:     }
298:     Ii += 3;
299:   }

301:   /* force the matrix to be exactly symmetric */
302:   for (i = 0; i < 81; i++) {
303:     for (j = 0; j < i; j++) Ke[i][j] = (Ke[i][j] + Ke[j][i]) / 2.0;
304:   }
305:   return 0;
306: }
307: /* -------------------------------------------------------------------- */
308: /*
309:   paulsetup20 - Sets up data structure for forming local elastic stiffness.
310:  */
311: PetscErrorCode paulsetup20(void)
312: {
313:   PetscInt  i, j, k, cnt;
314:   PetscReal x[4], w[4];
315:   PetscReal c;

317:   n_int = 3;
318:   nu    = 0.3;
319:   E     = 1.0;

321:   /* Assign integration points and weights for
322:        Gaussian quadrature formulae. */
323:   if (n_int == 2) {
324:     x[0] = (-0.577350269189626);
325:     x[1] = (0.577350269189626);
326:     w[0] = 1.0000000;
327:     w[1] = 1.0000000;
328:   } else if (n_int == 3) {
329:     x[0] = (-0.774596669241483);
330:     x[1] = 0.0000000;
331:     x[2] = 0.774596669241483;
332:     w[0] = 0.555555555555555;
333:     w[1] = 0.888888888888888;
334:     w[2] = 0.555555555555555;
335:   } else if (n_int == 4) {
336:     x[0] = (-0.861136311594053);
337:     x[1] = (-0.339981043584856);
338:     x[2] = 0.339981043584856;
339:     x[3] = 0.861136311594053;
340:     w[0] = 0.347854845137454;
341:     w[1] = 0.652145154862546;
342:     w[2] = 0.652145154862546;
343:     w[3] = 0.347854845137454;
344:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Value for n_int not supported");

346:   /* rst[][i] contains the location of the i-th integration point
347:       in the canonical (r,s,t) coordinate system.  weight[i] contains
348:       the Gaussian weighting factor. */

350:   cnt = 0;
351:   for (i = 0; i < n_int; i++) {
352:     for (j = 0; j < n_int; j++) {
353:       for (k = 0; k < n_int; k++) {
354:         rst[0][cnt] = x[i];
355:         rst[1][cnt] = x[j];
356:         rst[2][cnt] = x[k];
357:         weight[cnt] = w[i] * w[j] * w[k];
358:         ++cnt;
359:       }
360:     }
361:   }
362:   N_int = cnt;

364:   /* N[][j] is the interpolation vector, N[][j] .* xyz[] */
365:   /* yields the (x,y,z)  locations of the integration point. */
366:   /*  part_N[][][j] is the partials of the N function */
367:   /*  w.r.t. (r,s,t). */

369:   c = 1.0 / 8.0;
370:   for (j = 0; j < N_int; j++) {
371:     for (i = 0; i < 20; i++) {
372:       if (i == 0 || i == 2 || i == 5 || i == 7 || i == 12 || i == 14 || i == 17 || i == 19) {
373:         N[i][j]         = c * (1.0 + r2[i] * rst[0][j]) * (1.0 + s2[i] * rst[1][j]) * (1.0 + t2[i] * rst[2][j]) * (-2.0 + r2[i] * rst[0][j] + s2[i] * rst[1][j] + t2[i] * rst[2][j]);
374:         part_N[0][i][j] = c * r2[i] * (1 + s2[i] * rst[1][j]) * (1 + t2[i] * rst[2][j]) * (-1.0 + 2.0 * r2[i] * rst[0][j] + s2[i] * rst[1][j] + t2[i] * rst[2][j]);
375:         part_N[1][i][j] = c * s2[i] * (1 + r2[i] * rst[0][j]) * (1 + t2[i] * rst[2][j]) * (-1.0 + r2[i] * rst[0][j] + 2.0 * s2[i] * rst[1][j] + t2[i] * rst[2][j]);
376:         part_N[2][i][j] = c * t2[i] * (1 + r2[i] * rst[0][j]) * (1 + s2[i] * rst[1][j]) * (-1.0 + r2[i] * rst[0][j] + s2[i] * rst[1][j] + 2.0 * t2[i] * rst[2][j]);
377:       } else if (i == 1 || i == 6 || i == 13 || i == 18) {
378:         N[i][j]         = .25 * (1.0 - rst[0][j] * rst[0][j]) * (1.0 + s2[i] * rst[1][j]) * (1.0 + t2[i] * rst[2][j]);
379:         part_N[0][i][j] = -.5 * rst[0][j] * (1 + s2[i] * rst[1][j]) * (1 + t2[i] * rst[2][j]);
380:         part_N[1][i][j] = .25 * s2[i] * (1 + t2[i] * rst[2][j]) * (1.0 - rst[0][j] * rst[0][j]);
381:         part_N[2][i][j] = .25 * t2[i] * (1.0 - rst[0][j] * rst[0][j]) * (1 + s2[i] * rst[1][j]);
382:       } else if (i == 3 || i == 4 || i == 15 || i == 16) {
383:         N[i][j]         = .25 * (1.0 - rst[1][j] * rst[1][j]) * (1.0 + r2[i] * rst[0][j]) * (1.0 + t2[i] * rst[2][j]);
384:         part_N[0][i][j] = .25 * r2[i] * (1 + t2[i] * rst[2][j]) * (1.0 - rst[1][j] * rst[1][j]);
385:         part_N[1][i][j] = -.5 * rst[1][j] * (1 + r2[i] * rst[0][j]) * (1 + t2[i] * rst[2][j]);
386:         part_N[2][i][j] = .25 * t2[i] * (1.0 - rst[1][j] * rst[1][j]) * (1 + r2[i] * rst[0][j]);
387:       } else if (i == 8 || i == 9 || i == 10 || i == 11) {
388:         N[i][j]         = .25 * (1.0 - rst[2][j] * rst[2][j]) * (1.0 + r2[i] * rst[0][j]) * (1.0 + s2[i] * rst[1][j]);
389:         part_N[0][i][j] = .25 * r2[i] * (1 + s2[i] * rst[1][j]) * (1.0 - rst[2][j] * rst[2][j]);
390:         part_N[1][i][j] = .25 * s2[i] * (1.0 - rst[2][j] * rst[2][j]) * (1 + r2[i] * rst[0][j]);
391:         part_N[2][i][j] = -.5 * rst[2][j] * (1 + r2[i] * rst[0][j]) * (1 + s2[i] * rst[1][j]);
392:       }
393:     }
394:   }
395:   return 0;
396: }
397: /* -------------------------------------------------------------------- */
398: /*
399:    paulintegrate20 - Does actual numerical integration on 20 node element.
400:  */
401: PetscErrorCode paulintegrate20(PetscReal K[60][60])
402: {
403:   PetscReal det_jac, jac[3][3], inv_jac[3][3];
404:   PetscReal B[6][60], B_temp[6][60], C[6][6];
405:   PetscReal temp;
406:   PetscInt  i, j, k, step;

408:   /* Zero out K, since we will accumulate the result here */
409:   for (i = 0; i < 60; i++) {
410:     for (j = 0; j < 60; j++) K[i][j] = 0.0;
411:   }

413:   /* Loop over integration points ... */
414:   for (step = 0; step < N_int; step++) {
415:     /* Compute the Jacobian, its determinant, and inverse. */
416:     for (i = 0; i < 3; i++) {
417:       for (j = 0; j < 3; j++) {
418:         jac[i][j] = 0;
419:         for (k = 0; k < 20; k++) jac[i][j] += part_N[i][k][step] * xyz[k][j];
420:       }
421:     }
422:     det_jac       = jac[0][0] * (jac[1][1] * jac[2][2] - jac[1][2] * jac[2][1]) + jac[0][1] * (jac[1][2] * jac[2][0] - jac[1][0] * jac[2][2]) + jac[0][2] * (jac[1][0] * jac[2][1] - jac[1][1] * jac[2][0]);
423:     inv_jac[0][0] = (jac[1][1] * jac[2][2] - jac[1][2] * jac[2][1]) / det_jac;
424:     inv_jac[0][1] = (jac[0][2] * jac[2][1] - jac[0][1] * jac[2][2]) / det_jac;
425:     inv_jac[0][2] = (jac[0][1] * jac[1][2] - jac[1][1] * jac[0][2]) / det_jac;
426:     inv_jac[1][0] = (jac[1][2] * jac[2][0] - jac[1][0] * jac[2][2]) / det_jac;
427:     inv_jac[1][1] = (jac[0][0] * jac[2][2] - jac[2][0] * jac[0][2]) / det_jac;
428:     inv_jac[1][2] = (jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) / det_jac;
429:     inv_jac[2][0] = (jac[1][0] * jac[2][1] - jac[1][1] * jac[2][0]) / det_jac;
430:     inv_jac[2][1] = (jac[0][1] * jac[2][0] - jac[0][0] * jac[2][1]) / det_jac;
431:     inv_jac[2][2] = (jac[0][0] * jac[1][1] - jac[1][0] * jac[0][1]) / det_jac;

433:     /* Compute the B matrix. */
434:     for (i = 0; i < 3; i++) {
435:       for (j = 0; j < 20; j++) {
436:         B_temp[i][j] = 0.0;
437:         for (k = 0; k < 3; k++) B_temp[i][j] += inv_jac[i][k] * part_N[k][j][step];
438:       }
439:     }
440:     for (i = 0; i < 6; i++) {
441:       for (j = 0; j < 60; j++) B[i][j] = 0.0;
442:     }

444:     /* Put values in correct places in B. */
445:     for (k = 0; k < 20; k++) {
446:       B[0][3 * k]     = B_temp[0][k];
447:       B[1][3 * k + 1] = B_temp[1][k];
448:       B[2][3 * k + 2] = B_temp[2][k];
449:       B[3][3 * k]     = B_temp[1][k];
450:       B[3][3 * k + 1] = B_temp[0][k];
451:       B[4][3 * k + 1] = B_temp[2][k];
452:       B[4][3 * k + 2] = B_temp[1][k];
453:       B[5][3 * k]     = B_temp[2][k];
454:       B[5][3 * k + 2] = B_temp[0][k];
455:     }

457:     /* Construct the C matrix, uses the constants "nu" and "E". */
458:     for (i = 0; i < 6; i++) {
459:       for (j = 0; j < 6; j++) C[i][j] = 0.0;
460:     }
461:     temp    = (1.0 + nu) * (1.0 - 2.0 * nu);
462:     temp    = E / temp;
463:     C[0][0] = temp * (1.0 - nu);
464:     C[1][1] = C[0][0];
465:     C[2][2] = C[0][0];
466:     C[3][3] = temp * (0.5 - nu);
467:     C[4][4] = C[3][3];
468:     C[5][5] = C[3][3];
469:     C[0][1] = temp * nu;
470:     C[0][2] = C[0][1];
471:     C[1][0] = C[0][1];
472:     C[1][2] = C[0][1];
473:     C[2][0] = C[0][1];
474:     C[2][1] = C[0][1];

476:     for (i = 0; i < 6; i++) {
477:       for (j = 0; j < 60; j++) {
478:         B_temp[i][j] = 0.0;
479:         for (k = 0; k < 6; k++) B_temp[i][j] += C[i][k] * B[k][j];
480:         B_temp[i][j] *= det_jac;
481:       }
482:     }

484:     /* Accumulate B'*C*B*det(J)*weight, as a function of (r,s,t), in K. */
485:     for (i = 0; i < 60; i++) {
486:       for (j = 0; j < 60; j++) {
487:         temp = 0.0;
488:         for (k = 0; k < 6; k++) temp += B[k][i] * B_temp[k][j];
489:         K[i][j] += temp * weight[step];
490:       }
491:     }
492:   } /* end of loop over integration points */
493:   return 0;
494: }

496: /*TEST

498:     test:
499:       args: -matconvert_type seqaij -ksp_monitor_short -ksp_rtol 1.e-2  -pc_type jacobi
500:       requires: x

502: TEST*/