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