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