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