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