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