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*/