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