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