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