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.
38: The initial setup follows the example in pg. 268 of "Nonlinear Finite Element Methods" by Peter Wriggers, but is a
39: 3D extension.
41: ------------------------------------------------------------------------F*/
43: #include <petscsnes.h>
44: #include <petscdm.h>
45: #include <petscdmda.h>
47: #define QP0 0.2113248654051871
48: #define QP1 0.7886751345948129
49: #define NQ 2
50: #define NB 2
51: #define NEB 8
52: #define NEQ 8
53: #define NPB 24
55: #define NVALS NEB *NEQ
56: const PetscReal pts[NQ] = {QP0, QP1};
57: const PetscReal wts[NQ] = {0.5, 0.5};
59: PetscScalar vals[NVALS];
60: PetscScalar grad[3 * NVALS];
62: typedef PetscScalar Field[3];
63: typedef PetscScalar CoordField[3];
65: typedef PetscScalar JacField[9];
67: PetscErrorCode FormFunctionLocal(DMDALocalInfo *, Field ***, Field ***, void *);
68: PetscErrorCode FormJacobianLocal(DMDALocalInfo *, Field ***, Mat, Mat, void *);
69: PetscErrorCode DisplayLine(SNES, Vec);
70: PetscErrorCode FormElements(void);
72: typedef struct {
73: PetscReal loading;
74: PetscReal mu;
75: PetscReal lambda;
76: PetscReal rad;
77: PetscReal height;
78: PetscReal width;
79: PetscReal arc;
80: PetscReal ploading;
81: } AppCtx;
83: PetscErrorCode InitialGuess(DM, AppCtx *, Vec);
84: PetscErrorCode FormRHS(DM, AppCtx *, Vec);
85: PetscErrorCode FormCoordinates(DM, AppCtx *);
86: extern PetscErrorCode NonlinearGS(SNES, Vec, Vec, void *);
88: int main(int argc, char **argv)
89: {
90: AppCtx user; /* user-defined work context */
91: PetscInt mx, my, its;
92: MPI_Comm comm;
93: SNES snes;
94: DM da;
95: Vec x, X, b;
96: PetscBool youngflg, poissonflg, muflg, lambdaflg, view = PETSC_FALSE, viewline = PETSC_FALSE;
97: PetscReal poisson = 0.2, young = 4e4;
98: char filename[PETSC_MAX_PATH_LEN] = "ex16.vts";
99: char filename_def[PETSC_MAX_PATH_LEN] = "ex16_def.vts";
101: PetscFunctionBeginUser;
102: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
103: PetscCall(FormElements());
104: comm = PETSC_COMM_WORLD;
105: PetscCall(SNESCreate(comm, &snes));
106: 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));
107: PetscCall(DMSetFromOptions(da));
108: PetscCall(DMSetUp(da));
109: PetscCall(SNESSetDM(snes, (DM)da));
111: PetscCall(SNESSetNGS(snes, NonlinearGS, &user));
113: 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));
114: user.loading = 0.0;
115: user.arc = PETSC_PI / 3.;
116: user.mu = 4.0;
117: user.lambda = 1.0;
118: user.rad = 100.0;
119: user.height = 3.;
120: user.width = 1.;
121: user.ploading = -5e3;
123: PetscCall(PetscOptionsGetReal(NULL, NULL, "-arc", &user.arc, NULL));
124: PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, &muflg));
125: PetscCall(PetscOptionsGetReal(NULL, NULL, "-lambda", &user.lambda, &lambdaflg));
126: PetscCall(PetscOptionsGetReal(NULL, NULL, "-rad", &user.rad, NULL));
127: PetscCall(PetscOptionsGetReal(NULL, NULL, "-height", &user.height, NULL));
128: PetscCall(PetscOptionsGetReal(NULL, NULL, "-width", &user.width, NULL));
129: PetscCall(PetscOptionsGetReal(NULL, NULL, "-loading", &user.loading, NULL));
130: PetscCall(PetscOptionsGetReal(NULL, NULL, "-ploading", &user.ploading, NULL));
131: PetscCall(PetscOptionsGetReal(NULL, NULL, "-poisson", &poisson, &poissonflg));
132: PetscCall(PetscOptionsGetReal(NULL, NULL, "-young", &young, &youngflg));
133: if ((youngflg || poissonflg) || !(muflg || lambdaflg)) {
134: /* set the lame' parameters based upon the poisson ratio and young's modulus */
135: user.lambda = poisson * young / ((1. + poisson) * (1. - 2. * poisson));
136: user.mu = young / (2. * (1. + poisson));
137: }
138: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view", &view, NULL));
139: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_line", &viewline, NULL));
141: PetscCall(DMDASetFieldName(da, 0, "x_disp"));
142: PetscCall(DMDASetFieldName(da, 1, "y_disp"));
143: PetscCall(DMDASetFieldName(da, 2, "z_disp"));
145: PetscCall(DMSetApplicationContext(da, &user));
146: PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode(*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user));
147: PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobian)FormJacobianLocal, &user));
148: PetscCall(SNESSetFromOptions(snes));
149: PetscCall(FormCoordinates(da, &user));
151: PetscCall(DMCreateGlobalVector(da, &x));
152: PetscCall(DMCreateGlobalVector(da, &b));
153: PetscCall(InitialGuess(da, &user, x));
154: PetscCall(FormRHS(da, &user, b));
156: PetscCall(PetscPrintf(comm, "lambda: %f mu: %f\n", (double)user.lambda, (double)user.mu));
158: /* show a cross-section of the initial state */
159: if (viewline) PetscCall(DisplayLine(snes, x));
161: /* get the loaded configuration */
162: PetscCall(SNESSolve(snes, b, x));
164: PetscCall(SNESGetIterationNumber(snes, &its));
165: PetscCall(PetscPrintf(comm, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
166: PetscCall(SNESGetSolution(snes, &X));
167: /* show a cross-section of the final state */
168: if (viewline) PetscCall(DisplayLine(snes, X));
170: if (view) {
171: PetscViewer viewer;
172: Vec coords;
173: PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename, FILE_MODE_WRITE, &viewer));
174: PetscCall(VecView(x, viewer));
175: PetscCall(PetscViewerDestroy(&viewer));
176: PetscCall(DMGetCoordinates(da, &coords));
177: PetscCall(VecAXPY(coords, 1.0, x));
178: PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename_def, FILE_MODE_WRITE, &viewer));
179: PetscCall(VecView(x, viewer));
180: PetscCall(PetscViewerDestroy(&viewer));
181: }
183: PetscCall(VecDestroy(&x));
184: PetscCall(VecDestroy(&b));
185: PetscCall(DMDestroy(&da));
186: PetscCall(SNESDestroy(&snes));
187: PetscCall(PetscFinalize());
188: return 0;
189: }
191: PetscInt OnBoundary(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz)
192: {
193: if ((i == 0 || i == mx - 1) && j == my - 1) return 1;
194: return 0;
195: }
197: void BoundaryValue(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz, PetscScalar *val, AppCtx *user)
198: {
199: /* reference coordinates */
200: PetscReal p_x = ((PetscReal)i) / (((PetscReal)(mx - 1)));
201: PetscReal p_y = ((PetscReal)j) / (((PetscReal)(my - 1)));
202: PetscReal p_z = ((PetscReal)k) / (((PetscReal)(mz - 1)));
203: PetscReal o_x = p_x;
204: PetscReal o_y = p_y;
205: PetscReal o_z = p_z;
206: val[0] = o_x - p_x;
207: val[1] = o_y - p_y;
208: val[2] = o_z - p_z;
209: }
211: void InvertTensor(PetscScalar *t, PetscScalar *ti, PetscReal *dett)
212: {
213: const PetscScalar a = t[0];
214: const PetscScalar b = t[1];
215: const PetscScalar c = t[2];
216: const PetscScalar d = t[3];
217: const PetscScalar e = t[4];
218: const PetscScalar f = t[5];
219: const PetscScalar g = t[6];
220: const PetscScalar h = t[7];
221: const PetscScalar i = t[8];
222: const PetscReal det = PetscRealPart(a * (e * i - f * h) - b * (i * d - f * g) + c * (d * h - e * g));
223: const PetscReal di = 1. / det;
224: if (ti) {
225: const PetscScalar A = (e * i - f * h);
226: const PetscScalar B = -(d * i - f * g);
227: const PetscScalar C = (d * h - e * g);
228: const PetscScalar D = -(b * i - c * h);
229: const PetscScalar E = (a * i - c * g);
230: const PetscScalar F = -(a * h - b * g);
231: const PetscScalar G = (b * f - c * e);
232: const PetscScalar H = -(a * f - c * d);
233: const PetscScalar II = (a * e - b * d);
234: ti[0] = di * A;
235: ti[1] = di * D;
236: ti[2] = di * G;
237: ti[3] = di * B;
238: ti[4] = di * E;
239: ti[5] = di * H;
240: ti[6] = di * C;
241: ti[7] = di * F;
242: ti[8] = di * II;
243: }
244: if (dett) *dett = det;
245: }
247: void TensorTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
248: {
249: PetscInt i, j, m;
250: for (i = 0; i < 3; i++) {
251: for (j = 0; j < 3; j++) {
252: c[i + 3 * j] = 0;
253: for (m = 0; m < 3; m++) c[i + 3 * j] += a[m + 3 * j] * b[i + 3 * m];
254: }
255: }
256: }
258: void TensorTransposeTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
259: {
260: PetscInt i, j, m;
261: for (i = 0; i < 3; i++) {
262: for (j = 0; j < 3; j++) {
263: c[i + 3 * j] = 0;
264: for (m = 0; m < 3; m++) c[i + 3 * j] += a[3 * m + j] * b[i + 3 * m];
265: }
266: }
267: }
269: void TensorVector(PetscScalar *rot, PetscScalar *vec, PetscScalar *tvec)
270: {
271: tvec[0] = rot[0] * vec[0] + rot[1] * vec[1] + rot[2] * vec[2];
272: tvec[1] = rot[3] * vec[0] + rot[4] * vec[1] + rot[5] * vec[2];
273: tvec[2] = rot[6] * vec[0] + rot[7] * vec[1] + rot[8] * vec[2];
274: }
276: void DeformationGradient(Field *ex, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *invJ, PetscScalar *F)
277: {
278: PetscInt ii, jj, kk, l;
279: for (l = 0; l < 9; l++) F[l] = 0.;
280: F[0] = 1.;
281: F[4] = 1.;
282: F[8] = 1.;
283: /* form the deformation gradient at this basis function -- loop over element unknowns */
284: for (kk = 0; kk < NB; kk++) {
285: for (jj = 0; jj < NB; jj++) {
286: for (ii = 0; ii < NB; ii++) {
287: PetscInt idx = ii + jj * NB + kk * NB * NB;
288: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
289: PetscScalar lgrad[3];
290: TensorVector(invJ, &grad[3 * bidx], lgrad);
291: F[0] += lgrad[0] * ex[idx][0];
292: F[1] += lgrad[1] * ex[idx][0];
293: F[2] += lgrad[2] * ex[idx][0];
294: F[3] += lgrad[0] * ex[idx][1];
295: F[4] += lgrad[1] * ex[idx][1];
296: F[5] += lgrad[2] * ex[idx][1];
297: F[6] += lgrad[0] * ex[idx][2];
298: F[7] += lgrad[1] * ex[idx][2];
299: F[8] += lgrad[2] * ex[idx][2];
300: }
301: }
302: }
303: }
305: void DeformationGradientJacobian(PetscInt qi, PetscInt qj, PetscInt qk, PetscInt ii, PetscInt jj, PetscInt kk, PetscInt fld, PetscScalar *invJ, PetscScalar *dF)
306: {
307: PetscInt l;
308: PetscScalar lgrad[3];
309: PetscInt idx = ii + jj * NB + kk * NB * NB;
310: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
311: for (l = 0; l < 9; l++) dF[l] = 0.;
312: /* form the deformation gradient at this basis function -- loop over element unknowns */
313: TensorVector(invJ, &grad[3 * bidx], lgrad);
314: dF[3 * fld] = lgrad[0];
315: dF[3 * fld + 1] = lgrad[1];
316: dF[3 * fld + 2] = lgrad[2];
317: }
319: void LagrangeGreenStrain(PetscScalar *F, PetscScalar *E)
320: {
321: PetscInt i, j, m;
322: for (i = 0; i < 3; i++) {
323: for (j = 0; j < 3; j++) {
324: E[i + 3 * j] = 0;
325: for (m = 0; m < 3; m++) E[i + 3 * j] += 0.5 * F[3 * m + j] * F[i + 3 * m];
326: }
327: }
328: for (i = 0; i < 3; i++) E[i + 3 * i] -= 0.5;
329: }
331: void SaintVenantKirchoff(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *S)
332: {
333: PetscInt i, j;
334: PetscScalar E[9];
335: PetscScalar trE = 0;
336: LagrangeGreenStrain(F, E);
337: for (i = 0; i < 3; i++) trE += E[i + 3 * i];
338: for (i = 0; i < 3; i++) {
339: for (j = 0; j < 3; j++) {
340: S[i + 3 * j] = 2. * mu * E[i + 3 * j];
341: if (i == j) S[i + 3 * j] += trE * lambda;
342: }
343: }
344: }
346: void SaintVenantKirchoffJacobian(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *dF, PetscScalar *dS)
347: {
348: PetscScalar FtdF[9], dE[9];
349: PetscInt i, j;
350: PetscScalar dtrE = 0.;
351: TensorTransposeTensor(dF, F, dE);
352: TensorTransposeTensor(F, dF, FtdF);
353: for (i = 0; i < 9; i++) dE[i] += FtdF[i];
354: for (i = 0; i < 9; i++) dE[i] *= 0.5;
355: for (i = 0; i < 3; i++) dtrE += dE[i + 3 * i];
356: for (i = 0; i < 3; i++) {
357: for (j = 0; j < 3; j++) {
358: dS[i + 3 * j] = 2. * mu * dE[i + 3 * j];
359: if (i == j) dS[i + 3 * j] += dtrE * lambda;
360: }
361: }
362: }
364: PetscErrorCode FormElements()
365: {
366: PetscInt i, j, k, ii, jj, kk;
367: PetscReal bx, by, bz, dbx, dby, dbz;
369: PetscFunctionBegin;
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: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
641: PetscCall(DMGetCoordinateDM(info->da, &cda));
642: PetscCall(DMGetCoordinatesLocal(info->da, &C));
643: PetscCall(DMDAVecGetArray(cda, C, &c));
644: PetscCall(MatScale(jac, 0.0));
646: xes = xs;
647: yes = ys;
648: zes = zs;
649: xee = xs + xm;
650: yee = ys + ym;
651: zee = zs + zm;
652: if (xs > 0) xes = xs - 1;
653: if (ys > 0) yes = ys - 1;
654: if (zs > 0) zes = zs - 1;
655: if (xs + xm == mx) xee = xs + xm - 1;
656: if (ys + ym == my) yee = ys + ym - 1;
657: if (zs + zm == mz) zee = zs + zm - 1;
658: for (k = zes; k < zee; k++) {
659: for (j = yes; j < yee; j++) {
660: for (i = xes; i < xee; i++) {
661: GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
662: FormElementJacobian(ex, ec, NULL, ej, user);
663: ApplyBCsElement(mx, my, mz, i, j, k, ej);
664: nrows = 0.;
665: for (kk = 0; kk < NB; kk++) {
666: for (jj = 0; jj < NB; jj++) {
667: for (ii = 0; ii < NB; ii++) {
668: PetscInt idx = ii + jj * 2 + kk * 4;
669: for (m = 0; m < 3; m++) {
670: col[3 * idx + m].i = i + ii;
671: col[3 * idx + m].j = j + jj;
672: col[3 * idx + m].k = k + kk;
673: col[3 * idx + m].c = m;
674: if (i + ii >= xs && i + ii < xm + xs && j + jj >= ys && j + jj < ys + ym && k + kk >= zs && k + kk < zs + zm) {
675: row[nrows].i = i + ii;
676: row[nrows].j = j + jj;
677: row[nrows].k = k + kk;
678: row[nrows].c = m;
679: for (l = 0; l < NPB; l++) vals[NPB * nrows + l] = ej[NPB * (3 * idx + m) + l];
680: nrows++;
681: }
682: }
683: }
684: }
685: }
686: PetscCall(MatSetValuesStencil(jac, nrows, row, NPB, col, vals, ADD_VALUES));
687: }
688: }
689: }
691: PetscCall(MatAssemblyBegin(jac, MAT_FLUSH_ASSEMBLY));
692: PetscCall(MatAssemblyEnd(jac, MAT_FLUSH_ASSEMBLY));
694: /* set the diagonal */
695: v[0] = 1.;
696: v[1] = 0.;
697: v[2] = 0.;
698: v[3] = 0.;
699: v[4] = 1.;
700: v[5] = 0.;
701: v[6] = 0.;
702: v[7] = 0.;
703: v[8] = 1.;
704: for (k = zs; k < zs + zm; k++) {
705: for (j = ys; j < ys + ym; j++) {
706: for (i = xs; i < xs + xm; i++) {
707: if (OnBoundary(i, j, k, mx, my, mz)) {
708: for (m = 0; m < 3; m++) {
709: col[m].i = i;
710: col[m].j = j;
711: col[m].k = k;
712: col[m].c = m;
713: }
714: PetscCall(MatSetValuesStencil(jac, 3, col, 3, col, v, INSERT_VALUES));
715: }
716: }
717: }
718: }
720: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
721: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
723: PetscCall(DMDAVecRestoreArray(cda, C, &c));
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field ***x, Field ***f, void *ptr)
728: {
729: /* values for each basis function at each quadrature point */
730: AppCtx *user = (AppCtx *)ptr;
731: PetscInt i, j, k, l;
732: PetscInt ii, jj, kk;
734: Field ef[NEB];
735: Field ex[NEB];
736: CoordField ec[NEB];
738: PetscInt xs = info->xs, ys = info->ys, zs = info->zs;
739: PetscInt xm = info->xm, ym = info->ym, zm = info->zm;
740: PetscInt xes, yes, zes, xee, yee, zee;
741: PetscInt mx = info->mx, my = info->my, mz = info->mz;
742: DM cda;
743: CoordField ***c;
744: Vec C;
746: PetscFunctionBegin;
747: PetscCall(DMGetCoordinateDM(info->da, &cda));
748: PetscCall(DMGetCoordinatesLocal(info->da, &C));
749: PetscCall(DMDAVecGetArray(cda, C, &c));
750: PetscCall(DMDAGetInfo(info->da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
751: PetscCall(DMDAGetCorners(info->da, &xs, &ys, &zs, &xm, &ym, &zm));
753: /* loop over elements */
754: for (k = zs; k < zs + zm; k++) {
755: for (j = ys; j < ys + ym; j++) {
756: for (i = xs; i < xs + xm; i++) {
757: for (l = 0; l < 3; l++) f[k][j][i][l] = 0.;
758: }
759: }
760: }
761: /* element starts and ends */
762: xes = xs;
763: yes = ys;
764: zes = zs;
765: xee = xs + xm;
766: yee = ys + ym;
767: zee = zs + zm;
768: if (xs > 0) xes = xs - 1;
769: if (ys > 0) yes = ys - 1;
770: if (zs > 0) zes = zs - 1;
771: if (xs + xm == mx) xee = xs + xm - 1;
772: if (ys + ym == my) yee = ys + ym - 1;
773: if (zs + zm == mz) zee = zs + zm - 1;
774: for (k = zes; k < zee; k++) {
775: for (j = yes; j < yee; j++) {
776: for (i = xes; i < xee; i++) {
777: GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
778: FormElementJacobian(ex, ec, ef, NULL, user);
779: /* put this element's additions into the residuals */
780: for (kk = 0; kk < NB; kk++) {
781: for (jj = 0; jj < NB; jj++) {
782: for (ii = 0; ii < NB; ii++) {
783: PetscInt idx = ii + jj * NB + kk * NB * NB;
784: if (k + kk >= zs && j + jj >= ys && i + ii >= xs && k + kk < zs + zm && j + jj < ys + ym && i + ii < xs + xm) {
785: if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
786: 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];
787: } else {
788: for (l = 0; l < 3; l++) f[k + kk][j + jj][i + ii][l] += ef[idx][l];
789: }
790: }
791: }
792: }
793: }
794: }
795: }
796: }
797: PetscCall(DMDAVecRestoreArray(cda, C, &c));
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: PetscErrorCode NonlinearGS(SNES snes, Vec X, Vec B, void *ptr)
802: {
803: /* values for each basis function at each quadrature point */
804: AppCtx *user = (AppCtx *)ptr;
805: PetscInt i, j, k, l, m, n, s;
806: PetscInt pi, pj, pk;
807: Field ef[1];
808: Field ex[8];
809: PetscScalar ej[9];
810: CoordField ec[8];
811: PetscScalar pjac[9], pjinv[9];
812: PetscScalar pf[3], py[3];
813: PetscInt xs, ys, zs;
814: PetscInt xm, ym, zm;
815: PetscInt mx, my, mz;
816: DM cda;
817: CoordField ***c;
818: Vec C;
819: DM da;
820: Vec Xl, Bl;
821: Field ***x, ***b;
822: PetscInt sweeps, its;
823: PetscReal atol, rtol, stol;
824: PetscReal fnorm0 = 0.0, fnorm, ynorm, xnorm = 0.0;
826: PetscFunctionBegin;
827: PetscCall(SNESNGSGetSweeps(snes, &sweeps));
828: PetscCall(SNESNGSGetTolerances(snes, &atol, &rtol, &stol, &its));
830: PetscCall(SNESGetDM(snes, &da));
831: PetscCall(DMGetLocalVector(da, &Xl));
832: if (B) PetscCall(DMGetLocalVector(da, &Bl));
833: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xl));
834: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xl));
835: if (B) {
836: PetscCall(DMGlobalToLocalBegin(da, B, INSERT_VALUES, Bl));
837: PetscCall(DMGlobalToLocalEnd(da, B, INSERT_VALUES, Bl));
838: }
839: PetscCall(DMDAVecGetArray(da, Xl, &x));
840: if (B) PetscCall(DMDAVecGetArray(da, Bl, &b));
842: PetscCall(DMGetCoordinateDM(da, &cda));
843: PetscCall(DMGetCoordinatesLocal(da, &C));
844: PetscCall(DMDAVecGetArray(cda, C, &c));
845: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
846: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
848: for (s = 0; s < sweeps; s++) {
849: for (k = zs; k < zs + zm; k++) {
850: for (j = ys; j < ys + ym; j++) {
851: for (i = xs; i < xs + xm; i++) {
852: if (OnBoundary(i, j, k, mx, my, mz)) {
853: BoundaryValue(i, j, k, mx, my, mz, x[k][j][i], user);
854: } else {
855: for (n = 0; n < its; n++) {
856: for (m = 0; m < 9; m++) pjac[m] = 0.;
857: for (m = 0; m < 3; m++) pf[m] = 0.;
858: /* gather the elements for this point */
859: for (pk = -1; pk < 1; pk++) {
860: for (pj = -1; pj < 1; pj++) {
861: for (pi = -1; pi < 1; pi++) {
862: /* check that this element exists */
863: if (i + pi >= 0 && i + pi < mx - 1 && j + pj >= 0 && j + pj < my - 1 && k + pk >= 0 && k + pk < mz - 1) {
864: /* create the element function and jacobian */
865: GatherElementData(mx, my, mz, x, c, i + pi, j + pj, k + pk, ex, ec, user);
866: FormPBJacobian(-pi, -pj, -pk, ex, ec, ef, ej, user);
867: /* extract the point named by i,j,k from the whole element jacobian and function */
868: for (l = 0; l < 3; l++) {
869: pf[l] += ef[0][l];
870: for (m = 0; m < 3; m++) pjac[3 * m + l] += ej[3 * m + l];
871: }
872: }
873: }
874: }
875: }
876: /* invert */
877: InvertTensor(pjac, pjinv, NULL);
878: /* apply */
879: if (B)
880: for (m = 0; m < 3; m++) pf[m] -= b[k][j][i][m];
881: TensorVector(pjinv, pf, py);
882: xnorm = 0.;
883: for (m = 0; m < 3; m++) {
884: x[k][j][i][m] -= py[m];
885: xnorm += PetscRealPart(x[k][j][i][m] * x[k][j][i][m]);
886: }
887: fnorm = PetscRealPart(pf[0] * pf[0] + pf[1] * pf[1] + pf[2] * pf[2]);
888: if (n == 0) fnorm0 = fnorm;
889: ynorm = PetscRealPart(py[0] * py[0] + py[1] * py[1] + py[2] * py[2]);
890: if (fnorm < atol * atol || fnorm < rtol * rtol * fnorm0 || ynorm < stol * stol * xnorm) break;
891: }
892: }
893: }
894: }
895: }
896: }
897: PetscCall(DMDAVecRestoreArray(da, Xl, &x));
898: PetscCall(DMLocalToGlobalBegin(da, Xl, INSERT_VALUES, X));
899: PetscCall(DMLocalToGlobalEnd(da, Xl, INSERT_VALUES, X));
900: PetscCall(DMRestoreLocalVector(da, &Xl));
901: if (B) {
902: PetscCall(DMDAVecRestoreArray(da, Bl, &b));
903: PetscCall(DMRestoreLocalVector(da, &Bl));
904: }
905: PetscCall(DMDAVecRestoreArray(cda, C, &c));
906: PetscFunctionReturn(PETSC_SUCCESS);
907: }
909: PetscErrorCode FormCoordinates(DM da, AppCtx *user)
910: {
911: Vec coords;
912: DM cda;
913: PetscInt mx, my, mz;
914: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
915: CoordField ***x;
917: PetscFunctionBegin;
918: PetscCall(DMGetCoordinateDM(da, &cda));
919: PetscCall(DMCreateGlobalVector(cda, &coords));
920: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
921: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
922: PetscCall(DMDAVecGetArray(da, coords, &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: PetscReal cx = ((PetscReal)i) / (((PetscReal)(mx - 1)));
927: PetscReal cy = ((PetscReal)j) / (((PetscReal)(my - 1)));
928: PetscReal cz = ((PetscReal)k) / (((PetscReal)(mz - 1)));
929: PetscReal rad = user->rad + cy * user->height;
930: PetscReal ang = (cx - 0.5) * user->arc;
931: x[k][j][i][0] = rad * PetscSinReal(ang);
932: x[k][j][i][1] = rad * PetscCosReal(ang) - (user->rad + 0.5 * user->height) * PetscCosReal(-0.5 * user->arc);
933: x[k][j][i][2] = user->width * (cz - 0.5);
934: }
935: }
936: }
937: PetscCall(DMDAVecRestoreArray(da, coords, &x));
938: PetscCall(DMSetCoordinates(da, coords));
939: PetscCall(VecDestroy(&coords));
940: PetscFunctionReturn(PETSC_SUCCESS);
941: }
943: PetscErrorCode InitialGuess(DM da, AppCtx *user, Vec X)
944: {
945: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
946: PetscInt mx, my, mz;
947: Field ***x;
949: PetscFunctionBegin;
950: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
951: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
952: PetscCall(DMDAVecGetArray(da, X, &x));
954: for (k = zs; k < zs + zm; k++) {
955: for (j = ys; j < ys + ym; j++) {
956: for (i = xs; i < xs + xm; i++) {
957: /* reference coordinates */
958: PetscReal p_x = ((PetscReal)i) / (((PetscReal)(mx - 1)));
959: PetscReal p_y = ((PetscReal)j) / (((PetscReal)(my - 1)));
960: PetscReal p_z = ((PetscReal)k) / (((PetscReal)(mz - 1)));
961: PetscReal o_x = p_x;
962: PetscReal o_y = p_y;
963: PetscReal o_z = p_z;
964: x[k][j][i][0] = o_x - p_x;
965: x[k][j][i][1] = o_y - p_y;
966: x[k][j][i][2] = o_z - p_z;
967: }
968: }
969: }
970: PetscCall(DMDAVecRestoreArray(da, X, &x));
971: PetscFunctionReturn(PETSC_SUCCESS);
972: }
974: PetscErrorCode FormRHS(DM da, AppCtx *user, Vec X)
975: {
976: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
977: PetscInt mx, my, mz;
978: Field ***x;
980: PetscFunctionBegin;
981: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
982: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
983: PetscCall(DMDAVecGetArray(da, X, &x));
985: for (k = zs; k < zs + zm; k++) {
986: for (j = ys; j < ys + ym; j++) {
987: for (i = xs; i < xs + xm; i++) {
988: x[k][j][i][0] = 0.;
989: x[k][j][i][1] = 0.;
990: x[k][j][i][2] = 0.;
991: if (i == (mx - 1) / 2 && j == (my - 1)) x[k][j][i][1] = user->ploading / (mz - 1);
992: }
993: }
994: }
995: PetscCall(DMDAVecRestoreArray(da, X, &x));
996: PetscFunctionReturn(PETSC_SUCCESS);
997: }
999: PetscErrorCode DisplayLine(SNES snes, Vec X)
1000: {
1001: PetscInt r, i, j = 0, k = 0, xs, xm, ys, ym, zs, zm, mx, my, mz;
1002: Field ***x;
1003: CoordField ***c;
1004: DM da, cda;
1005: Vec C;
1006: PetscMPIInt size, rank;
1008: PetscFunctionBegin;
1009: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1010: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1011: PetscCall(SNESGetDM(snes, &da));
1012: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1013: PetscCall(DMGetCoordinateDM(da, &cda));
1014: PetscCall(DMGetCoordinates(da, &C));
1015: j = my / 2;
1016: k = mz / 2;
1017: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
1018: PetscCall(DMDAVecGetArray(da, X, &x));
1019: PetscCall(DMDAVecGetArray(cda, C, &c));
1020: for (r = 0; r < size; r++) {
1021: if (rank == r) {
1022: if (j >= ys && j < ys + ym && k >= zs && k < zs + zm) {
1023: for (i = xs; i < xs + xm; i++) {
1024: 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])));
1025: }
1026: }
1027: }
1028: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
1029: }
1030: PetscCall(DMDAVecRestoreArray(da, X, &x));
1031: PetscCall(DMDAVecRestoreArray(cda, C, &c));
1032: PetscFunctionReturn(PETSC_SUCCESS);
1033: }
1035: /*TEST
1037: test:
1038: nsize: 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 -snes_max_it 7
1040: requires: !single
1041: timeoutfactor: 3
1043: test:
1044: suffix: 2
1045: 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
1046: requires: !single
1048: test:
1049: suffix: 3
1050: 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
1051: requires: !single
1053: TEST*/