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 ierr = FormElements();CHKERRQ(ierr); 110 comm = PETSC_COMM_WORLD; 111 ierr = SNESCreate(comm,&snes);CHKERRQ(ierr); 112 ierr = 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);CHKERRQ(ierr); 113 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 114 ierr = DMSetUp(da);CHKERRQ(ierr); 115 ierr = SNESSetDM(snes,(DM)da);CHKERRQ(ierr); 116 117 ierr = SNESSetNGS(snes,NonlinearGS,&user);CHKERRQ(ierr); 118 119 ierr = 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);CHKERRQ(ierr); 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 ierr = PetscOptionsGetReal(NULL,NULL,"-arc",&user.arc,NULL);CHKERRQ(ierr); 130 ierr = PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,&muflg);CHKERRQ(ierr); 131 ierr = PetscOptionsGetReal(NULL,NULL,"-lambda",&user.lambda,&lambdaflg);CHKERRQ(ierr); 132 ierr = PetscOptionsGetReal(NULL,NULL,"-rad",&user.rad,NULL);CHKERRQ(ierr); 133 ierr = PetscOptionsGetReal(NULL,NULL,"-height",&user.height,NULL);CHKERRQ(ierr); 134 ierr = PetscOptionsGetReal(NULL,NULL,"-width",&user.width,NULL);CHKERRQ(ierr); 135 ierr = PetscOptionsGetReal(NULL,NULL,"-loading",&user.loading,NULL);CHKERRQ(ierr); 136 ierr = PetscOptionsGetReal(NULL,NULL,"-ploading",&user.ploading,NULL);CHKERRQ(ierr); 137 ierr = PetscOptionsGetReal(NULL,NULL,"-poisson",&poisson,&poissonflg);CHKERRQ(ierr); 138 ierr = PetscOptionsGetReal(NULL,NULL,"-young",&young,&youngflg);CHKERRQ(ierr); 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 ierr = PetscOptionsGetBool(NULL,NULL,"-view",&view,NULL);CHKERRQ(ierr); 145 ierr = PetscOptionsGetBool(NULL,NULL,"-view_line",&viewline,NULL);CHKERRQ(ierr); 146 147 ierr = DMDASetFieldName(da,0,"x_disp");CHKERRQ(ierr); 148 ierr = DMDASetFieldName(da,1,"y_disp");CHKERRQ(ierr); 149 ierr = DMDASetFieldName(da,2,"z_disp");CHKERRQ(ierr); 150 151 ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr); 152 ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);CHKERRQ(ierr); 153 ierr = DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);CHKERRQ(ierr); 154 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 155 ierr = FormCoordinates(da,&user);CHKERRQ(ierr); 156 157 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr); 158 ierr = DMCreateGlobalVector(da,&b);CHKERRQ(ierr); 159 ierr = InitialGuess(da,&user,x);CHKERRQ(ierr); 160 ierr = FormRHS(da,&user,b);CHKERRQ(ierr); 161 162 ierr = PetscPrintf(comm,"lambda: %f mu: %f\n",(double)user.lambda,(double)user.mu);CHKERRQ(ierr); 163 164 /* show a cross-section of the initial state */ 165 if (viewline) { 166 ierr = DisplayLine(snes,x);CHKERRQ(ierr); 167 } 168 169 /* get the loaded configuration */ 170 ierr = SNESSolve(snes,b,x);CHKERRQ(ierr); 171 172 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 173 ierr = PetscPrintf(comm,"Number of SNES iterations = %D\n", its);CHKERRQ(ierr); 174 ierr = SNESGetSolution(snes,&X);CHKERRQ(ierr); 175 /* show a cross-section of the final state */ 176 if (viewline) { 177 ierr = DisplayLine(snes,X);CHKERRQ(ierr); 178 } 179 180 if (view) { 181 PetscViewer viewer; 182 Vec coords; 183 ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 184 ierr = VecView(x,viewer);CHKERRQ(ierr); 185 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 186 ierr = DMGetCoordinates(da,&coords);CHKERRQ(ierr); 187 ierr = VecAXPY(coords,1.0,x);CHKERRQ(ierr); 188 ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,filename_def,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 189 ierr = VecView(x,viewer);CHKERRQ(ierr); 190 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 191 } 192 193 ierr = VecDestroy(&x);CHKERRQ(ierr); 194 ierr = VecDestroy(&b);CHKERRQ(ierr); 195 ierr = DMDestroy(&da);CHKERRQ(ierr); 196 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 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 PetscErrorCode ierr; 638 PetscInt xs=info->xs,ys=info->ys,zs=info->zs; 639 PetscInt xm=info->xm,ym=info->ym,zm=info->zm; 640 PetscInt xes,yes,zes,xee,yee,zee; 641 PetscInt mx=info->mx,my=info->my,mz=info->mz; 642 DM cda; 643 CoordField ***c; 644 Vec C; 645 PetscInt nrows; 646 MatStencil col[NPB],row[NPB]; 647 PetscScalar v[9]; 648 649 PetscFunctionBegin; 650 ierr = DMGetCoordinateDM(info->da,&cda);CHKERRQ(ierr); 651 ierr = DMGetCoordinatesLocal(info->da,&C);CHKERRQ(ierr); 652 ierr = DMDAVecGetArray(cda,C,&c);CHKERRQ(ierr); 653 ierr = MatScale(jac,0.0);CHKERRQ(ierr); 654 655 xes = xs; 656 yes = ys; 657 zes = zs; 658 xee = xs+xm; 659 yee = ys+ym; 660 zee = zs+zm; 661 if (xs > 0) xes = xs-1; 662 if (ys > 0) yes = ys-1; 663 if (zs > 0) zes = zs-1; 664 if (xs+xm == mx) xee = xs+xm-1; 665 if (ys+ym == my) yee = ys+ym-1; 666 if (zs+zm == mz) zee = zs+zm-1; 667 for (k=zes; k<zee; k++) { 668 for (j=yes; j<yee; j++) { 669 for (i=xes; i<xee; i++) { 670 GatherElementData(mx,my,mz,x,c,i,j,k,ex,ec,user); 671 FormElementJacobian(ex,ec,NULL,ej,user); 672 ApplyBCsElement(mx,my,mz,i,j,k,ej); 673 nrows = 0.; 674 for (kk=0;kk<NB;kk++) { 675 for (jj=0;jj<NB;jj++) { 676 for (ii=0;ii<NB;ii++) { 677 PetscInt idx = ii + jj*2 + kk*4; 678 for (m=0;m<3;m++) { 679 col[3*idx+m].i = i+ii; 680 col[3*idx+m].j = j+jj; 681 col[3*idx+m].k = k+kk; 682 col[3*idx+m].c = m; 683 if (i+ii >= xs && i+ii < xm+xs && j+jj >= ys && j+jj < ys+ym && k+kk >= zs && k+kk < zs+zm) { 684 row[nrows].i = i+ii; 685 row[nrows].j = j+jj; 686 row[nrows].k = k+kk; 687 row[nrows].c = m; 688 for (l=0;l<NPB;l++) vals[NPB*nrows + l] = ej[NPB*(3*idx+m) + l]; 689 nrows++; 690 } 691 } 692 } 693 } 694 } 695 ierr = MatSetValuesStencil(jac,nrows,row,NPB,col,vals,ADD_VALUES);CHKERRQ(ierr); 696 } 697 } 698 } 699 700 ierr = MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 701 ierr = MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 702 703 /* set the diagonal */ 704 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.; 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 ierr = MatSetValuesStencil(jac,3,col,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 716 } 717 } 718 } 719 } 720 721 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 722 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 723 724 ierr = DMDAVecRestoreArray(cda,C,&c);CHKERRQ(ierr); 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 PetscErrorCode ierr; 740 PetscInt xs=info->xs,ys=info->ys,zs=info->zs; 741 PetscInt xm=info->xm,ym=info->ym,zm=info->zm; 742 PetscInt xes,yes,zes,xee,yee,zee; 743 PetscInt mx=info->mx,my=info->my,mz=info->mz; 744 DM cda; 745 CoordField ***c; 746 Vec C; 747 748 PetscFunctionBegin; 749 ierr = DMGetCoordinateDM(info->da,&cda);CHKERRQ(ierr); 750 ierr = DMGetCoordinatesLocal(info->da,&C);CHKERRQ(ierr); 751 ierr = DMDAVecGetArray(cda,C,&c);CHKERRQ(ierr); 752 ierr = DMDAGetInfo(info->da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 753 ierr = DMDAGetCorners(info->da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 754 755 /* loop over elements */ 756 for (k=zs; k<zs+zm; k++) { 757 for (j=ys; j<ys+ym; j++) { 758 for (i=xs; i<xs+xm; i++) { 759 for (l=0;l<3;l++) { 760 f[k][j][i][l] = 0.; 761 } 762 } 763 } 764 } 765 /* element starts and ends */ 766 xes = xs; 767 yes = ys; 768 zes = zs; 769 xee = xs+xm; 770 yee = ys+ym; 771 zee = zs+zm; 772 if (xs > 0) xes = xs - 1; 773 if (ys > 0) yes = ys - 1; 774 if (zs > 0) zes = zs - 1; 775 if (xs+xm == mx) xee = xs+xm-1; 776 if (ys+ym == my) yee = ys+ym-1; 777 if (zs+zm == mz) zee = zs+zm-1; 778 for (k=zes; k<zee; k++) { 779 for (j=yes; j<yee; j++) { 780 for (i=xes; i<xee; i++) { 781 GatherElementData(mx,my,mz,x,c,i,j,k,ex,ec,user); 782 FormElementJacobian(ex,ec,ef,NULL,user); 783 /* put this element's additions into the residuals */ 784 for (kk=0;kk<NB;kk++) { 785 for (jj=0;jj<NB;jj++) { 786 for (ii=0;ii<NB;ii++) { 787 PetscInt idx = ii + jj*NB + kk*NB*NB; 788 if (k+kk >= zs && j+jj >= ys && i+ii >= xs && k+kk < zs+zm && j+jj < ys+ym && i+ii < xs+xm) { 789 if (OnBoundary(i+ii,j+jj,k+kk,mx,my,mz)) { 790 for (l=0;l<3;l++) 791 f[k+kk][j+jj][i+ii][l] = x[k+kk][j+jj][i+ii][l] - ex[idx][l]; 792 } else { 793 for (l=0;l<3;l++) 794 f[k+kk][j+jj][i+ii][l] += ef[idx][l]; 795 } 796 } 797 } 798 } 799 } 800 } 801 } 802 } 803 ierr = DMDAVecRestoreArray(cda,C,&c);CHKERRQ(ierr); 804 PetscFunctionReturn(0); 805 } 806 807 PetscErrorCode NonlinearGS(SNES snes,Vec X,Vec B,void *ptr) 808 { 809 /* values for each basis function at each quadrature point */ 810 AppCtx *user = (AppCtx*)ptr; 811 PetscInt i,j,k,l,m,n,s; 812 PetscInt pi,pj,pk; 813 Field ef[1]; 814 Field ex[8]; 815 PetscScalar ej[9]; 816 CoordField ec[8]; 817 PetscScalar pjac[9],pjinv[9]; 818 PetscScalar pf[3],py[3]; 819 PetscErrorCode ierr; 820 PetscInt xs,ys,zs; 821 PetscInt xm,ym,zm; 822 PetscInt mx,my,mz; 823 DM cda; 824 CoordField ***c; 825 Vec C; 826 DM da; 827 Vec Xl,Bl; 828 Field ***x,***b; 829 PetscInt sweeps,its; 830 PetscReal atol,rtol,stol; 831 PetscReal fnorm0 = 0.0,fnorm,ynorm,xnorm = 0.0; 832 833 PetscFunctionBegin; 834 ierr = SNESNGSGetSweeps(snes,&sweeps);CHKERRQ(ierr); 835 ierr = SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);CHKERRQ(ierr); 836 837 ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 838 ierr = DMGetLocalVector(da,&Xl);CHKERRQ(ierr); 839 if (B) { 840 ierr = DMGetLocalVector(da,&Bl);CHKERRQ(ierr); 841 } 842 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xl);CHKERRQ(ierr); 843 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xl);CHKERRQ(ierr); 844 if (B) { 845 ierr = DMGlobalToLocalBegin(da,B,INSERT_VALUES,Bl);CHKERRQ(ierr); 846 ierr = DMGlobalToLocalEnd(da,B,INSERT_VALUES,Bl);CHKERRQ(ierr); 847 } 848 ierr = DMDAVecGetArray(da,Xl,&x);CHKERRQ(ierr); 849 if (B) ierr = DMDAVecGetArray(da,Bl,&b);CHKERRQ(ierr); 850 851 ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr); 852 ierr = DMGetCoordinatesLocal(da,&C);CHKERRQ(ierr); 853 ierr = DMDAVecGetArray(cda,C,&c);CHKERRQ(ierr); 854 ierr = DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 855 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 856 857 for (s=0;s<sweeps;s++) { 858 for (k=zs; k<zs+zm; k++) { 859 for (j=ys; j<ys+ym; j++) { 860 for (i=xs; i<xs+xm; i++) { 861 if (OnBoundary(i,j,k,mx,my,mz)) { 862 BoundaryValue(i,j,k,mx,my,mz,x[k][j][i],user); 863 } else { 864 for (n=0;n<its;n++) { 865 for (m=0;m<9;m++) pjac[m] = 0.; 866 for (m=0;m<3;m++) pf[m] = 0.; 867 /* gather the elements for this point */ 868 for (pk=-1; pk<1; pk++) { 869 for (pj=-1; pj<1; pj++) { 870 for (pi=-1; pi<1; pi++) { 871 /* check that this element exists */ 872 if (i+pi >= 0 && i+pi < mx-1 && j+pj >= 0 && j+pj < my-1 && k+pk >= 0 && k+pk < mz-1) { 873 /* create the element function and jacobian */ 874 GatherElementData(mx,my,mz,x,c,i+pi,j+pj,k+pk,ex,ec,user); 875 FormPBJacobian(-pi,-pj,-pk,ex,ec,ef,ej,user); 876 /* extract the point named by i,j,k from the whole element jacobian and function */ 877 for (l=0;l<3;l++) { 878 pf[l] += ef[0][l]; 879 for (m=0;m<3;m++) { 880 pjac[3*m+l] += ej[3*m+l]; 881 } 882 } 883 } 884 } 885 } 886 } 887 /* invert */ 888 InvertTensor(pjac,pjinv,NULL); 889 /* apply */ 890 if (B) for (m=0;m<3;m++) { 891 pf[m] -= b[k][j][i][m]; 892 } 893 TensorVector(pjinv,pf,py); 894 xnorm=0.; 895 for (m=0;m<3;m++) { 896 x[k][j][i][m] -= py[m]; 897 xnorm += PetscRealPart(x[k][j][i][m]*x[k][j][i][m]); 898 } 899 fnorm = PetscRealPart(pf[0]*pf[0]+pf[1]*pf[1]+pf[2]*pf[2]); 900 if (n==0) fnorm0 = fnorm; 901 ynorm = PetscRealPart(py[0]*py[0]+py[1]*py[1]+py[2]*py[2]); 902 if (fnorm < atol*atol || fnorm < rtol*rtol*fnorm0 || ynorm < stol*stol*xnorm) break; 903 } 904 } 905 } 906 } 907 } 908 } 909 ierr = DMDAVecRestoreArray(da,Xl,&x);CHKERRQ(ierr); 910 ierr = DMLocalToGlobalBegin(da,Xl,INSERT_VALUES,X);CHKERRQ(ierr); 911 ierr = DMLocalToGlobalEnd(da,Xl,INSERT_VALUES,X);CHKERRQ(ierr); 912 ierr = DMRestoreLocalVector(da,&Xl);CHKERRQ(ierr); 913 if (B) { 914 ierr = DMDAVecRestoreArray(da,Bl,&b);CHKERRQ(ierr); 915 ierr = DMRestoreLocalVector(da,&Bl);CHKERRQ(ierr); 916 } 917 ierr = DMDAVecRestoreArray(cda,C,&c);CHKERRQ(ierr); 918 PetscFunctionReturn(0); 919 } 920 921 PetscErrorCode FormCoordinates(DM da,AppCtx *user) 922 { 923 PetscErrorCode ierr; 924 Vec coords; 925 DM cda; 926 PetscInt mx,my,mz; 927 PetscInt i,j,k,xs,ys,zs,xm,ym,zm; 928 CoordField ***x; 929 930 PetscFunctionBegin; 931 ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr); 932 ierr = DMCreateGlobalVector(cda,&coords);CHKERRQ(ierr); 933 ierr = DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 934 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 935 ierr = DMDAVecGetArray(da,coords,&x);CHKERRQ(ierr); 936 for (k=zs; k<zs+zm; k++) { 937 for (j=ys; j<ys+ym; j++) { 938 for (i=xs; i<xs+xm; i++) { 939 PetscReal cx = ((PetscReal)i) / (((PetscReal)(mx-1))); 940 PetscReal cy = ((PetscReal)j) / (((PetscReal)(my-1))); 941 PetscReal cz = ((PetscReal)k) / (((PetscReal)(mz-1))); 942 PetscReal rad = user->rad + cy*user->height; 943 PetscReal ang = (cx - 0.5)*user->arc; 944 x[k][j][i][0] = rad*PetscSinReal(ang); 945 x[k][j][i][1] = rad*PetscCosReal(ang) - (user->rad + 0.5*user->height)*PetscCosReal(-0.5*user->arc); 946 x[k][j][i][2] = user->width*(cz - 0.5); 947 } 948 } 949 } 950 ierr = DMDAVecRestoreArray(da,coords,&x);CHKERRQ(ierr); 951 ierr = DMSetCoordinates(da,coords);CHKERRQ(ierr); 952 ierr = VecDestroy(&coords);CHKERRQ(ierr); 953 PetscFunctionReturn(0); 954 } 955 956 PetscErrorCode InitialGuess(DM da,AppCtx *user,Vec X) 957 { 958 PetscInt i,j,k,xs,ys,zs,xm,ym,zm; 959 PetscInt mx,my,mz; 960 PetscErrorCode ierr; 961 Field ***x; 962 963 PetscFunctionBegin; 964 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 965 ierr = DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 966 ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr); 967 968 for (k=zs; k<zs+zm; k++) { 969 for (j=ys; j<ys+ym; j++) { 970 for (i=xs; i<xs+xm; i++) { 971 /* reference coordinates */ 972 PetscReal p_x = ((PetscReal)i) / (((PetscReal)(mx-1))); 973 PetscReal p_y = ((PetscReal)j) / (((PetscReal)(my-1))); 974 PetscReal p_z = ((PetscReal)k) / (((PetscReal)(mz-1))); 975 PetscReal o_x = p_x; 976 PetscReal o_y = p_y; 977 PetscReal o_z = p_z; 978 x[k][j][i][0] = o_x - p_x; 979 x[k][j][i][1] = o_y - p_y; 980 x[k][j][i][2] = o_z - p_z; 981 } 982 } 983 } 984 ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr); 985 PetscFunctionReturn(0); 986 } 987 988 PetscErrorCode FormRHS(DM da,AppCtx *user,Vec X) 989 { 990 PetscInt i,j,k,xs,ys,zs,xm,ym,zm; 991 PetscInt mx,my,mz; 992 PetscErrorCode ierr; 993 Field ***x; 994 995 PetscFunctionBegin; 996 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 997 ierr = DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 998 ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr); 999 1000 for (k=zs; k<zs+zm; k++) { 1001 for (j=ys; j<ys+ym; j++) { 1002 for (i=xs; i<xs+xm; i++) { 1003 x[k][j][i][0] = 0.; 1004 x[k][j][i][1] = 0.; 1005 x[k][j][i][2] = 0.; 1006 if (i == (mx-1)/2 && j == (my-1)) x[k][j][i][1] = user->ploading/(mz-1); 1007 } 1008 } 1009 } 1010 ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr); 1011 PetscFunctionReturn(0); 1012 } 1013 1014 PetscErrorCode DisplayLine(SNES snes,Vec X) 1015 { 1016 PetscInt r,i,j=0,k=0,xs,xm,ys,ym,zs,zm,mx,my,mz; 1017 PetscErrorCode ierr; 1018 Field ***x; 1019 CoordField ***c; 1020 DM da,cda; 1021 Vec C; 1022 PetscMPIInt size,rank; 1023 1024 PetscFunctionBegin; 1025 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 1026 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 1027 ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 1028 ierr = DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1029 ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr); 1030 ierr = DMGetCoordinates(da,&C);CHKERRQ(ierr); 1031 j = my / 2; 1032 k = mz / 2; 1033 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 1034 ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr); 1035 ierr = DMDAVecGetArray(cda,C,&c);CHKERRQ(ierr); 1036 for (r=0;r<size;r++) { 1037 if (rank == r) { 1038 if (j >= ys && j < ys+ym && k >= zs && k < zs+zm) { 1039 for (i=xs; i<xs+xm; i++) { 1040 ierr = 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]));CHKERRQ(ierr); 1041 } 1042 } 1043 } 1044 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 1045 } 1046 ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr); 1047 ierr = DMDAVecRestoreArray(cda,C,&c);CHKERRQ(ierr); 1048 PetscFunctionReturn(0); 1049 } 1050 1051 /*TEST 1052 1053 test: 1054 nsize: 2 1055 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 1056 requires: !single 1057 timeoutfactor: 3 1058 1059 test: 1060 suffix: 2 1061 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 1062 requires: !single 1063 1064 test: 1065 suffix: 3 1066 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 1067 requires: !single 1068 1069 TEST*/ 1070