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