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 ArcLengthScaling(DM, AppCtx *, Vec);
88: PetscErrorCode FormRHS(DM, AppCtx *, Vec);
89: PetscErrorCode FormCoordinates(DM, AppCtx *);
90: PetscErrorCode TangentLoad(SNES, Vec, Vec, void *);
92: int main(int argc, char **argv)
93: {
94: AppCtx user; /* user-defined work context */
95: PetscInt mx, my, its;
96: MPI_Comm comm;
97: SNES snes;
98: DM da;
99: Vec x, X, b;
100: PetscBool youngflg, poissonflg, muflg, lambdaflg, alflg, view = PETSC_FALSE, viewline = PETSC_FALSE, al_rescale = PETSC_FALSE;
101: PetscReal poisson = 0.2, young = 4e4;
102: char filename[PETSC_MAX_PATH_LEN] = "ex16.vts";
103: char filename_def[PETSC_MAX_PATH_LEN] = "ex16_def.vts";
105: PetscFunctionBeginUser;
106: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
107: PetscCall(FormElements());
108: comm = PETSC_COMM_WORLD;
109: PetscCall(SNESCreate(comm, &snes));
110: 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));
111: PetscCall(DMSetFromOptions(da));
112: PetscCall(DMSetUp(da));
113: PetscCall(SNESSetDM(snes, (DM)da));
115: 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));
116: user.loading = 0.0;
117: user.arc = PETSC_PI / 3.;
118: user.mu = 4.0;
119: user.lambda = 1.0;
120: user.rad = 100.0;
121: user.height = 3.;
122: user.width = 1.;
123: user.ploading = -5e3;
124: user.load_factor = 1.0;
126: PetscCall(PetscOptionsGetReal(NULL, NULL, "-arc", &user.arc, NULL));
127: PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, &muflg));
128: PetscCall(PetscOptionsGetReal(NULL, NULL, "-lambda", &user.lambda, &lambdaflg));
129: PetscCall(PetscOptionsGetReal(NULL, NULL, "-rad", &user.rad, NULL));
130: PetscCall(PetscOptionsGetReal(NULL, NULL, "-height", &user.height, NULL));
131: PetscCall(PetscOptionsGetReal(NULL, NULL, "-width", &user.width, NULL));
132: PetscCall(PetscOptionsGetReal(NULL, NULL, "-loading", &user.loading, NULL));
133: PetscCall(PetscOptionsGetReal(NULL, NULL, "-ploading", &user.ploading, NULL));
134: PetscCall(PetscOptionsGetReal(NULL, NULL, "-poisson", &poisson, &poissonflg));
135: PetscCall(PetscOptionsGetReal(NULL, NULL, "-young", &young, &youngflg));
136: if (youngflg || poissonflg || !(muflg || lambdaflg)) {
137: /* set the lame' parameters based upon the poisson ratio and young's modulus */
138: user.lambda = poisson * young / ((1. + poisson) * (1. - 2. * poisson));
139: user.mu = young / (2. * (1. + poisson));
140: }
141: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view", &view, NULL));
142: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_line", &viewline, NULL));
143: PetscCall(PetscOptionsGetBool(NULL, NULL, "-rescale", &al_rescale, NULL));
145: PetscCall(DMDASetFieldName(da, 0, "x_disp"));
146: PetscCall(DMDASetFieldName(da, 1, "y_disp"));
147: PetscCall(DMDASetFieldName(da, 2, "z_disp"));
149: PetscCall(DMSetApplicationContext(da, &user));
150: PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunctionFn *)FormFunctionLocal, &user));
151: PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)FormJacobianLocal, &user));
152: PetscCall(SNESSetFromOptions(snes));
153: PetscCall(SNESNewtonALSetFunction(snes, TangentLoad, &user));
154: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONAL, &alflg));
155: if (alflg) {
156: user.load_factor = 0.0;
158: if (al_rescale) {
159: Vec scaling;
161: PetscCall(DMCreateGlobalVector(da, &scaling));
162: PetscCall(ArcLengthScaling(da, &user, scaling));
163: PetscCall(VecPointwiseMult(scaling, scaling, scaling));
164: PetscCall(SNESNewtonALSetDiagonalScaling(snes, scaling));
165: PetscCall(VecDestroy(&scaling));
166: }
167: }
169: PetscCall(FormCoordinates(da, &user));
171: PetscCall(DMCreateGlobalVector(da, &x));
172: PetscCall(DMCreateGlobalVector(da, &b));
173: PetscCall(InitialGuess(da, &user, x));
174: PetscCall(FormRHS(da, &user, b));
176: PetscCall(PetscPrintf(comm, "lambda: %f mu: %f\n", (double)user.lambda, (double)user.mu));
178: /* show a cross-section of the initial state */
179: if (viewline) PetscCall(DisplayLine(snes, x));
181: /* get the loaded configuration */
182: PetscCall(SNESSolve(snes, b, x));
184: PetscCall(SNESGetIterationNumber(snes, &its));
185: PetscCall(PetscPrintf(comm, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
186: PetscCall(SNESGetSolution(snes, &X));
187: /* show a cross-section of the final state */
188: if (viewline) PetscCall(DisplayLine(snes, X));
190: if (view) {
191: PetscViewer viewer;
192: Vec coords;
193: PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename, FILE_MODE_WRITE, &viewer));
194: PetscCall(VecView(x, viewer));
195: PetscCall(PetscViewerDestroy(&viewer));
196: PetscCall(DMGetCoordinates(da, &coords));
197: PetscCall(VecAXPY(coords, 1.0, x));
198: PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename_def, FILE_MODE_WRITE, &viewer));
199: PetscCall(VecView(x, viewer));
200: PetscCall(PetscViewerDestroy(&viewer));
201: }
203: PetscCall(VecDestroy(&x));
204: PetscCall(VecDestroy(&b));
205: PetscCall(DMDestroy(&da));
206: PetscCall(SNESDestroy(&snes));
207: PetscCall(PetscFinalize());
208: return 0;
209: }
211: PetscInt OnBoundary(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz)
212: {
213: if ((i == 0 || i == mx - 1) && j == my - 1) return 1;
214: return 0;
215: }
217: void BoundaryValue(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz, PetscScalar *val, AppCtx *user)
218: {
219: /* reference coordinates */
220: PetscReal p_x = ((PetscReal)i) / ((PetscReal)(mx - 1));
221: PetscReal p_y = ((PetscReal)j) / ((PetscReal)(my - 1));
222: PetscReal p_z = ((PetscReal)k) / ((PetscReal)(mz - 1));
223: PetscReal o_x = p_x;
224: PetscReal o_y = p_y;
225: PetscReal o_z = p_z;
226: val[0] = o_x - p_x;
227: val[1] = o_y - p_y;
228: val[2] = o_z - p_z;
229: }
231: void InvertTensor(PetscScalar *t, PetscScalar *ti, PetscReal *dett)
232: {
233: const PetscScalar a = t[0];
234: const PetscScalar b = t[1];
235: const PetscScalar c = t[2];
236: const PetscScalar d = t[3];
237: const PetscScalar e = t[4];
238: const PetscScalar f = t[5];
239: const PetscScalar g = t[6];
240: const PetscScalar h = t[7];
241: const PetscScalar i = t[8];
242: const PetscReal det = PetscRealPart(a * (e * i - f * h) - b * (i * d - f * g) + c * (d * h - e * g));
243: const PetscReal di = 1. / det;
244: if (ti) {
245: const PetscScalar A = (e * i - f * h);
246: const PetscScalar B = -(d * i - f * g);
247: const PetscScalar C = (d * h - e * g);
248: const PetscScalar D = -(b * i - c * h);
249: const PetscScalar E = (a * i - c * g);
250: const PetscScalar F = -(a * h - b * g);
251: const PetscScalar G = (b * f - c * e);
252: const PetscScalar H = -(a * f - c * d);
253: const PetscScalar II = (a * e - b * d);
254: ti[0] = di * A;
255: ti[1] = di * D;
256: ti[2] = di * G;
257: ti[3] = di * B;
258: ti[4] = di * E;
259: ti[5] = di * H;
260: ti[6] = di * C;
261: ti[7] = di * F;
262: ti[8] = di * II;
263: }
264: if (dett) *dett = det;
265: }
267: void TensorTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
268: {
269: PetscInt i, j, m;
270: for (i = 0; i < 3; i++) {
271: for (j = 0; j < 3; j++) {
272: c[i + 3 * j] = 0;
273: for (m = 0; m < 3; m++) c[i + 3 * j] += a[m + 3 * j] * b[i + 3 * m];
274: }
275: }
276: }
278: void TensorTransposeTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
279: {
280: PetscInt i, j, m;
281: for (i = 0; i < 3; i++) {
282: for (j = 0; j < 3; j++) {
283: c[i + 3 * j] = 0;
284: for (m = 0; m < 3; m++) c[i + 3 * j] += a[3 * m + j] * b[i + 3 * m];
285: }
286: }
287: }
289: void TensorVector(PetscScalar *rot, PetscScalar *vec, PetscScalar *tvec)
290: {
291: tvec[0] = rot[0] * vec[0] + rot[1] * vec[1] + rot[2] * vec[2];
292: tvec[1] = rot[3] * vec[0] + rot[4] * vec[1] + rot[5] * vec[2];
293: tvec[2] = rot[6] * vec[0] + rot[7] * vec[1] + rot[8] * vec[2];
294: }
296: void DeformationGradient(Field *ex, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *invJ, PetscScalar *F)
297: {
298: PetscInt ii, jj, kk, l;
299: for (l = 0; l < 9; l++) F[l] = 0.;
300: F[0] = 1.;
301: F[4] = 1.;
302: F[8] = 1.;
303: /* form the deformation gradient at this basis function -- loop over element unknowns */
304: for (kk = 0; kk < NB; kk++) {
305: for (jj = 0; jj < NB; jj++) {
306: for (ii = 0; ii < NB; ii++) {
307: PetscInt idx = ii + jj * NB + kk * NB * NB;
308: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
309: PetscScalar lgrad[3];
310: TensorVector(invJ, &grad[3 * bidx], lgrad);
311: F[0] += lgrad[0] * ex[idx][0];
312: F[1] += lgrad[1] * ex[idx][0];
313: F[2] += lgrad[2] * ex[idx][0];
314: F[3] += lgrad[0] * ex[idx][1];
315: F[4] += lgrad[1] * ex[idx][1];
316: F[5] += lgrad[2] * ex[idx][1];
317: F[6] += lgrad[0] * ex[idx][2];
318: F[7] += lgrad[1] * ex[idx][2];
319: F[8] += lgrad[2] * ex[idx][2];
320: }
321: }
322: }
323: }
325: void DeformationGradientJacobian(PetscInt qi, PetscInt qj, PetscInt qk, PetscInt ii, PetscInt jj, PetscInt kk, PetscInt fld, PetscScalar *invJ, PetscScalar *dF)
326: {
327: PetscScalar lgrad[3];
328: PetscInt idx = ii + jj * NB + kk * NB * NB;
329: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
330: for (PetscInt l = 0; l < 9; l++) dF[l] = 0.;
331: /* form the deformation gradient at this basis function -- loop over element unknowns */
332: TensorVector(invJ, &grad[3 * bidx], lgrad);
333: dF[3 * fld] = lgrad[0];
334: dF[3 * fld + 1] = lgrad[1];
335: dF[3 * fld + 2] = lgrad[2];
336: }
338: void LagrangeGreenStrain(PetscScalar *F, PetscScalar *E)
339: {
340: PetscInt i, j, m;
341: for (i = 0; i < 3; i++) {
342: for (j = 0; j < 3; j++) {
343: E[i + 3 * j] = 0;
344: for (m = 0; m < 3; m++) E[i + 3 * j] += 0.5 * F[3 * m + j] * F[i + 3 * m];
345: }
346: }
347: for (i = 0; i < 3; i++) E[i + 3 * i] -= 0.5;
348: }
350: void SaintVenantKirchoff(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *S)
351: {
352: PetscInt i, j;
353: PetscScalar E[9];
354: PetscScalar trE = 0;
355: LagrangeGreenStrain(F, E);
356: for (i = 0; i < 3; i++) trE += E[i + 3 * i];
357: for (i = 0; i < 3; i++) {
358: for (j = 0; j < 3; j++) {
359: S[i + 3 * j] = 2. * mu * E[i + 3 * j];
360: if (i == j) S[i + 3 * j] += trE * lambda;
361: }
362: }
363: }
365: void SaintVenantKirchoffJacobian(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *dF, PetscScalar *dS)
366: {
367: PetscScalar FtdF[9], dE[9];
368: PetscInt i, j;
369: PetscScalar dtrE = 0.;
370: TensorTransposeTensor(dF, F, dE);
371: TensorTransposeTensor(F, dF, FtdF);
372: for (i = 0; i < 9; i++) dE[i] += FtdF[i];
373: for (i = 0; i < 9; i++) dE[i] *= 0.5;
374: for (i = 0; i < 3; i++) dtrE += dE[i + 3 * i];
375: for (i = 0; i < 3; i++) {
376: for (j = 0; j < 3; j++) {
377: dS[i + 3 * j] = 2. * mu * dE[i + 3 * j];
378: if (i == j) dS[i + 3 * j] += dtrE * lambda;
379: }
380: }
381: }
383: PetscErrorCode FormElements(void)
384: {
385: PetscInt i, j, k, ii, jj, kk;
386: PetscReal bx, by, bz, dbx, dby, dbz;
388: PetscFunctionBeginUser;
389: /* construct the basis function values and derivatives */
390: for (k = 0; k < NB; k++) {
391: for (j = 0; j < NB; j++) {
392: for (i = 0; i < NB; i++) {
393: /* loop over the quadrature points */
394: for (kk = 0; kk < NQ; kk++) {
395: for (jj = 0; jj < NQ; jj++) {
396: for (ii = 0; ii < NQ; ii++) {
397: PetscInt idx = ii + NQ * jj + NQ * NQ * kk + NEQ * i + NEQ * NB * j + NEQ * NB * NB * k;
398: bx = pts[ii];
399: by = pts[jj];
400: bz = pts[kk];
401: dbx = 1.;
402: dby = 1.;
403: dbz = 1.;
404: if (i == 0) {
405: bx = 1. - bx;
406: dbx = -1;
407: }
408: if (j == 0) {
409: by = 1. - by;
410: dby = -1;
411: }
412: if (k == 0) {
413: bz = 1. - bz;
414: dbz = -1;
415: }
416: vals[idx] = bx * by * bz;
417: grad[3 * idx + 0] = dbx * by * bz;
418: grad[3 * idx + 1] = dby * bx * bz;
419: grad[3 * idx + 2] = dbz * bx * by;
420: }
421: }
422: }
423: }
424: }
425: }
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: void GatherElementData(PetscInt mx, PetscInt my, PetscInt mz, Field ***x, CoordField ***c, PetscInt i, PetscInt j, PetscInt k, Field *ex, CoordField *ec, AppCtx *user)
430: {
431: PetscInt m;
432: PetscInt ii, jj, kk;
433: /* gather the data -- loop over element unknowns */
434: for (kk = 0; kk < NB; kk++) {
435: for (jj = 0; jj < NB; jj++) {
436: for (ii = 0; ii < NB; ii++) {
437: PetscInt idx = ii + jj * NB + kk * NB * NB;
438: /* decouple the boundary nodes for the displacement variables */
439: if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
440: BoundaryValue(i + ii, j + jj, k + kk, mx, my, mz, ex[idx], user);
441: } else {
442: for (m = 0; m < 3; m++) ex[idx][m] = x[k + kk][j + jj][i + ii][m];
443: }
444: for (m = 0; m < 3; m++) ec[idx][m] = c[k + kk][j + jj][i + ii][m];
445: }
446: }
447: }
448: }
450: void QuadraturePointGeometricJacobian(CoordField *ec, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *J)
451: {
452: PetscInt ii, jj, kk;
453: /* construct the gradient at the given quadrature point named by i,j,k */
454: for (ii = 0; ii < 9; ii++) J[ii] = 0;
455: for (kk = 0; kk < NB; kk++) {
456: for (jj = 0; jj < NB; jj++) {
457: for (ii = 0; ii < NB; ii++) {
458: PetscInt idx = ii + jj * NB + kk * NB * NB;
459: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
460: J[0] += grad[3 * bidx + 0] * ec[idx][0];
461: J[1] += grad[3 * bidx + 1] * ec[idx][0];
462: J[2] += grad[3 * bidx + 2] * ec[idx][0];
463: J[3] += grad[3 * bidx + 0] * ec[idx][1];
464: J[4] += grad[3 * bidx + 1] * ec[idx][1];
465: J[5] += grad[3 * bidx + 2] * ec[idx][1];
466: J[6] += grad[3 * bidx + 0] * ec[idx][2];
467: J[7] += grad[3 * bidx + 1] * ec[idx][2];
468: J[8] += grad[3 * bidx + 2] * ec[idx][2];
469: }
470: }
471: }
472: }
474: void FormElementJacobian(Field *ex, CoordField *ec, Field *ef, Field *eq, PetscScalar *ej, AppCtx *user)
475: {
476: PetscReal vol;
477: PetscScalar J[9];
478: PetscScalar invJ[9];
479: PetscScalar F[9], S[9], dF[9], dS[9], dFS[9], FdS[9], FS[9];
480: PetscReal scl;
481: PetscInt i, j, k, l, ii, jj, kk, ll, qi, qj, qk, m;
483: if (ej)
484: for (i = 0; i < NPB * NPB; i++) ej[i] = 0.;
485: if (ef)
486: for (i = 0; i < NEB; i++) {
487: ef[i][0] = 0.;
488: ef[i][1] = 0.;
489: ef[i][2] = 0.;
490: }
491: if (eq)
492: for (i = 0; i < NEB; i++) {
493: eq[i][0] = 0.;
494: eq[i][1] = 0.;
495: eq[i][2] = 0.;
496: }
497: /* loop over quadrature */
498: for (qk = 0; qk < NQ; qk++) {
499: for (qj = 0; qj < NQ; qj++) {
500: for (qi = 0; qi < NQ; qi++) {
501: QuadraturePointGeometricJacobian(ec, qi, qj, qk, J);
502: InvertTensor(J, invJ, &vol);
503: scl = vol * wts[qi] * wts[qj] * wts[qk];
504: DeformationGradient(ex, qi, qj, qk, invJ, F);
505: SaintVenantKirchoff(user->lambda, user->mu, F, S);
506: /* form the function */
507: if (ef) {
508: TensorTensor(F, S, FS);
509: for (kk = 0; kk < NB; kk++) {
510: for (jj = 0; jj < NB; jj++) {
511: for (ii = 0; ii < NB; ii++) {
512: PetscInt idx = ii + jj * NB + kk * NB * NB;
513: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
514: PetscScalar lgrad[3];
515: TensorVector(invJ, &grad[3 * bidx], lgrad);
516: /* mu*F : grad phi_{u,v,w} */
517: 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]);
518: ef[idx][1] -= user->load_factor * scl * user->loading * vals[bidx];
519: }
520: }
521: }
522: }
523: if (eq) {
524: for (kk = 0; kk < NB; kk++) {
525: for (jj = 0; jj < NB; jj++) {
526: for (ii = 0; ii < NB; ii++) {
527: PetscInt idx = ii + jj * NB + kk * NB * NB;
528: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
529: /* external force vector */
530: eq[idx][1] += scl * user->loading * vals[bidx];
531: }
532: }
533: }
534: }
535: /* form the jacobian */
536: if (ej) {
537: /* loop over trialfunctions */
538: for (k = 0; k < NB; k++) {
539: for (j = 0; j < NB; j++) {
540: for (i = 0; i < NB; i++) {
541: for (l = 0; l < 3; l++) {
542: PetscInt tridx = l + 3 * (i + j * NB + k * NB * NB);
543: DeformationGradientJacobian(qi, qj, qk, i, j, k, l, invJ, dF);
544: SaintVenantKirchoffJacobian(user->lambda, user->mu, F, dF, dS);
545: TensorTensor(dF, S, dFS);
546: TensorTensor(F, dS, FdS);
547: for (m = 0; m < 9; m++) dFS[m] += FdS[m];
548: /* loop over testfunctions */
549: for (kk = 0; kk < NB; kk++) {
550: for (jj = 0; jj < NB; jj++) {
551: for (ii = 0; ii < NB; ii++) {
552: PetscInt idx = ii + jj * NB + kk * NB * NB;
553: PetscInt bidx = 8 * idx + qi + NQ * qj + NQ * NQ * qk;
554: PetscScalar lgrad[3];
555: TensorVector(invJ, &grad[3 * bidx], lgrad);
556: for (ll = 0; ll < 3; ll++) {
557: PetscInt teidx = ll + 3 * (ii + jj * NB + kk * NB * NB);
558: ej[teidx + NPB * tridx] += scl * (lgrad[0] * dFS[3 * ll + 0] + lgrad[1] * dFS[3 * ll + 1] + lgrad[2] * dFS[3 * ll + 2]);
559: }
560: }
561: }
562: } /* end of testfunctions */
563: }
564: }
565: }
566: } /* end of trialfunctions */
567: }
568: }
569: }
570: } /* end of quadrature points */
571: }
573: void ApplyBCsElement(PetscInt mx, PetscInt my, PetscInt mz, PetscInt i, PetscInt j, PetscInt k, PetscScalar *jacobian)
574: {
575: PetscInt ii, jj, kk, ll, ei, ej, ek, el;
576: for (kk = 0; kk < NB; kk++) {
577: for (jj = 0; jj < NB; jj++) {
578: for (ii = 0; ii < NB; ii++) {
579: for (ll = 0; ll < 3; ll++) {
580: PetscInt tridx = ll + 3 * (ii + jj * NB + kk * NB * NB);
581: for (ek = 0; ek < NB; ek++) {
582: for (ej = 0; ej < NB; ej++) {
583: for (ei = 0; ei < NB; ei++) {
584: for (el = 0; el < 3; el++) {
585: if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz) || OnBoundary(i + ei, j + ej, k + ek, mx, my, mz)) {
586: PetscInt teidx = el + 3 * (ei + ej * NB + ek * NB * NB);
587: if (teidx == tridx) {
588: jacobian[tridx + NPB * teidx] = 1.;
589: } else {
590: jacobian[tridx + NPB * teidx] = 0.;
591: }
592: }
593: }
594: }
595: }
596: }
597: }
598: }
599: }
600: }
601: }
603: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field ***x, Mat jacpre, Mat jac, void *ptr)
604: {
605: /* values for each basis function at each quadrature point */
606: AppCtx *user = (AppCtx *)ptr;
607: PetscInt i, j, k, m, l;
608: PetscInt ii, jj, kk;
609: PetscScalar ej[NPB * NPB];
610: PetscScalar vals[NPB * NPB];
611: Field ex[NEB];
612: CoordField ec[NEB];
614: PetscInt xs = info->xs, ys = info->ys, zs = info->zs;
615: PetscInt xm = info->xm, ym = info->ym, zm = info->zm;
616: PetscInt xes, yes, zes, xee, yee, zee;
617: PetscInt mx = info->mx, my = info->my, mz = info->mz;
618: DM cda;
619: CoordField ***c;
620: Vec C;
621: PetscInt nrows;
622: MatStencil col[NPB], row[NPB];
623: PetscScalar v[9];
625: PetscFunctionBeginUser;
626: PetscCall(DMGetCoordinateDM(info->da, &cda));
627: PetscCall(DMGetCoordinatesLocal(info->da, &C));
628: PetscCall(DMDAVecGetArray(cda, C, &c));
629: PetscCall(MatScale(jac, 0.0));
631: xes = xs;
632: yes = ys;
633: zes = zs;
634: xee = xs + xm;
635: yee = ys + ym;
636: zee = zs + zm;
637: if (xs > 0) xes = xs - 1;
638: if (ys > 0) yes = ys - 1;
639: if (zs > 0) zes = zs - 1;
640: if (xs + xm == mx) xee = xs + xm - 1;
641: if (ys + ym == my) yee = ys + ym - 1;
642: if (zs + zm == mz) zee = zs + zm - 1;
643: for (k = zes; k < zee; k++) {
644: for (j = yes; j < yee; j++) {
645: for (i = xes; i < xee; i++) {
646: GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
647: FormElementJacobian(ex, ec, NULL, NULL, ej, user);
648: ApplyBCsElement(mx, my, mz, i, j, k, ej);
649: nrows = 0.;
650: for (kk = 0; kk < NB; kk++) {
651: for (jj = 0; jj < NB; jj++) {
652: for (ii = 0; ii < NB; ii++) {
653: PetscInt idx = ii + jj * 2 + kk * 4;
654: for (m = 0; m < 3; m++) {
655: col[3 * idx + m].i = i + ii;
656: col[3 * idx + m].j = j + jj;
657: col[3 * idx + m].k = k + kk;
658: col[3 * idx + m].c = m;
659: if (i + ii >= xs && i + ii < xm + xs && j + jj >= ys && j + jj < ys + ym && k + kk >= zs && k + kk < zs + zm) {
660: row[nrows].i = i + ii;
661: row[nrows].j = j + jj;
662: row[nrows].k = k + kk;
663: row[nrows].c = m;
664: for (l = 0; l < NPB; l++) vals[NPB * nrows + l] = ej[NPB * (3 * idx + m) + l];
665: nrows++;
666: }
667: }
668: }
669: }
670: }
671: PetscCall(MatSetValuesStencil(jac, nrows, row, NPB, col, vals, ADD_VALUES));
672: }
673: }
674: }
676: PetscCall(MatAssemblyBegin(jac, MAT_FLUSH_ASSEMBLY));
677: PetscCall(MatAssemblyEnd(jac, MAT_FLUSH_ASSEMBLY));
679: /* set the diagonal */
680: v[0] = 1.;
681: v[1] = 0.;
682: v[2] = 0.;
683: v[3] = 0.;
684: v[4] = 1.;
685: v[5] = 0.;
686: v[6] = 0.;
687: v[7] = 0.;
688: v[8] = 1.;
689: for (k = zs; k < zs + zm; k++) {
690: for (j = ys; j < ys + ym; j++) {
691: for (i = xs; i < xs + xm; i++) {
692: if (OnBoundary(i, j, k, mx, my, mz)) {
693: for (m = 0; m < 3; m++) {
694: col[m].i = i;
695: col[m].j = j;
696: col[m].k = k;
697: col[m].c = m;
698: }
699: PetscCall(MatSetValuesStencil(jac, 3, col, 3, col, v, INSERT_VALUES));
700: }
701: }
702: }
703: }
705: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
706: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
708: PetscCall(DMDAVecRestoreArray(cda, C, &c));
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
712: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field ***x, Field ***f, void *ptr)
713: {
714: /* values for each basis function at each quadrature point */
715: AppCtx *user = (AppCtx *)ptr;
716: PetscInt i, j, k, l;
717: PetscInt ii, jj, kk;
719: Field ef[NEB];
720: Field ex[NEB];
721: CoordField ec[NEB];
723: PetscInt xs = info->xs, ys = info->ys, zs = info->zs;
724: PetscInt xm = info->xm, ym = info->ym, zm = info->zm;
725: PetscInt xes, yes, zes, xee, yee, zee;
726: PetscInt mx = info->mx, my = info->my, mz = info->mz;
727: DM cda;
728: CoordField ***c;
729: Vec C;
731: PetscFunctionBeginUser;
732: PetscCall(DMGetCoordinateDM(info->da, &cda));
733: PetscCall(DMGetCoordinatesLocal(info->da, &C));
734: PetscCall(DMDAVecGetArray(cda, C, &c));
735: PetscCall(DMDAGetInfo(info->da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
736: PetscCall(DMDAGetCorners(info->da, &xs, &ys, &zs, &xm, &ym, &zm));
738: /* loop over elements */
739: for (k = zs; k < zs + zm; k++) {
740: for (j = ys; j < ys + ym; j++) {
741: for (i = xs; i < xs + xm; i++) {
742: for (l = 0; l < 3; l++) f[k][j][i][l] = 0.;
743: }
744: }
745: }
746: /* element starts and ends */
747: xes = xs;
748: yes = ys;
749: zes = zs;
750: xee = xs + xm;
751: yee = ys + ym;
752: zee = zs + zm;
753: if (xs > 0) xes = xs - 1;
754: if (ys > 0) yes = ys - 1;
755: if (zs > 0) zes = zs - 1;
756: if (xs + xm == mx) xee = xs + xm - 1;
757: if (ys + ym == my) yee = ys + ym - 1;
758: if (zs + zm == mz) zee = zs + zm - 1;
759: for (k = zes; k < zee; k++) {
760: for (j = yes; j < yee; j++) {
761: for (i = xes; i < xee; i++) {
762: GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
763: FormElementJacobian(ex, ec, ef, NULL, NULL, user);
764: /* put this element's additions into the residuals */
765: for (kk = 0; kk < NB; kk++) {
766: for (jj = 0; jj < NB; jj++) {
767: for (ii = 0; ii < NB; ii++) {
768: PetscInt idx = ii + jj * NB + kk * NB * NB;
769: if (k + kk >= zs && j + jj >= ys && i + ii >= xs && k + kk < zs + zm && j + jj < ys + ym && i + ii < xs + xm) {
770: if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
771: 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];
772: } else {
773: for (l = 0; l < 3; l++) f[k + kk][j + jj][i + ii][l] += ef[idx][l];
774: }
775: }
776: }
777: }
778: }
779: }
780: }
781: }
782: PetscCall(DMDAVecRestoreArray(cda, C, &c));
783: PetscFunctionReturn(PETSC_SUCCESS);
784: }
786: PetscErrorCode TangentLoad(SNES snes, Vec X, Vec Q, void *ptr)
787: {
788: /* values for each basis function at each quadrature point */
789: AppCtx *user = (AppCtx *)ptr;
790: PetscInt xs, ys, zs;
791: PetscInt xm, ym, zm;
792: PetscInt mx, my, mz;
793: DM da;
794: Vec Xl, Ql;
795: Field ***x, ***q;
796: PetscInt i, j, k, l;
797: PetscInt ii, jj, kk;
799: Field eq[NEB];
800: Field ex[NEB];
801: CoordField ec[NEB];
803: PetscInt xes, yes, zes, xee, yee, zee;
804: DM cda;
805: CoordField ***c;
806: Vec C;
808: PetscFunctionBeginUser;
809: /* update user context with current load parameter */
810: PetscCall(SNESNewtonALGetLoadParameter(snes, &user->load_factor));
812: PetscCall(SNESGetDM(snes, &da));
813: PetscCall(DMGetLocalVector(da, &Xl));
814: PetscCall(DMGetLocalVector(da, &Ql));
815: PetscCall(DMGlobalToLocal(da, X, INSERT_VALUES, Xl));
817: PetscCall(DMDAVecGetArray(da, Xl, &x));
818: PetscCall(DMDAVecGetArray(da, Ql, &q));
820: PetscCall(DMGetCoordinateDM(da, &cda));
821: PetscCall(DMGetCoordinatesLocal(da, &C));
822: PetscCall(DMDAVecGetArray(cda, C, &c));
823: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
824: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
826: /* loop over elements */
827: for (k = zs; k < zs + zm; k++) {
828: for (j = ys; j < ys + ym; j++) {
829: for (i = xs; i < xs + xm; i++) {
830: for (l = 0; l < 3; l++) q[k][j][i][l] = 0.;
831: }
832: }
833: }
834: /* element starts and ends */
835: xes = xs;
836: yes = ys;
837: zes = zs;
838: xee = xs + xm;
839: yee = ys + ym;
840: zee = zs + zm;
841: if (xs > 0) xes = xs - 1;
842: if (ys > 0) yes = ys - 1;
843: if (zs > 0) zes = zs - 1;
844: if (xs + xm == mx) xee = xs + xm - 1;
845: if (ys + ym == my) yee = ys + ym - 1;
846: if (zs + zm == mz) zee = zs + zm - 1;
847: for (k = zes; k < zee; k++) {
848: for (j = yes; j < yee; j++) {
849: for (i = xes; i < xee; i++) {
850: GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
851: FormElementJacobian(ex, ec, NULL, eq, NULL, user);
852: /* put this element's additions into the residuals */
853: for (kk = 0; kk < NB; kk++) {
854: for (jj = 0; jj < NB; jj++) {
855: for (ii = 0; ii < NB; ii++) {
856: PetscInt idx = ii + jj * NB + kk * NB * NB;
857: if (k + kk >= zs && j + jj >= ys && i + ii >= xs && k + kk < zs + zm && j + jj < ys + ym && i + ii < xs + xm) {
858: if (!OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
859: for (l = 0; l < 3; l++) q[k + kk][j + jj][i + ii][l] += eq[idx][l];
860: }
861: }
862: }
863: }
864: }
865: }
866: }
867: }
868: PetscCall(DMDAVecRestoreArray(da, Xl, &x));
869: PetscCall(DMDAVecRestoreArray(da, Ql, &q));
870: PetscCall(VecZeroEntries(Q));
871: PetscCall(DMLocalToGlobal(da, Ql, INSERT_VALUES, Q));
872: PetscCall(DMRestoreLocalVector(da, &Ql));
873: PetscCall(DMRestoreLocalVector(da, &Xl));
874: PetscCall(DMDAVecRestoreArray(cda, C, &c));
875: PetscFunctionReturn(PETSC_SUCCESS);
876: }
878: PetscErrorCode FormCoordinates(DM da, AppCtx *user)
879: {
880: Vec coords;
881: DM cda;
882: PetscInt mx, my, mz;
883: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
884: CoordField ***x;
886: PetscFunctionBeginUser;
887: PetscCall(DMGetCoordinateDM(da, &cda));
888: PetscCall(DMCreateGlobalVector(cda, &coords));
889: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
890: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
891: PetscCall(DMDAVecGetArray(da, coords, &x));
892: for (k = zs; k < zs + zm; k++) {
893: for (j = ys; j < ys + ym; j++) {
894: for (i = xs; i < xs + xm; i++) {
895: PetscReal cx = ((PetscReal)i) / ((PetscReal)(mx - 1));
896: PetscReal cy = ((PetscReal)j) / ((PetscReal)(my - 1));
897: PetscReal cz = ((PetscReal)k) / ((PetscReal)(mz - 1));
898: PetscReal rad = user->rad + cy * user->height;
899: PetscReal ang = (cx - 0.5) * user->arc;
900: x[k][j][i][0] = rad * PetscSinReal(ang);
901: x[k][j][i][1] = rad * PetscCosReal(ang) - (user->rad + 0.5 * user->height) * PetscCosReal(-0.5 * user->arc);
902: x[k][j][i][2] = user->width * (cz - 0.5);
903: }
904: }
905: }
906: PetscCall(DMDAVecRestoreArray(da, coords, &x));
907: PetscCall(DMSetCoordinates(da, coords));
908: PetscCall(VecDestroy(&coords));
909: PetscFunctionReturn(PETSC_SUCCESS);
910: }
912: PetscErrorCode InitialGuess(DM da, AppCtx *user, Vec X)
913: {
914: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
915: PetscInt mx, my, mz;
916: Field ***x;
918: PetscFunctionBeginUser;
919: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
920: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
921: PetscCall(DMDAVecGetArray(da, X, &x));
923: for (k = zs; k < zs + zm; k++) {
924: for (j = ys; j < ys + ym; j++) {
925: for (i = xs; i < xs + xm; i++) {
926: /* reference coordinates */
927: PetscReal p_x = ((PetscReal)i) / ((PetscReal)(mx - 1));
928: PetscReal p_y = ((PetscReal)j) / ((PetscReal)(my - 1));
929: PetscReal p_z = ((PetscReal)k) / ((PetscReal)(mz - 1));
930: PetscReal o_x = p_x;
931: PetscReal o_y = p_y;
932: PetscReal o_z = p_z;
933: x[k][j][i][0] = o_x - p_x;
934: x[k][j][i][1] = o_y - p_y;
935: x[k][j][i][2] = o_z - p_z;
936: }
937: }
938: }
939: PetscCall(DMDAVecRestoreArray(da, X, &x));
940: PetscFunctionReturn(PETSC_SUCCESS);
941: }
943: PetscErrorCode ArcLengthScaling(DM da, AppCtx *user, Vec V)
944: {
945: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
946: PetscInt mx, my, mz;
947: Field ***v;
948: PetscReal rad = user->rad + user->height;
950: PetscFunctionBeginUser;
951: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
952: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
953: PetscCall(DMDAVecGetArray(da, V, &v));
955: for (k = zs; k < zs + zm; k++) {
956: for (j = ys; j < ys + ym; j++) {
957: for (i = xs; i < xs + xm; i++) {
958: v[k][j][i][0] = 2. / rad;
959: v[k][j][i][1] = 2. / rad;
960: v[k][j][i][2] = 2. / user->width;
961: }
962: }
963: }
964: PetscCall(DMDAVecRestoreArray(da, V, &v));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: PetscErrorCode FormRHS(DM da, AppCtx *user, Vec X)
969: {
970: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
971: PetscInt mx, my, mz;
972: Field ***x;
974: PetscFunctionBeginUser;
975: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
976: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
977: PetscCall(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: PetscCall(DMDAVecRestoreArray(da, X, &x));
990: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
1003: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1004: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1005: PetscCall(SNESGetDM(snes, &da));
1006: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1007: PetscCall(DMGetCoordinateDM(da, &cda));
1008: PetscCall(DMGetCoordinates(da, &C));
1009: j = my / 2;
1010: k = mz / 2;
1011: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
1012: PetscCall(DMDAVecGetArray(da, X, &x));
1013: PetscCall(DMDAVecGetArray(cda, C, &c));
1014: for (r = 0; r < size; r++) {
1015: if (rank == r) {
1016: if (j >= ys && j < ys + ym && k >= zs && k < zs + zm) {
1017: for (i = xs; i < xs + xm; i++) {
1018: 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])));
1019: }
1020: }
1021: }
1022: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
1023: }
1024: PetscCall(DMDAVecRestoreArray(da, X, &x));
1025: PetscCall(DMDAVecRestoreArray(cda, C, &c));
1026: PetscFunctionReturn(PETSC_SUCCESS);
1027: }
1029: /*TEST
1031: test:
1032: nsize: 2
1033: 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
1034: requires: !single
1035: timeoutfactor: 3
1037: test:
1038: suffix: 2
1039: 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 -npc_fas_levels_snes_linesearch_maxlambda 2.0
1040: requires: !single
1042: test:
1043: suffix: 3
1044: 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
1045: requires: !single
1047: test:
1048: suffix: 4
1049: 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
1050: requires: !single defined(PETSC_USE_INFO)
1051: filter: grep -h -e "KSP Residual norm" -e "SNES Function norm" -e "Number of SNES iterations" -e "mu:"
1053: test:
1054: suffix: 5
1055: 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
1056: requires: !single
1058: test:
1059: suffix: 6
1060: 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
1061: requires: !single
1063: test:
1064: nsize: {{1 2}}
1065: suffix: 7
1066: args: -dm_vec_type kokkos -dm_mat_type aijkokkos -da_refine 1 -rad 10.0 -young 10. -ploading -1. -loading -1. -rescale -pc_type mg -mg_levels_ksp_max_it 2 -snes_monitor_short -snes_type newtonal -snes_newtonal_step_size 10 -snes_newtonal_correction_type exact -ksp_rtol 1e-4
1067: requires: kokkos_kernels !complex !single
1068: timeoutfactor: 3
1070: test:
1071: nsize: {{1 2}}
1072: suffix: 8
1073: args: -dm_vec_type cuda -dm_mat_type aijcusparse -da_refine 1 -rad 10.0 -young 10. -ploading -1. -loading -1. -rescale -pc_type mg -mg_levels_ksp_max_it 2 -snes_monitor_short -snes_type newtonal -snes_newtonal_step_size 10 -snes_newtonal_correction_type exact -ksp_rtol 1e-4
1074: requires: cuda !complex !single
1075: timeoutfactor: 3
1077: test:
1078: nsize: {{1 2}}
1079: suffix: 9
1080: args: -da_refine 2 -rad 10.0 -young 10. -ploading -1. -loading -1. -rescale -pc_type mg -mg_levels_ksp_max_it 2 -snes_monitor_short -snes_type newtonal -snes_newtonal_step_size 10 -snes_newtonal_correction_type normal -ksp_rtol 1e-4
1081: requires: !single
1083: TEST*/