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