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;
31: PetscFunctionBeginUser;
32: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
33: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
34: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
35: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
37: /* Form matrix */
38: PetscCall(GetElasticityMatrix(m, &mat));
40: /* Generate vectors */
41: PetscCall(MatGetSize(mat, &rdim, &cdim));
42: PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
43: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
44: PetscCall(VecSetSizes(u, PETSC_DECIDE, rdim));
45: PetscCall(VecSetFromOptions(u));
46: PetscCall(VecDuplicate(u, &b));
47: PetscCall(VecDuplicate(b, &x));
48: for (i = rstart; i < rend; i++) {
49: v = (PetscScalar)(i - rstart + 100 * rank);
50: PetscCall(VecSetValues(u, 1, &i, &v, INSERT_VALUES));
51: }
52: PetscCall(VecAssemblyBegin(u));
53: PetscCall(VecAssemblyEnd(u));
55: /* Compute right-hand-side */
56: PetscCall(MatMult(mat, u, b));
58: /* Solve linear system */
59: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
60: PetscCall(KSPSetOperators(ksp, mat, mat));
61: PetscCall(KSPSetFromOptions(ksp));
62: PetscCall(KSPSolve(ksp, b, x));
63: PetscCall(KSPGetIterationNumber(ksp, &its));
64: /* Check error */
65: PetscCall(VecAXPY(x, neg1, u));
66: PetscCall(VecNorm(x, NORM_2, &norm));
68: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of residual %g Number of iterations %" PetscInt_FMT "\n", (double)norm, its));
70: /* Free work space */
71: PetscCall(KSPDestroy(&ksp));
72: PetscCall(VecDestroy(&u));
73: PetscCall(VecDestroy(&x));
74: PetscCall(VecDestroy(&b));
75: PetscCall(MatDestroy(&mat));
77: PetscCall(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: PetscCall(PetscPrintf(PETSC_COMM_SELF, "m = %" PetscInt_FMT ", N=%" PetscInt_FMT "\n", m, N));
96: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 80, NULL, &mat));
98: /* Form stiffness for element */
99: PetscCall(PetscMalloc1(81, &K));
100: for (i = 0; i < 81; i++) PetscCall(PetscMalloc1(81, &K[i]));
101: PetscCall(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: PetscCall(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++) PetscCall(PetscFree(K[i]));
135: PetscCall(PetscFree(K));
137: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
138: PetscCall(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: PetscCall(PetscMalloc1(N - nstart, &rowkeep));
144: for (i = nstart; i < N; i++) {
145: PetscCall(MatGetRow(mat, i, &nz, 0, 0));
146: if (nz) rowkeep[ict++] = i;
147: PetscCall(MatRestoreRow(mat, i, &nz, 0, 0));
148: }
149: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ict, rowkeep, PETSC_COPY_VALUES, &iskeep));
150: PetscCall(MatCreateSubMatrices(mat, 1, &iskeep, &iskeep, MAT_INITIAL_MATRIX, &submatb));
151: submat = *submatb;
152: PetscCall(PetscFree(submatb));
153: PetscCall(PetscFree(rowkeep));
154: PetscCall(ISDestroy(&iskeep));
155: PetscCall(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: PetscCall(MatConvert(submat, type, MAT_INITIAL_MATRIX, newmat));
161: PetscCall(MatDestroy(&submat));
163: PetscCall(MatNorm(*newmat, NORM_1, &norm));
164: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "matrix 1 norm = %g\n", (double)norm));
166: return PETSC_SUCCESS;
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: PetscCall(MatSetValues(mat, 1, &row, 1, &col, &val, ADD_VALUES));
185: row = r2 + l2;
186: col = r1 + l1;
187: PetscCall(MatSetValues(mat, 1, &row, 1, &col, &val, ADD_VALUES));
188: }
189: }
190: }
191: return PETSC_SUCCESS;
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: PetscCall(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: PetscCall(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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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*/