Actual source code: ex10.c

  1: static char help[] = "Linear elastiticty with dimensions using 20 node serendipity elements.\n\
  2: This also demonstrates use of  block\n\
  3: diagonal data structure.  Input arguments are:\n\
  4:   -m : problem size\n\n";

  6: #include <petscksp.h>

  8: /* This code is not intended as an efficient implementation, it is only
  9:    here to produce an interesting sparse matrix quickly.

 11:    PLEASE DO NOT BASE ANY OF YOUR CODES ON CODE LIKE THIS, THERE ARE MUCH
 12:    BETTER WAYS TO DO THIS. */

 14: extern PetscErrorCode GetElasticityMatrix(PetscInt, Mat *);
 15: extern PetscErrorCode Elastic20Stiff(PetscReal **);
 16: extern PetscErrorCode AddElement(Mat, PetscInt, PetscInt, PetscReal **, PetscInt, PetscInt);
 17: extern PetscErrorCode paulsetup20(void);
 18: extern PetscErrorCode paulintegrate20(PetscReal K[60][60]);

 20: int main(int argc, char **args)
 21: {
 22:   Mat         mat;
 23:   PetscInt    i, its, m = 3, rdim, cdim, rstart, rend;
 24:   PetscMPIInt rank, size;
 25:   PetscScalar v, neg1 = -1.0;
 26:   Vec         u, x, b;
 27:   KSP         ksp;
 28:   PetscReal   norm;

 30:   PetscFunctionBeginUser;
 31:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 32:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 33:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 34:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 36:   /* Form matrix */
 37:   PetscCall(GetElasticityMatrix(m, &mat));

 39:   /* Generate vectors */
 40:   PetscCall(MatGetSize(mat, &rdim, &cdim));
 41:   PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
 42:   PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
 43:   PetscCall(VecSetSizes(u, PETSC_DECIDE, rdim));
 44:   PetscCall(VecSetFromOptions(u));
 45:   PetscCall(VecDuplicate(u, &b));
 46:   PetscCall(VecDuplicate(b, &x));
 47:   for (i = rstart; i < rend; i++) {
 48:     v = (PetscScalar)(i - rstart + 100 * rank);
 49:     PetscCall(VecSetValues(u, 1, &i, &v, INSERT_VALUES));
 50:   }
 51:   PetscCall(VecAssemblyBegin(u));
 52:   PetscCall(VecAssemblyEnd(u));

 54:   /* Compute right-hand side */
 55:   PetscCall(MatMult(mat, u, b));

 57:   /* Solve linear system */
 58:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 59:   PetscCall(KSPSetOperators(ksp, mat, mat));
 60:   PetscCall(KSPSetFromOptions(ksp));
 61:   PetscCall(KSPSolve(ksp, b, x));
 62:   PetscCall(KSPGetIterationNumber(ksp, &its));
 63:   /* Check error */
 64:   PetscCall(VecAXPY(x, neg1, u));
 65:   PetscCall(VecNorm(x, NORM_2, &norm));

 67:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of residual %g Number of iterations %" PetscInt_FMT "\n", (double)norm, its));

 69:   /* Free work space */
 70:   PetscCall(KSPDestroy(&ksp));
 71:   PetscCall(VecDestroy(&u));
 72:   PetscCall(VecDestroy(&x));
 73:   PetscCall(VecDestroy(&b));
 74:   PetscCall(MatDestroy(&mat));

 76:   PetscCall(PetscFinalize());
 77:   return 0;
 78: }
 79: /* -------------------------------------------------------------------- */
 80: /*
 81:   GetElasticityMatrix - Forms 3D linear elasticity matrix.
 82:  */
 83: PetscErrorCode GetElasticityMatrix(PetscInt m, Mat *newmat)
 84: {
 85:   PetscInt    i, j, k, i1, i2, j_1, j2, k1, k2, h1, h2, shiftx, shifty, shiftz;
 86:   PetscInt    ict, nz, base, r1, r2, N, *rowkeep, nstart;
 87:   IS          iskeep;
 88:   PetscReal **K, norm;
 89:   Mat         mat, submat = 0, *submatb;
 90:   MatType     type = MATSEQBAIJ;

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

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

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

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

136:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
137:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

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

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

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

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

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

214:   PetscCall(paulsetup20());

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

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

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

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

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

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

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

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

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

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

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

432:     /* Compute the B matrix. */
433:     for (i = 0; i < 3; i++) {
434:       for (j = 0; j < 20; j++) {
435:         B_temp[i][j] = 0.0;
436:         for (k = 0; k < 3; k++) B_temp[i][j] += inv_jac[i][k] * part_N[k][j][step];
437:       }
438:     }
439:     for (i = 0; i < 6; i++) {
440:       for (j = 0; j < 60; j++) B[i][j] = 0.0;
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++) C[i][j] = 0.0;
459:     }
460:     temp    = (1.0 + nu) * (1.0 - 2.0 * nu);
461:     temp    = E / temp;
462:     C[0][0] = temp * (1.0 - nu);
463:     C[1][1] = C[0][0];
464:     C[2][2] = C[0][0];
465:     C[3][3] = temp * (0.5 - nu);
466:     C[4][4] = C[3][3];
467:     C[5][5] = C[3][3];
468:     C[0][1] = temp * nu;
469:     C[0][2] = C[0][1];
470:     C[1][0] = C[0][1];
471:     C[1][2] = C[0][1];
472:     C[2][0] = C[0][1];
473:     C[2][1] = C[0][1];

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

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

495: /*TEST

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

501: TEST*/