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