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