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";
5
6 #include <petscksp.h>
7
8 /* This code is not intended as an efficient implementation, it is only
9 here to produce an interesting sparse matrix quickly.
10
11 PLEASE DO NOT BASE ANY OF YOUR CODES ON CODE LIKE THIS, THERE ARE MUCH
12 BETTER WAYS TO DO THIS. */
13
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]);
19
main(int argc,char ** args)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;
29
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));
35
36 /* Form matrix */
37 PetscCall(GetElasticityMatrix(m, &mat));
38
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));
53
54 /* Compute right-hand side */
55 PetscCall(MatMult(mat, u, b));
56
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));
66
67 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of residual %g Number of iterations %" PetscInt_FMT "\n", (double)norm, its));
68
69 /* Free work space */
70 PetscCall(KSPDestroy(&ksp));
71 PetscCall(VecDestroy(&u));
72 PetscCall(VecDestroy(&x));
73 PetscCall(VecDestroy(&b));
74 PetscCall(MatDestroy(&mat));
75
76 PetscCall(PetscFinalize());
77 return 0;
78 }
79 /* -------------------------------------------------------------------- */
80 /*
81 GetElasticityMatrix - Forms 3D linear elasticity matrix.
82 */
GetElasticityMatrix(PetscInt m,Mat * newmat)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;
91
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));
96
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));
101
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 }
132
133 for (i = 0; i < 81; i++) PetscCall(PetscFree(K[i]));
134 PetscCall(PetscFree(K));
135
136 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
137 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
138
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));
155
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));
161
162 PetscCall(MatNorm(*newmat, NORM_1, &norm));
163 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "matrix 1 norm = %g\n", (double)norm));
164
165 return PETSC_SUCCESS;
166 }
167 /* -------------------------------------------------------------------- */
AddElement(Mat mat,PetscInt r1,PetscInt r2,PetscReal ** K,PetscInt h1,PetscInt h2)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;
172
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 */
Elastic20Stiff(PetscReal ** Ke)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;
213
214 PetscCall(paulsetup20());
215
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));
283
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 }
299
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 */
paulsetup20(void)310 PetscErrorCode paulsetup20(void)
311 {
312 PetscInt i, j, k, cnt;
313 PetscReal x[4], w[4];
314 PetscReal c;
315
316 n_int = 3;
317 nu = 0.3;
318 E = 1.0;
319
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");
344
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. */
348
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;
362
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). */
367
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 */
paulintegrate20(PetscReal K[60][60])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;
406
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 }
411
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;
431
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 }
442
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 }
455
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];
474
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 }
482
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 }
494
495 /*TEST
496
497 testset:
498 args: -matconvert_type seqaij -ksp_monitor_short -ksp_rtol 1.e-2 -pc_type jacobi
499 requires: x
500 output_file: output/ex10_1.out
501
502 test:
503 suffix: 1
504
505 test:
506 # Test MatMult_SeqAIJ OpenMP kernels with or without inode
507 suffix: 2
508 args: -omp_num_threads 4 -mat_no_inode {{0 1}}
509 requires: PETSC_HAVE_OPENMP_KERNELS
510
511 TEST*/
512