Actual source code: ex16.c
2: static char help[] = "Large-deformation Elasticity Buckling Example";
4: /*F-----------------------------------------------------------------------
6: This example solves the 3D large deformation elasticity problem
8: \begin{equation}
9: \int_{\Omega}F \cdot S : \nabla v d\Omega + \int_{\Omega} (loading)\mathbf{e}_y \cdot v d\Omega = 0
10: \end{equation}
12: F is the deformation gradient, and S is the second Piola-Kirchhoff tensor from the Saint Venant-Kirchhoff model of
13: hyperelasticity. \Omega is a (arc) angle subsection of a cylindrical shell of thickness (height), inner radius
14: (rad) and width (width). The problem is discretized using Q1 finite elements on a logically structured grid.
15: Homogeneous Dirichlet boundary conditions are applied at the centers of the ends of the sphere.
17: This example is tunable with the following options:
18: -rad : the radius of the circle
19: -arc : set the angle of the arch represented
20: -loading : set the bulk loading (the mass)
21: -ploading : set the point loading at the top
22: -height : set the height of the arch
23: -width : set the width of the arch
24: -view_line : print initial and final offsets of the centerline of the
25: beam along the x direction
27: The material properties may be modified using either:
28: -mu : the first lame' parameter
29: -lambda : the second lame' parameter
31: Or:
32: -young : the Young's modulus
33: -poisson : the poisson ratio
35: This example is meant to show the strain placed upon the nonlinear solvers when trying to "snap through" the arch
36: using the loading. Under certain parameter regimes, the arch will invert under the load, and the number of Newton
37: steps will jump considerably. Composed nonlinear solvers may be used to mitigate this difficulty.
39: The initial setup follows the example in pg. 268 of "Nonlinear Finite Element Methods" by Peter Wriggers, but is a
40: 3D extension.
42: ------------------------------------------------------------------------F*/
44: #include <petscsnes.h>
45: #include <petscdm.h>
46: #include <petscdmda.h>
48: #define QP0 0.2113248654051871
49: #define QP1 0.7886751345948129
50: #define NQ 2
51: #define NB 2
52: #define NEB 8
53: #define NEQ 8
54: #define NPB 24
56: #define NVALS NEB *NEQ
57: const PetscReal pts[NQ] = {QP0, QP1};
58: const PetscReal wts[NQ] = {0.5, 0.5};
60: PetscScalar vals[NVALS];
61: PetscScalar grad[3 * NVALS];
63: typedef PetscScalar Field[3];
64: typedef PetscScalar CoordField[3];
66: typedef PetscScalar JacField[9];
68: PetscErrorCode FormFunctionLocal(DMDALocalInfo *, Field ***, Field ***, void *);
69: PetscErrorCode FormJacobianLocal(DMDALocalInfo *, Field ***, Mat, Mat, void *);
70: PetscErrorCode DisplayLine(SNES, Vec);
71: PetscErrorCode FormElements(void);
73: typedef struct {
74: PetscReal loading;
75: PetscReal mu;
76: PetscReal lambda;
77: PetscReal rad;
78: PetscReal height;
79: PetscReal width;
80: PetscReal arc;
81: PetscReal ploading;
82: } AppCtx;
84: PetscErrorCode InitialGuess(DM, AppCtx *, Vec);
85: PetscErrorCode FormRHS(DM, AppCtx *, Vec);
86: PetscErrorCode FormCoordinates(DM, AppCtx *);
87: extern PetscErrorCode NonlinearGS(SNES, Vec, Vec, void *);
89: int main(int argc, char **argv)
90: {
91: AppCtx user; /* user-defined work context */
92: PetscInt mx, my, its;
93: MPI_Comm comm;
94: SNES snes;
95: DM da;
96: Vec x, X, b;
97: PetscBool youngflg, poissonflg, muflg, lambdaflg, view = PETSC_FALSE, viewline = PETSC_FALSE;
98: PetscReal poisson = 0.2, young = 4e4;
99: char filename[PETSC_MAX_PATH_LEN] = "ex16.vts";
100: char filename_def[PETSC_MAX_PATH_LEN] = "ex16_def.vts";
103: PetscInitialize(&argc, &argv, (char *)0, help);
104: FormElements();
105: comm = PETSC_COMM_WORLD;
106: SNESCreate(comm, &snes);
107: DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 11, 2, 2, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, NULL, NULL, &da);
108: DMSetFromOptions(da);
109: DMSetUp(da);
110: SNESSetDM(snes, (DM)da);
112: SNESSetNGS(snes, NonlinearGS, &user);
114: DMDAGetInfo(da, 0, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE);
115: user.loading = 0.0;
116: user.arc = PETSC_PI / 3.;
117: user.mu = 4.0;
118: user.lambda = 1.0;
119: user.rad = 100.0;
120: user.height = 3.;
121: user.width = 1.;
122: user.ploading = -5e3;
124: PetscOptionsGetReal(NULL, NULL, "-arc", &user.arc, NULL);
125: PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, &muflg);
126: PetscOptionsGetReal(NULL, NULL, "-lambda", &user.lambda, &lambdaflg);
127: PetscOptionsGetReal(NULL, NULL, "-rad", &user.rad, NULL);
128: PetscOptionsGetReal(NULL, NULL, "-height", &user.height, NULL);
129: PetscOptionsGetReal(NULL, NULL, "-width", &user.width, NULL);
130: PetscOptionsGetReal(NULL, NULL, "-loading", &user.loading, NULL);
131: PetscOptionsGetReal(NULL, NULL, "-ploading", &user.ploading, NULL);
132: PetscOptionsGetReal(NULL, NULL, "-poisson", &poisson, &poissonflg);
133: PetscOptionsGetReal(NULL, NULL, "-young", &young, &youngflg);
134: if ((youngflg || poissonflg) || !(muflg || lambdaflg)) {
135: /* set the lame' parameters based upon the poisson ratio and young's modulus */
136: user.lambda = poisson * young / ((1. + poisson) * (1. - 2. * poisson));
137: user.mu = young / (2. * (1. + poisson));
138: }
139: PetscOptionsGetBool(NULL, NULL, "-view", &view, NULL);
140: PetscOptionsGetBool(NULL, NULL, "-view_line", &viewline, NULL);
142: DMDASetFieldName(da, 0, "x_disp");
143: DMDASetFieldName(da, 1, "y_disp");
144: DMDASetFieldName(da, 2, "z_disp");
146: DMSetApplicationContext(da, &user);
147: DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode(*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user);
148: DMDASNESSetJacobianLocal(da, (DMDASNESJacobian)FormJacobianLocal, &user);
149: SNESSetFromOptions(snes);
150: FormCoordinates(da, &user);
152: DMCreateGlobalVector(da, &x);
153: DMCreateGlobalVector(da, &b);
154: InitialGuess(da, &user, x);
155: FormRHS(da, &user, b);
157: PetscPrintf(comm, "lambda: %f mu: %f\n", (double)user.lambda, (double)user.mu);
159: /* show a cross-section of the initial state */
160: if (viewline) DisplayLine(snes, x);
162: /* get the loaded configuration */
163: SNESSolve(snes, b, x);
165: SNESGetIterationNumber(snes, &its);
166: PetscPrintf(comm, "Number of SNES iterations = %" PetscInt_FMT "\n", its);
167: SNESGetSolution(snes, &X);
168: /* show a cross-section of the final state */
169: if (viewline) DisplayLine(snes, X);
171: if (view) {
172: PetscViewer viewer;
173: Vec coords;
174: PetscViewerVTKOpen(PETSC_COMM_WORLD, filename, FILE_MODE_WRITE, &viewer);
175: VecView(x, viewer);
176: PetscViewerDestroy(&viewer);
177: DMGetCoordinates(da, &coords);
178: VecAXPY(coords, 1.0, x);
179: PetscViewerVTKOpen(PETSC_COMM_WORLD, filename_def, FILE_MODE_WRITE, &viewer);
180: VecView(x, viewer);
181: PetscViewerDestroy(&viewer);
182: }
184: VecDestroy(&x);
185: VecDestroy(&b);
186: DMDestroy(&da);
187: SNESDestroy(&snes);
188: PetscFinalize();
189: return 0;
190: }
192: PetscInt OnBoundary(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz)
193: {
194: if ((i == 0 || i == mx - 1) && j == my - 1) return 1;
195: return 0;
196: }
198: void BoundaryValue(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz, PetscScalar *val, AppCtx *user)
199: {
200: /* reference coordinates */
201: PetscReal p_x = ((PetscReal)i) / (((PetscReal)(mx - 1)));
202: PetscReal p_y = ((PetscReal)j) / (((PetscReal)(my - 1)));
203: PetscReal p_z = ((PetscReal)k) / (((PetscReal)(mz - 1)));
204: PetscReal o_x = p_x;
205: PetscReal o_y = p_y;
206: PetscReal o_z = p_z;
207: val[0] = o_x - p_x;
208: val[1] = o_y - p_y;
209: val[2] = o_z - p_z;
210: }
212: void InvertTensor(PetscScalar *t, PetscScalar *ti, PetscReal *dett)
213: {
214: const PetscScalar a = t[0];
215: const PetscScalar b = t[1];
216: const PetscScalar c = t[2];
217: const PetscScalar d = t[3];
218: const PetscScalar e = t[4];
219: const PetscScalar f = t[5];
220: const PetscScalar g = t[6];
221: const PetscScalar h = t[7];
222: const PetscScalar i = t[8];
223: const PetscReal det = PetscRealPart(a * (e * i - f * h) - b * (i * d - f * g) + c * (d * h - e * g));
224: const PetscReal di = 1. / det;
225: if (ti) {
226: const PetscScalar A = (e * i - f * h);
227: const PetscScalar B = -(d * i - f * g);
228: const PetscScalar C = (d * h - e * g);
229: const PetscScalar D = -(b * i - c * h);
230: const PetscScalar E = (a * i - c * g);
231: const PetscScalar F = -(a * h - b * g);
232: const PetscScalar G = (b * f - c * e);
233: const PetscScalar H = -(a * f - c * d);
234: const PetscScalar II = (a * e - b * d);
235: ti[0] = di * A;
236: ti[1] = di * D;
237: ti[2] = di * G;
238: ti[3] = di * B;
239: ti[4] = di * E;
240: ti[5] = di * H;
241: ti[6] = di * C;
242: ti[7] = di * F;
243: ti[8] = di * II;
244: }
245: if (dett) *dett = det;
246: }
248: void TensorTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
249: {
250: PetscInt i, j, m;
251: for (i = 0; i < 3; i++) {
252: for (j = 0; j < 3; j++) {
253: c[i + 3 * j] = 0;
254: for (m = 0; m < 3; m++) c[i + 3 * j] += a[m + 3 * j] * b[i + 3 * m];
255: }
256: }
257: }
259: void TensorTransposeTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
260: {
261: PetscInt i, j, m;
262: for (i = 0; i < 3; i++) {
263: for (j = 0; j < 3; j++) {
264: c[i + 3 * j] = 0;
265: for (m = 0; m < 3; m++) c[i + 3 * j] += a[3 * m + j] * b[i + 3 * m];
266: }
267: }
268: }
270: void TensorVector(PetscScalar *rot, PetscScalar *vec, PetscScalar *tvec)
271: {
272: tvec[0] = rot[0] * vec[0] + rot[1] * vec[1] + rot[2] * vec[2];
273: tvec[1] = rot[3] * vec[0] + rot[4] * vec[1] + rot[5] * vec[2];
274: tvec[2] = rot[6] * vec[0] + rot[7] * vec[1] + rot[8] * vec[2];
275: }
277: void DeformationGradient(Field *ex, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *invJ, PetscScalar *F)
278: {
279: PetscInt ii, jj, kk, l;
280: for (l = 0; l < 9; l++) F[l] = 0.;
281: F[0] = 1.;
282: F[4] = 1.;
283: F[8] = 1.;
284: /* form the deformation gradient at this basis function -- loop over element unknowns */
285: for (kk = 0; kk < NB; kk++) {
286: for (jj = 0; jj < NB; jj++) {
287: for (ii = 0; ii < NB; ii++) {
288: PetscInt idx = ii + jj * NB + kk * NB * NB;
289: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
290: PetscScalar lgrad[3];
291: TensorVector(invJ, &grad[3 * bidx], lgrad);
292: F[0] += lgrad[0] * ex[idx][0];
293: F[1] += lgrad[1] * ex[idx][0];
294: F[2] += lgrad[2] * ex[idx][0];
295: F[3] += lgrad[0] * ex[idx][1];
296: F[4] += lgrad[1] * ex[idx][1];
297: F[5] += lgrad[2] * ex[idx][1];
298: F[6] += lgrad[0] * ex[idx][2];
299: F[7] += lgrad[1] * ex[idx][2];
300: F[8] += lgrad[2] * ex[idx][2];
301: }
302: }
303: }
304: }
306: void DeformationGradientJacobian(PetscInt qi, PetscInt qj, PetscInt qk, PetscInt ii, PetscInt jj, PetscInt kk, PetscInt fld, PetscScalar *invJ, PetscScalar *dF)
307: {
308: PetscInt l;
309: PetscScalar lgrad[3];
310: PetscInt idx = ii + jj * NB + kk * NB * NB;
311: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
312: for (l = 0; l < 9; l++) dF[l] = 0.;
313: /* form the deformation gradient at this basis function -- loop over element unknowns */
314: TensorVector(invJ, &grad[3 * bidx], lgrad);
315: dF[3 * fld] = lgrad[0];
316: dF[3 * fld + 1] = lgrad[1];
317: dF[3 * fld + 2] = lgrad[2];
318: }
320: void LagrangeGreenStrain(PetscScalar *F, PetscScalar *E)
321: {
322: PetscInt i, j, m;
323: for (i = 0; i < 3; i++) {
324: for (j = 0; j < 3; j++) {
325: E[i + 3 * j] = 0;
326: for (m = 0; m < 3; m++) E[i + 3 * j] += 0.5 * F[3 * m + j] * F[i + 3 * m];
327: }
328: }
329: for (i = 0; i < 3; i++) E[i + 3 * i] -= 0.5;
330: }
332: void SaintVenantKirchoff(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *S)
333: {
334: PetscInt i, j;
335: PetscScalar E[9];
336: PetscScalar trE = 0;
337: LagrangeGreenStrain(F, E);
338: for (i = 0; i < 3; i++) trE += E[i + 3 * i];
339: for (i = 0; i < 3; i++) {
340: for (j = 0; j < 3; j++) {
341: S[i + 3 * j] = 2. * mu * E[i + 3 * j];
342: if (i == j) S[i + 3 * j] += trE * lambda;
343: }
344: }
345: }
347: void SaintVenantKirchoffJacobian(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *dF, PetscScalar *dS)
348: {
349: PetscScalar FtdF[9], dE[9];
350: PetscInt i, j;
351: PetscScalar dtrE = 0.;
352: TensorTransposeTensor(dF, F, dE);
353: TensorTransposeTensor(F, dF, FtdF);
354: for (i = 0; i < 9; i++) dE[i] += FtdF[i];
355: for (i = 0; i < 9; i++) dE[i] *= 0.5;
356: for (i = 0; i < 3; i++) dtrE += dE[i + 3 * i];
357: for (i = 0; i < 3; i++) {
358: for (j = 0; j < 3; j++) {
359: dS[i + 3 * j] = 2. * mu * dE[i + 3 * j];
360: if (i == j) dS[i + 3 * j] += dtrE * lambda;
361: }
362: }
363: }
365: PetscErrorCode FormElements()
366: {
367: PetscInt i, j, k, ii, jj, kk;
368: PetscReal bx, by, bz, dbx, dby, dbz;
370: /* construct the basis function values and derivatives */
371: for (k = 0; k < NB; k++) {
372: for (j = 0; j < NB; j++) {
373: for (i = 0; i < NB; i++) {
374: /* loop over the quadrature points */
375: for (kk = 0; kk < NQ; kk++) {
376: for (jj = 0; jj < NQ; jj++) {
377: for (ii = 0; ii < NQ; ii++) {
378: PetscInt idx = ii + NQ * jj + NQ * NQ * kk + NEQ * i + NEQ * NB * j + NEQ * NB * NB * k;
379: bx = pts[ii];
380: by = pts[jj];
381: bz = pts[kk];
382: dbx = 1.;
383: dby = 1.;
384: dbz = 1.;
385: if (i == 0) {
386: bx = 1. - bx;
387: dbx = -1;
388: }
389: if (j == 0) {
390: by = 1. - by;
391: dby = -1;
392: }
393: if (k == 0) {
394: bz = 1. - bz;
395: dbz = -1;
396: }
397: vals[idx] = bx * by * bz;
398: grad[3 * idx + 0] = dbx * by * bz;
399: grad[3 * idx + 1] = dby * bx * bz;
400: grad[3 * idx + 2] = dbz * bx * by;
401: }
402: }
403: }
404: }
405: }
406: }
407: return 0;
408: }
410: void GatherElementData(PetscInt mx, PetscInt my, PetscInt mz, Field ***x, CoordField ***c, PetscInt i, PetscInt j, PetscInt k, Field *ex, CoordField *ec, AppCtx *user)
411: {
412: PetscInt m;
413: PetscInt ii, jj, kk;
414: /* gather the data -- loop over element unknowns */
415: for (kk = 0; kk < NB; kk++) {
416: for (jj = 0; jj < NB; jj++) {
417: for (ii = 0; ii < NB; ii++) {
418: PetscInt idx = ii + jj * NB + kk * NB * NB;
419: /* decouple the boundary nodes for the displacement variables */
420: if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
421: BoundaryValue(i + ii, j + jj, k + kk, mx, my, mz, ex[idx], user);
422: } else {
423: for (m = 0; m < 3; m++) ex[idx][m] = x[k + kk][j + jj][i + ii][m];
424: }
425: for (m = 0; m < 3; m++) ec[idx][m] = c[k + kk][j + jj][i + ii][m];
426: }
427: }
428: }
429: }
431: void QuadraturePointGeometricJacobian(CoordField *ec, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *J)
432: {
433: PetscInt ii, jj, kk;
434: /* construct the gradient at the given quadrature point named by i,j,k */
435: for (ii = 0; ii < 9; ii++) J[ii] = 0;
436: for (kk = 0; kk < NB; kk++) {
437: for (jj = 0; jj < NB; jj++) {
438: for (ii = 0; ii < NB; ii++) {
439: PetscInt idx = ii + jj * NB + kk * NB * NB;
440: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
441: J[0] += grad[3 * bidx + 0] * ec[idx][0];
442: J[1] += grad[3 * bidx + 1] * ec[idx][0];
443: J[2] += grad[3 * bidx + 2] * ec[idx][0];
444: J[3] += grad[3 * bidx + 0] * ec[idx][1];
445: J[4] += grad[3 * bidx + 1] * ec[idx][1];
446: J[5] += grad[3 * bidx + 2] * ec[idx][1];
447: J[6] += grad[3 * bidx + 0] * ec[idx][2];
448: J[7] += grad[3 * bidx + 1] * ec[idx][2];
449: J[8] += grad[3 * bidx + 2] * ec[idx][2];
450: }
451: }
452: }
453: }
455: void FormElementJacobian(Field *ex, CoordField *ec, Field *ef, PetscScalar *ej, AppCtx *user)
456: {
457: PetscReal vol;
458: PetscScalar J[9];
459: PetscScalar invJ[9];
460: PetscScalar F[9], S[9], dF[9], dS[9], dFS[9], FdS[9], FS[9];
461: PetscReal scl;
462: PetscInt i, j, k, l, ii, jj, kk, ll, qi, qj, qk, m;
464: if (ej)
465: for (i = 0; i < NPB * NPB; i++) ej[i] = 0.;
466: if (ef)
467: for (i = 0; i < NEB; i++) {
468: ef[i][0] = 0.;
469: ef[i][1] = 0.;
470: ef[i][2] = 0.;
471: }
472: /* loop over quadrature */
473: for (qk = 0; qk < NQ; qk++) {
474: for (qj = 0; qj < NQ; qj++) {
475: for (qi = 0; qi < NQ; qi++) {
476: QuadraturePointGeometricJacobian(ec, qi, qj, qk, J);
477: InvertTensor(J, invJ, &vol);
478: scl = vol * wts[qi] * wts[qj] * wts[qk];
479: DeformationGradient(ex, qi, qj, qk, invJ, F);
480: SaintVenantKirchoff(user->lambda, user->mu, F, S);
481: /* form the function */
482: if (ef) {
483: TensorTensor(F, S, FS);
484: for (kk = 0; kk < NB; kk++) {
485: for (jj = 0; jj < NB; jj++) {
486: for (ii = 0; ii < NB; ii++) {
487: PetscInt idx = ii + jj * NB + kk * NB * NB;
488: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
489: PetscScalar lgrad[3];
490: TensorVector(invJ, &grad[3 * bidx], lgrad);
491: /* mu*F : grad phi_{u,v,w} */
492: for (m = 0; m < 3; m++) ef[idx][m] += scl * (lgrad[0] * FS[3 * m + 0] + lgrad[1] * FS[3 * m + 1] + lgrad[2] * FS[3 * m + 2]);
493: ef[idx][1] -= scl * user->loading * vals[bidx];
494: }
495: }
496: }
497: }
498: /* form the jacobian */
499: if (ej) {
500: /* loop over trialfunctions */
501: for (k = 0; k < NB; k++) {
502: for (j = 0; j < NB; j++) {
503: for (i = 0; i < NB; i++) {
504: for (l = 0; l < 3; l++) {
505: PetscInt tridx = l + 3 * (i + j * NB + k * NB * NB);
506: DeformationGradientJacobian(qi, qj, qk, i, j, k, l, invJ, dF);
507: SaintVenantKirchoffJacobian(user->lambda, user->mu, F, dF, dS);
508: TensorTensor(dF, S, dFS);
509: TensorTensor(F, dS, FdS);
510: for (m = 0; m < 9; m++) dFS[m] += FdS[m];
511: /* loop over testfunctions */
512: for (kk = 0; kk < NB; kk++) {
513: for (jj = 0; jj < NB; jj++) {
514: for (ii = 0; ii < NB; ii++) {
515: PetscInt idx = ii + jj * NB + kk * NB * NB;
516: PetscInt bidx = 8 * idx + qi + NQ * qj + NQ * NQ * qk;
517: PetscScalar lgrad[3];
518: TensorVector(invJ, &grad[3 * bidx], lgrad);
519: for (ll = 0; ll < 3; ll++) {
520: PetscInt teidx = ll + 3 * (ii + jj * NB + kk * NB * NB);
521: ej[teidx + NPB * tridx] += scl * (lgrad[0] * dFS[3 * ll + 0] + lgrad[1] * dFS[3 * ll + 1] + lgrad[2] * dFS[3 * ll + 2]);
522: }
523: }
524: }
525: } /* end of testfunctions */
526: }
527: }
528: }
529: } /* end of trialfunctions */
530: }
531: }
532: }
533: } /* end of quadrature points */
534: }
536: void FormPBJacobian(PetscInt i, PetscInt j, PetscInt k, Field *ex, CoordField *ec, Field *ef, PetscScalar *ej, AppCtx *user)
537: {
538: PetscReal vol;
539: PetscScalar J[9];
540: PetscScalar invJ[9];
541: PetscScalar F[9], S[9], dF[9], dS[9], dFS[9], FdS[9], FS[9];
542: PetscReal scl;
543: PetscInt l, ll, qi, qj, qk, m;
544: PetscInt idx = i + j * NB + k * NB * NB;
545: PetscScalar lgrad[3];
547: if (ej)
548: for (l = 0; l < 9; l++) ej[l] = 0.;
549: if (ef)
550: for (l = 0; l < 1; l++) {
551: ef[l][0] = 0.;
552: ef[l][1] = 0.;
553: ef[l][2] = 0.;
554: }
555: /* loop over quadrature */
556: for (qk = 0; qk < NQ; qk++) {
557: for (qj = 0; qj < NQ; qj++) {
558: for (qi = 0; qi < NQ; qi++) {
559: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
560: QuadraturePointGeometricJacobian(ec, qi, qj, qk, J);
561: InvertTensor(J, invJ, &vol);
562: TensorVector(invJ, &grad[3 * bidx], lgrad);
563: scl = vol * wts[qi] * wts[qj] * wts[qk];
564: DeformationGradient(ex, qi, qj, qk, invJ, F);
565: SaintVenantKirchoff(user->lambda, user->mu, F, S);
566: /* form the function */
567: if (ef) {
568: TensorTensor(F, S, FS);
569: for (m = 0; m < 3; m++) ef[0][m] += scl * (lgrad[0] * FS[3 * m + 0] + lgrad[1] * FS[3 * m + 1] + lgrad[2] * FS[3 * m + 2]);
570: ef[0][1] -= scl * user->loading * vals[bidx];
571: }
572: /* form the jacobian */
573: if (ej) {
574: for (l = 0; l < 3; l++) {
575: DeformationGradientJacobian(qi, qj, qk, i, j, k, l, invJ, dF);
576: SaintVenantKirchoffJacobian(user->lambda, user->mu, F, dF, dS);
577: TensorTensor(dF, S, dFS);
578: TensorTensor(F, dS, FdS);
579: for (m = 0; m < 9; m++) dFS[m] += FdS[m];
580: for (ll = 0; ll < 3; ll++) ej[ll + 3 * l] += scl * (lgrad[0] * dFS[3 * ll + 0] + lgrad[1] * dFS[3 * ll + 1] + lgrad[2] * dFS[3 * ll + 2]);
581: }
582: }
583: }
584: } /* end of quadrature points */
585: }
586: }
588: void ApplyBCsElement(PetscInt mx, PetscInt my, PetscInt mz, PetscInt i, PetscInt j, PetscInt k, PetscScalar *jacobian)
589: {
590: PetscInt ii, jj, kk, ll, ei, ej, ek, el;
591: for (kk = 0; kk < NB; kk++) {
592: for (jj = 0; jj < NB; jj++) {
593: for (ii = 0; ii < NB; ii++) {
594: for (ll = 0; ll < 3; ll++) {
595: PetscInt tridx = ll + 3 * (ii + jj * NB + kk * NB * NB);
596: for (ek = 0; ek < NB; ek++) {
597: for (ej = 0; ej < NB; ej++) {
598: for (ei = 0; ei < NB; ei++) {
599: for (el = 0; el < 3; el++) {
600: if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz) || OnBoundary(i + ei, j + ej, k + ek, mx, my, mz)) {
601: PetscInt teidx = el + 3 * (ei + ej * NB + ek * NB * NB);
602: if (teidx == tridx) {
603: jacobian[tridx + NPB * teidx] = 1.;
604: } else {
605: jacobian[tridx + NPB * teidx] = 0.;
606: }
607: }
608: }
609: }
610: }
611: }
612: }
613: }
614: }
615: }
616: }
618: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field ***x, Mat jacpre, Mat jac, void *ptr)
619: {
620: /* values for each basis function at each quadrature point */
621: AppCtx *user = (AppCtx *)ptr;
622: PetscInt i, j, k, m, l;
623: PetscInt ii, jj, kk;
624: PetscScalar ej[NPB * NPB];
625: PetscScalar vals[NPB * NPB];
626: Field ex[NEB];
627: CoordField ec[NEB];
629: PetscInt xs = info->xs, ys = info->ys, zs = info->zs;
630: PetscInt xm = info->xm, ym = info->ym, zm = info->zm;
631: PetscInt xes, yes, zes, xee, yee, zee;
632: PetscInt mx = info->mx, my = info->my, mz = info->mz;
633: DM cda;
634: CoordField ***c;
635: Vec C;
636: PetscInt nrows;
637: MatStencil col[NPB], row[NPB];
638: PetscScalar v[9];
640: DMGetCoordinateDM(info->da, &cda);
641: DMGetCoordinatesLocal(info->da, &C);
642: DMDAVecGetArray(cda, C, &c);
643: MatScale(jac, 0.0);
645: xes = xs;
646: yes = ys;
647: zes = zs;
648: xee = xs + xm;
649: yee = ys + ym;
650: zee = zs + zm;
651: if (xs > 0) xes = xs - 1;
652: if (ys > 0) yes = ys - 1;
653: if (zs > 0) zes = zs - 1;
654: if (xs + xm == mx) xee = xs + xm - 1;
655: if (ys + ym == my) yee = ys + ym - 1;
656: if (zs + zm == mz) zee = zs + zm - 1;
657: for (k = zes; k < zee; k++) {
658: for (j = yes; j < yee; j++) {
659: for (i = xes; i < xee; i++) {
660: GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
661: FormElementJacobian(ex, ec, NULL, ej, user);
662: ApplyBCsElement(mx, my, mz, i, j, k, ej);
663: nrows = 0.;
664: for (kk = 0; kk < NB; kk++) {
665: for (jj = 0; jj < NB; jj++) {
666: for (ii = 0; ii < NB; ii++) {
667: PetscInt idx = ii + jj * 2 + kk * 4;
668: for (m = 0; m < 3; m++) {
669: col[3 * idx + m].i = i + ii;
670: col[3 * idx + m].j = j + jj;
671: col[3 * idx + m].k = k + kk;
672: col[3 * idx + m].c = m;
673: if (i + ii >= xs && i + ii < xm + xs && j + jj >= ys && j + jj < ys + ym && k + kk >= zs && k + kk < zs + zm) {
674: row[nrows].i = i + ii;
675: row[nrows].j = j + jj;
676: row[nrows].k = k + kk;
677: row[nrows].c = m;
678: for (l = 0; l < NPB; l++) vals[NPB * nrows + l] = ej[NPB * (3 * idx + m) + l];
679: nrows++;
680: }
681: }
682: }
683: }
684: }
685: MatSetValuesStencil(jac, nrows, row, NPB, col, vals, ADD_VALUES);
686: }
687: }
688: }
690: MatAssemblyBegin(jac, MAT_FLUSH_ASSEMBLY);
691: MatAssemblyEnd(jac, MAT_FLUSH_ASSEMBLY);
693: /* set the diagonal */
694: v[0] = 1.;
695: v[1] = 0.;
696: v[2] = 0.;
697: v[3] = 0.;
698: v[4] = 1.;
699: v[5] = 0.;
700: v[6] = 0.;
701: v[7] = 0.;
702: v[8] = 1.;
703: for (k = zs; k < zs + zm; k++) {
704: for (j = ys; j < ys + ym; j++) {
705: for (i = xs; i < xs + xm; i++) {
706: if (OnBoundary(i, j, k, mx, my, mz)) {
707: for (m = 0; m < 3; m++) {
708: col[m].i = i;
709: col[m].j = j;
710: col[m].k = k;
711: col[m].c = m;
712: }
713: MatSetValuesStencil(jac, 3, col, 3, col, v, INSERT_VALUES);
714: }
715: }
716: }
717: }
719: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
720: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
722: DMDAVecRestoreArray(cda, C, &c);
723: return 0;
724: }
726: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field ***x, Field ***f, void *ptr)
727: {
728: /* values for each basis function at each quadrature point */
729: AppCtx *user = (AppCtx *)ptr;
730: PetscInt i, j, k, l;
731: PetscInt ii, jj, kk;
733: Field ef[NEB];
734: Field ex[NEB];
735: CoordField ec[NEB];
737: PetscInt xs = info->xs, ys = info->ys, zs = info->zs;
738: PetscInt xm = info->xm, ym = info->ym, zm = info->zm;
739: PetscInt xes, yes, zes, xee, yee, zee;
740: PetscInt mx = info->mx, my = info->my, mz = info->mz;
741: DM cda;
742: CoordField ***c;
743: Vec C;
745: DMGetCoordinateDM(info->da, &cda);
746: DMGetCoordinatesLocal(info->da, &C);
747: DMDAVecGetArray(cda, C, &c);
748: DMDAGetInfo(info->da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
749: DMDAGetCorners(info->da, &xs, &ys, &zs, &xm, &ym, &zm);
751: /* loop over elements */
752: for (k = zs; k < zs + zm; k++) {
753: for (j = ys; j < ys + ym; j++) {
754: for (i = xs; i < xs + xm; i++) {
755: for (l = 0; l < 3; l++) f[k][j][i][l] = 0.;
756: }
757: }
758: }
759: /* element starts and ends */
760: xes = xs;
761: yes = ys;
762: zes = zs;
763: xee = xs + xm;
764: yee = ys + ym;
765: zee = zs + zm;
766: if (xs > 0) xes = xs - 1;
767: if (ys > 0) yes = ys - 1;
768: if (zs > 0) zes = zs - 1;
769: if (xs + xm == mx) xee = xs + xm - 1;
770: if (ys + ym == my) yee = ys + ym - 1;
771: if (zs + zm == mz) zee = zs + zm - 1;
772: for (k = zes; k < zee; k++) {
773: for (j = yes; j < yee; j++) {
774: for (i = xes; i < xee; i++) {
775: GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
776: FormElementJacobian(ex, ec, ef, NULL, user);
777: /* put this element's additions into the residuals */
778: for (kk = 0; kk < NB; kk++) {
779: for (jj = 0; jj < NB; jj++) {
780: for (ii = 0; ii < NB; ii++) {
781: PetscInt idx = ii + jj * NB + kk * NB * NB;
782: if (k + kk >= zs && j + jj >= ys && i + ii >= xs && k + kk < zs + zm && j + jj < ys + ym && i + ii < xs + xm) {
783: if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
784: for (l = 0; l < 3; l++) f[k + kk][j + jj][i + ii][l] = x[k + kk][j + jj][i + ii][l] - ex[idx][l];
785: } else {
786: for (l = 0; l < 3; l++) f[k + kk][j + jj][i + ii][l] += ef[idx][l];
787: }
788: }
789: }
790: }
791: }
792: }
793: }
794: }
795: DMDAVecRestoreArray(cda, C, &c);
796: return 0;
797: }
799: PetscErrorCode NonlinearGS(SNES snes, Vec X, Vec B, void *ptr)
800: {
801: /* values for each basis function at each quadrature point */
802: AppCtx *user = (AppCtx *)ptr;
803: PetscInt i, j, k, l, m, n, s;
804: PetscInt pi, pj, pk;
805: Field ef[1];
806: Field ex[8];
807: PetscScalar ej[9];
808: CoordField ec[8];
809: PetscScalar pjac[9], pjinv[9];
810: PetscScalar pf[3], py[3];
811: PetscInt xs, ys, zs;
812: PetscInt xm, ym, zm;
813: PetscInt mx, my, mz;
814: DM cda;
815: CoordField ***c;
816: Vec C;
817: DM da;
818: Vec Xl, Bl;
819: Field ***x, ***b;
820: PetscInt sweeps, its;
821: PetscReal atol, rtol, stol;
822: PetscReal fnorm0 = 0.0, fnorm, ynorm, xnorm = 0.0;
824: SNESNGSGetSweeps(snes, &sweeps);
825: SNESNGSGetTolerances(snes, &atol, &rtol, &stol, &its);
827: SNESGetDM(snes, &da);
828: DMGetLocalVector(da, &Xl);
829: if (B) DMGetLocalVector(da, &Bl);
830: DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xl);
831: DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xl);
832: if (B) {
833: DMGlobalToLocalBegin(da, B, INSERT_VALUES, Bl);
834: DMGlobalToLocalEnd(da, B, INSERT_VALUES, Bl);
835: }
836: DMDAVecGetArray(da, Xl, &x);
837: if (B) DMDAVecGetArray(da, Bl, &b);
839: DMGetCoordinateDM(da, &cda);
840: DMGetCoordinatesLocal(da, &C);
841: DMDAVecGetArray(cda, C, &c);
842: DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
843: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
845: for (s = 0; s < sweeps; s++) {
846: for (k = zs; k < zs + zm; k++) {
847: for (j = ys; j < ys + ym; j++) {
848: for (i = xs; i < xs + xm; i++) {
849: if (OnBoundary(i, j, k, mx, my, mz)) {
850: BoundaryValue(i, j, k, mx, my, mz, x[k][j][i], user);
851: } else {
852: for (n = 0; n < its; n++) {
853: for (m = 0; m < 9; m++) pjac[m] = 0.;
854: for (m = 0; m < 3; m++) pf[m] = 0.;
855: /* gather the elements for this point */
856: for (pk = -1; pk < 1; pk++) {
857: for (pj = -1; pj < 1; pj++) {
858: for (pi = -1; pi < 1; pi++) {
859: /* check that this element exists */
860: if (i + pi >= 0 && i + pi < mx - 1 && j + pj >= 0 && j + pj < my - 1 && k + pk >= 0 && k + pk < mz - 1) {
861: /* create the element function and jacobian */
862: GatherElementData(mx, my, mz, x, c, i + pi, j + pj, k + pk, ex, ec, user);
863: FormPBJacobian(-pi, -pj, -pk, ex, ec, ef, ej, user);
864: /* extract the point named by i,j,k from the whole element jacobian and function */
865: for (l = 0; l < 3; l++) {
866: pf[l] += ef[0][l];
867: for (m = 0; m < 3; m++) pjac[3 * m + l] += ej[3 * m + l];
868: }
869: }
870: }
871: }
872: }
873: /* invert */
874: InvertTensor(pjac, pjinv, NULL);
875: /* apply */
876: if (B)
877: for (m = 0; m < 3; m++) pf[m] -= b[k][j][i][m];
878: TensorVector(pjinv, pf, py);
879: xnorm = 0.;
880: for (m = 0; m < 3; m++) {
881: x[k][j][i][m] -= py[m];
882: xnorm += PetscRealPart(x[k][j][i][m] * x[k][j][i][m]);
883: }
884: fnorm = PetscRealPart(pf[0] * pf[0] + pf[1] * pf[1] + pf[2] * pf[2]);
885: if (n == 0) fnorm0 = fnorm;
886: ynorm = PetscRealPart(py[0] * py[0] + py[1] * py[1] + py[2] * py[2]);
887: if (fnorm < atol * atol || fnorm < rtol * rtol * fnorm0 || ynorm < stol * stol * xnorm) break;
888: }
889: }
890: }
891: }
892: }
893: }
894: DMDAVecRestoreArray(da, Xl, &x);
895: DMLocalToGlobalBegin(da, Xl, INSERT_VALUES, X);
896: DMLocalToGlobalEnd(da, Xl, INSERT_VALUES, X);
897: DMRestoreLocalVector(da, &Xl);
898: if (B) {
899: DMDAVecRestoreArray(da, Bl, &b);
900: DMRestoreLocalVector(da, &Bl);
901: }
902: DMDAVecRestoreArray(cda, C, &c);
903: return 0;
904: }
906: PetscErrorCode FormCoordinates(DM da, AppCtx *user)
907: {
908: Vec coords;
909: DM cda;
910: PetscInt mx, my, mz;
911: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
912: CoordField ***x;
914: DMGetCoordinateDM(da, &cda);
915: DMCreateGlobalVector(cda, &coords);
916: DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
917: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
918: DMDAVecGetArray(da, coords, &x);
919: for (k = zs; k < zs + zm; k++) {
920: for (j = ys; j < ys + ym; j++) {
921: for (i = xs; i < xs + xm; i++) {
922: PetscReal cx = ((PetscReal)i) / (((PetscReal)(mx - 1)));
923: PetscReal cy = ((PetscReal)j) / (((PetscReal)(my - 1)));
924: PetscReal cz = ((PetscReal)k) / (((PetscReal)(mz - 1)));
925: PetscReal rad = user->rad + cy * user->height;
926: PetscReal ang = (cx - 0.5) * user->arc;
927: x[k][j][i][0] = rad * PetscSinReal(ang);
928: x[k][j][i][1] = rad * PetscCosReal(ang) - (user->rad + 0.5 * user->height) * PetscCosReal(-0.5 * user->arc);
929: x[k][j][i][2] = user->width * (cz - 0.5);
930: }
931: }
932: }
933: DMDAVecRestoreArray(da, coords, &x);
934: DMSetCoordinates(da, coords);
935: VecDestroy(&coords);
936: return 0;
937: }
939: PetscErrorCode InitialGuess(DM da, AppCtx *user, Vec X)
940: {
941: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
942: PetscInt mx, my, mz;
943: Field ***x;
945: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
946: DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
947: DMDAVecGetArray(da, X, &x);
949: for (k = zs; k < zs + zm; k++) {
950: for (j = ys; j < ys + ym; j++) {
951: for (i = xs; i < xs + xm; i++) {
952: /* reference coordinates */
953: PetscReal p_x = ((PetscReal)i) / (((PetscReal)(mx - 1)));
954: PetscReal p_y = ((PetscReal)j) / (((PetscReal)(my - 1)));
955: PetscReal p_z = ((PetscReal)k) / (((PetscReal)(mz - 1)));
956: PetscReal o_x = p_x;
957: PetscReal o_y = p_y;
958: PetscReal o_z = p_z;
959: x[k][j][i][0] = o_x - p_x;
960: x[k][j][i][1] = o_y - p_y;
961: x[k][j][i][2] = o_z - p_z;
962: }
963: }
964: }
965: DMDAVecRestoreArray(da, X, &x);
966: return 0;
967: }
969: PetscErrorCode FormRHS(DM da, AppCtx *user, Vec X)
970: {
971: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
972: PetscInt mx, my, mz;
973: Field ***x;
975: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
976: DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
977: DMDAVecGetArray(da, X, &x);
979: for (k = zs; k < zs + zm; k++) {
980: for (j = ys; j < ys + ym; j++) {
981: for (i = xs; i < xs + xm; i++) {
982: x[k][j][i][0] = 0.;
983: x[k][j][i][1] = 0.;
984: x[k][j][i][2] = 0.;
985: if (i == (mx - 1) / 2 && j == (my - 1)) x[k][j][i][1] = user->ploading / (mz - 1);
986: }
987: }
988: }
989: DMDAVecRestoreArray(da, X, &x);
990: return 0;
991: }
993: PetscErrorCode DisplayLine(SNES snes, Vec X)
994: {
995: PetscInt r, i, j = 0, k = 0, xs, xm, ys, ym, zs, zm, mx, my, mz;
996: Field ***x;
997: CoordField ***c;
998: DM da, cda;
999: Vec C;
1000: PetscMPIInt size, rank;
1002: MPI_Comm_size(PETSC_COMM_WORLD, &size);
1003: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1004: SNESGetDM(snes, &da);
1005: DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
1006: DMGetCoordinateDM(da, &cda);
1007: DMGetCoordinates(da, &C);
1008: j = my / 2;
1009: k = mz / 2;
1010: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
1011: DMDAVecGetArray(da, X, &x);
1012: DMDAVecGetArray(cda, C, &c);
1013: for (r = 0; r < size; r++) {
1014: if (rank == r) {
1015: if (j >= ys && j < ys + ym && k >= zs && k < zs + zm) {
1016: for (i = xs; i < xs + xm; i++) {
1017: PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %d %d: %f %f %f\n", i, 0, 0, (double)PetscRealPart(c[k][j][i][0] + x[k][j][i][0]), (double)PetscRealPart(c[k][j][i][1] + x[k][j][i][1]), (double)PetscRealPart(c[k][j][i][2] + x[k][j][i][2]));
1018: }
1019: }
1020: }
1021: MPI_Barrier(PETSC_COMM_WORLD);
1022: }
1023: DMDAVecRestoreArray(da, X, &x);
1024: DMDAVecRestoreArray(cda, C, &c);
1025: return 0;
1026: }
1028: /*TEST
1030: test:
1031: nsize: 2
1032: args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading 0.0 -loading -1. -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -snes_max_it 7
1033: requires: !single
1034: timeoutfactor: 3
1036: test:
1037: suffix: 2
1038: args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading 0.0 -loading -1. -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -npc_snes_type fas -npc_fas_levels_snes_type ncg -npc_fas_levels_snes_max_it 3 -npc_snes_monitor_short -snes_max_it 2
1039: requires: !single
1041: test:
1042: suffix: 3
1043: args: -da_refine 1 -da_overlap 3 -da_local_subdomains 4 -snes_type aspin -rad 10.0 -young 10. -ploading 0.0 -loading -0.5 -snes_monitor_short -ksp_monitor_short -npc_sub_snes_rtol 1e-2 -ksp_rtol 1e-2 -ksp_max_it 14 -snes_converged_reason -snes_max_linear_solve_fail 100 -snes_max_it 4
1044: requires: !single
1046: TEST*/