Actual source code: ex56.c

petsc-master 2015-05-28
Report Typos and Errors
  1: static char help[] = "3D, tri-linear quadrilateral (Q1), displacement finite element formulation\n\
  2: of linear elasticity.  E=1.0, nu=0.25.\n\
  3: Unit square domain with Dirichelet boundary condition on the y=0 side only.\n\
  4: Load of 1.0 in x + 2y direction on all nodes (not a true uniform load).\n\
  5:   -ne <size>      : number of (square) quadrilateral elements in each dimension\n\
  6:   -alpha <v>      : scaling of material coeficient in embedded circle\n\n";

  8: #include <petscksp.h>

 10: static PetscBool log_stages = PETSC_TRUE;
 11: static PetscErrorCode MaybeLogStagePush(PetscLogStage stage) { return log_stages ? PetscLogStagePush(stage) : 0; }
 12: static PetscErrorCode MaybeLogStagePop() { return log_stages ? PetscLogStagePop() : 0; }
 13: PetscErrorCode elem_3d_elast_v_25(PetscScalar *);

 17: int main(int argc,char **args)
 18: {
 19:   Mat            Amat;
 21:   PetscInt       m,nn,M,Istart,Iend,i,j,k,ii,jj,kk,ic,ne=4,id;
 22:   PetscReal      x,y,z,h,*coords,soft_alpha=1.e-3;
 23:   PetscBool      two_solves=PETSC_FALSE,test_nonzero_cols=PETSC_FALSE,use_nearnullspace=PETSC_FALSE;
 24:   Vec            xx,bb;
 25:   KSP            ksp;
 26:   MPI_Comm       comm;
 27:   PetscMPIInt    npe,mype;
 28:   PC             pc;
 29:   PetscScalar    DD[24][24],DD2[24][24];
 30:   PetscLogStage  stage[6];
 31:   PetscScalar    DD1[24][24];
 32:   PCType         type;

 34:   PetscInitialize(&argc,&args,(char*)0,help);
 35:   comm = PETSC_COMM_WORLD;
 36:   MPI_Comm_rank(comm, &mype);
 37:   MPI_Comm_size(comm, &npe);

 39:   PetscOptionsBegin(comm,NULL,"3D bilinear Q1 elasticity options","");
 40:   {
 41:     char nestring[256];
 42:     PetscSNPrintf(nestring,sizeof nestring,"number of elements in each direction, ne+1 must be a multiple of %D (sizes^{1/3})",(PetscInt)(PetscPowReal((PetscReal)npe,1./3.) + .5));
 43:     PetscOptionsInt("-ne",nestring,"",ne,&ne,NULL);
 44:     PetscOptionsBool("-log_stages","Log stages of solve separately","",log_stages,&log_stages,NULL);
 45:     PetscOptionsReal("-alpha","material coefficient inside circle","",soft_alpha,&soft_alpha,NULL);
 46:     PetscOptionsBool("-two_solves","solve additional variant of the problem","",two_solves,&two_solves,NULL);
 47:     PetscOptionsBool("-test_nonzero_cols","nonzero test","",test_nonzero_cols,&test_nonzero_cols,NULL);
 48:     PetscOptionsBool("-use_mat_nearnullspace","MatNearNullSpace API test","",use_nearnullspace,&use_nearnullspace,NULL);
 49:   }
 50:   PetscOptionsEnd();

 52:   if (log_stages) {
 53:     PetscLogStageRegister("Setup", &stage[0]);
 54:     PetscLogStageRegister("Solve", &stage[1]);
 55:     PetscLogStageRegister("2nd Setup", &stage[2]);
 56:     PetscLogStageRegister("2nd Solve", &stage[3]);
 57:     PetscLogStageRegister("3rd Setup", &stage[4]);
 58:     PetscLogStageRegister("3rd Solve", &stage[5]);
 59:   } else {
 60:     for (i=0; i<(PetscInt)(sizeof(stage)/sizeof(stage[0])); i++) stage[i] = -1;
 61:   }

 63:   h = 1./ne; nn = ne+1;
 64:   /* ne*ne; number of global elements */
 65:   M = 3*nn*nn*nn; /* global number of equations */
 66:   if (npe==2) {
 67:     if (mype==1) m=0;
 68:     else m = nn*nn*nn;
 69:     npe = 1;
 70:   } else {
 71:     m = nn*nn*nn/npe;
 72:     if (mype==npe-1) m = nn*nn*nn - (npe-1)*m;
 73:   }
 74:   m *= 3; /* number of equations local*/
 75:   /* Setup solver, get PC type and pc */
 76:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 77:   KSPSetType(ksp, KSPCG);
 78:   KSPSetComputeSingularValues(ksp, PETSC_TRUE);
 79:   KSPGetPC(ksp, &pc);
 80:   PCSetType(pc, PCGAMG); /* default */
 81:   KSPSetFromOptions(ksp);
 82:   PCGetType(pc, &type);

 84:   {
 85:     /* configureation */
 86:     const PetscInt NP = (PetscInt)(PetscPowReal((PetscReal)npe,1./3.) + .5);
 87:     const PetscInt ipx = mype%NP, ipy = (mype%(NP*NP))/NP, ipz = mype/(NP*NP);
 88:     const PetscInt Ni0 = ipx*(nn/NP), Nj0 = ipy*(nn/NP), Nk0 = ipz*(nn/NP);
 89:     const PetscInt Ni1 = Ni0 + (m>0 ? (nn/NP) : 0), Nj1 = Nj0 + (nn/NP), Nk1 = Nk0 + (nn/NP);
 90:     const PetscInt NN  = nn/NP, id0 = ipz*nn*nn*NN + ipy*nn*NN*NN + ipx*NN*NN*NN;
 91:     PetscInt       *d_nnz, *o_nnz,osz[4]={0,9,15,19},nbc;
 92:     PetscScalar    vv[24], v2[24];
 93:     if (npe!=NP*NP*NP) SETERRQ1(comm,PETSC_ERR_ARG_WRONG, "npe=%d: npe^{1/3} must be integer",npe);
 94:     if (nn!=NP*(nn/NP)) SETERRQ1(comm,PETSC_ERR_ARG_WRONG, "-ne %d: (ne+1)%(npe^{1/3}) must equal zero",ne);

 96:     /* count nnz */
 97:     PetscMalloc1(m+1, &d_nnz);
 98:     PetscMalloc1(m+1, &o_nnz);
 99:     for (i=Ni0,ic=0; i<Ni1; i++) {
100:       for (j=Nj0; j<Nj1; j++) {
101:         for (k=Nk0; k<Nk1; k++) {
102:           nbc = 0;
103:           if (i==Ni0 || i==Ni1-1) nbc++;
104:           if (j==Nj0 || j==Nj1-1) nbc++;
105:           if (k==Nk0 || k==Nk1-1) nbc++;
106:           for (jj=0; jj<3; jj++,ic++) {
107:             d_nnz[ic] = 3*(27-osz[nbc]);
108:             o_nnz[ic] = 3*osz[nbc];
109:           }
110:         }
111:       }
112:     }
113:     if (ic != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ic %D does not equal m %D",ic,m);

115:     /* create stiffness matrix */
116:     MatCreate(comm,&Amat);
117:     MatSetSizes(Amat,m,m,M,M);
118:     MatSetBlockSize(Amat,3);
119:     MatSetType(Amat,MATAIJ);
120:     MatSeqAIJSetPreallocation(Amat,0,d_nnz);
121:     MatMPIAIJSetPreallocation(Amat,0,d_nnz,0,o_nnz);

123:     PetscFree(d_nnz);
124:     PetscFree(o_nnz);

126:     MatGetOwnershipRange(Amat,&Istart,&Iend);

128:     if (m != Iend - Istart) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"m %D does not equal Iend %D - Istart %D",m,Iend,Istart);
129:     /* Generate vectors */
130:     VecCreate(comm,&xx);
131:     VecSetSizes(xx,m,M);
132:     VecSetBlockSize(xx,3);
133:     VecSetFromOptions(xx);
134:     VecDuplicate(xx,&bb);
135:     VecSet(bb,.0);
136:     /* generate element matrices */
137:     {
138:       PetscBool hasData = PETSC_TRUE;
139:       if (!hasData) {
140:         PetscPrintf(PETSC_COMM_WORLD,"\t No data is provided\n");
141:         for (i=0; i<24; i++) {
142:           for (j=0; j<24; j++) {
143:             if (i==j) DD1[i][j] = 1.0;
144:             else DD1[i][j] = -.25;
145:           }
146:         }
147:       } else {
148:         /* Get array data */
149:         elem_3d_elast_v_25((PetscScalar*)DD1);
150:       }
151: 
152:       /* BC version of element */
153:       for (i=0; i<24; i++) {
154:         for (j=0; j<24; j++) {
155:           if (i<12 || (j < 12 && !test_nonzero_cols)) {
156:             if (i==j) DD2[i][j] = 0.1*DD1[i][j];
157:             else DD2[i][j] = 0.0;
158:           } else DD2[i][j] = DD1[i][j];
159:         }
160:       }
161:       /* element residual/load vector */
162:       for (i=0; i<24; i++) {
163:         if (i%3==0) vv[i] = h*h;
164:         else if (i%3==1) vv[i] = 2.0*h*h;
165:         else vv[i] = .0;
166:       }
167:       for (i=0; i<24; i++) {
168:         if (i%3==0 && i>=12) v2[i] = h*h;
169:         else if (i%3==1 && i>=12) v2[i] = 2.0*h*h;
170:         else v2[i] = .0;
171:       }
172:     }

174:     PetscMalloc1(m+1, &coords);
175:     coords[m] = -99.0;

177:     /* forms the element stiffness and coordinates */
178:     for (i=Ni0,ic=0,ii=0; i<Ni1; i++,ii++) {
179:       for (j=Nj0,jj=0; j<Nj1; j++,jj++) {
180:         for (k=Nk0,kk=0; k<Nk1; k++,kk++,ic++) {

182:           /* coords */
183:           x = coords[3*ic] = h*(PetscReal)i;
184:           y = coords[3*ic+1] = h*(PetscReal)j;
185:           z = coords[3*ic+2] = h*(PetscReal)k;
186:           /* matrix */
187:           id = id0 + ii + NN*jj + NN*NN*kk;

189:           if (i<ne && j<ne && k<ne) {
190:             /* radius */
191:             PetscReal radius = PetscSqrtReal((x-.5+h/2)*(x-.5+h/2)+(y-.5+h/2)*(y-.5+h/2)+(z-.5+h/2)*(z-.5+h/2));
192:             PetscReal alpha = 1.0;
193:             PetscInt  jx,ix,idx[8] = { id, id+1, id+NN+1, id+NN,
194:                                        id        + NN*NN, id+1    + NN*NN,
195:                                        id+NN+1 + NN*NN, id+NN + NN*NN };

197:             /* correct indices */
198:             if (i==Ni1-1 && Ni1!=nn) {
199:               idx[1] += NN*(NN*NN-1);
200:               idx[2] += NN*(NN*NN-1);
201:               idx[5] += NN*(NN*NN-1);
202:               idx[6] += NN*(NN*NN-1);
203:             }
204:             if (j==Nj1-1 && Nj1!=nn) {
205:               idx[2] += NN*NN*(nn-1);
206:               idx[3] += NN*NN*(nn-1);
207:               idx[6] += NN*NN*(nn-1);
208:               idx[7] += NN*NN*(nn-1);
209:             }
210:             if (k==Nk1-1 && Nk1!=nn) {
211:               idx[4] += NN*(nn*nn-NN*NN);
212:               idx[5] += NN*(nn*nn-NN*NN);
213:               idx[6] += NN*(nn*nn-NN*NN);
214:               idx[7] += NN*(nn*nn-NN*NN);
215:             }

217:             if (radius < 0.25) alpha = soft_alpha;

219:             for (ix=0; ix<24; ix++) {
220:               for (jx=0;jx<24;jx++) DD[ix][jx] = alpha*DD1[ix][jx];
221:             }
222:             if (k>0) {
223:               MatSetValuesBlocked(Amat,8,idx,8,idx,(const PetscScalar*)DD,ADD_VALUES);
224:               VecSetValuesBlocked(bb,8,idx,(const PetscScalar*)vv,ADD_VALUES);
225:             } else {
226:               /* a BC */
227:               for (ix=0;ix<24;ix++) {
228:                 for (jx=0;jx<24;jx++) DD[ix][jx] = alpha*DD2[ix][jx];
229:               }
230:               MatSetValuesBlocked(Amat,8,idx,8,idx,(const PetscScalar*)DD,ADD_VALUES);
231:               VecSetValuesBlocked(bb,8,idx,(const PetscScalar*)v2,ADD_VALUES);
232:             }
233:           }
234:         }
235:       }

237:     }
238:     MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
239:     MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
240:     VecAssemblyBegin(bb);
241:     VecAssemblyEnd(bb);
242:   }

244:   if (!PETSC_TRUE) {
245:     PetscViewer viewer;
246:     PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
247:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
248:     MatView(Amat,viewer);
249:     PetscViewerDestroy(&viewer);
250:   }

252:   /* finish KSP/PC setup */
253:   KSPSetOperators(ksp, Amat, Amat);
254:   if (use_nearnullspace) {
255:     MatNullSpace matnull;
256:     Vec          vec_coords;
257:     PetscScalar  *c;

259:     VecCreate(MPI_COMM_WORLD,&vec_coords);
260:     VecSetBlockSize(vec_coords,3);
261:     VecSetSizes(vec_coords,m,PETSC_DECIDE);
262:     VecSetUp(vec_coords);
263:     VecGetArray(vec_coords,&c);
264:     for (i=0; i<m; i++) c[i] = coords[i]; /* Copy since Scalar type might be Complex */
265:     VecRestoreArray(vec_coords,&c);
266:     MatNullSpaceCreateRigidBody(vec_coords,&matnull);
267:     MatSetNearNullSpace(Amat,matnull);
268:     MatNullSpaceDestroy(&matnull);
269:     VecDestroy(&vec_coords);
270:   } else {
271:     PCSetCoordinates(pc, 3, m/3, coords);
272:   }

274:   MaybeLogStagePush(stage[0]);

276:   /* PC setup basically */
277:   KSPSetUp(ksp);

279:   MaybeLogStagePop();
280:   MaybeLogStagePush(stage[1]);

282:   /* test BCs */
283:   if (test_nonzero_cols) {
284:     VecZeroEntries(xx);
285:     if (mype==0) VecSetValue(xx,0,1.0,INSERT_VALUES);
286:     VecAssemblyBegin(xx);
287:     VecAssemblyEnd(xx);
288:     KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
289:   }

291:   /* 1st solve */
292:   KSPSolve(ksp, bb, xx);

294:   MaybeLogStagePop();

296:   /* 2nd solve */
297:   if (two_solves) {
298:     PetscReal emax, emin;
299:     PetscReal norm,norm2;
300:     Vec       res;

302:     MaybeLogStagePush(stage[2]);
303:     /* PC setup basically */
304:     MatScale(Amat, 100000.0);
305:     KSPSetOperators(ksp, Amat, Amat);
306:     KSPSetUp(ksp);

308:     MaybeLogStagePop();
309:     MaybeLogStagePush(stage[3]);
310:     KSPSolve(ksp, bb, xx);
311:     KSPComputeExtremeSingularValues(ksp, &emax, &emin);

313:     MaybeLogStagePop();
314:     MaybeLogStagePush(stage[4]);

316:     /* 3rd solve */
317:     MatScale(Amat, 100000.0);
318:     KSPSetOperators(ksp, Amat, Amat);
319:     KSPSetUp(ksp);

321:     MaybeLogStagePop();
322:     MaybeLogStagePush(stage[5]);

324:     KSPSolve(ksp, bb, xx);

326:     MaybeLogStagePop();


329:     VecNorm(bb, NORM_2, &norm2);

331:     VecDuplicate(xx, &res);
332:     MatMult(Amat, xx, res);
333:     VecAXPY(bb, -1.0, res);
334:     VecDestroy(&res);
335:     VecNorm(bb, NORM_2, &norm);
336:     PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e, emax=%e\n",0,__FUNCT__,(double)(norm/norm2),(double)norm2,(double)emax);
337:     /*PetscViewerASCIIOpen(comm, "residual.m", &viewer);
338:      PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
339:      VecView(bb,viewer);
340:      PetscViewerDestroy(&viewer);*/


343:     /* PetscViewerASCIIOpen(comm, "rhs.m", &viewer); */
344:     /* PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB); */
345:     /*  */
346:     /* VecView(bb,viewer); */
347:     /* PetscViewerDestroy(&viewer); */

349:     /* PetscViewerASCIIOpen(comm, "solution.m", &viewer); */
350:     /* PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB); */
351:     /*  */
352:     /* VecView(xx, viewer); */
353:     /* PetscViewerDestroy(&viewer); */
354:   }

356:   /* Free work space */
357:   KSPDestroy(&ksp);
358:   VecDestroy(&xx);
359:   VecDestroy(&bb);
360:   MatDestroy(&Amat);
361:   PetscFree(coords);

363:   PetscFinalize();
364:   return 0;
365: }

367: /* Data was previously provided in the file data/elem_3d_elast_v_25.tx */
370: PetscErrorCode elem_3d_elast_v_25(PetscScalar *dd)
371: {
373:   PetscScalar    DD[] = {
374:   0.18981481481481474     ,
375:   5.27777777777777568E-002,
376:   5.27777777777777568E-002,
377:  -5.64814814814814659E-002,
378:  -1.38888888888889072E-002,
379:  -1.38888888888889089E-002,
380:  -8.24074074074073876E-002,
381:  -5.27777777777777429E-002,
382:   1.38888888888888725E-002,
383:   4.90740740740740339E-002,
384:   1.38888888888889124E-002,
385:   4.72222222222222071E-002,
386:   4.90740740740740339E-002,
387:   4.72222222222221932E-002,
388:   1.38888888888888968E-002,
389:  -8.24074074074073876E-002,
390:   1.38888888888888673E-002,
391:  -5.27777777777777429E-002,
392:  -7.87037037037036785E-002,
393:  -4.72222222222221932E-002,
394:  -4.72222222222222071E-002,
395:   1.20370370370370180E-002,
396:  -1.38888888888888742E-002,
397:  -1.38888888888888829E-002,
398:   5.27777777777777568E-002,
399:   0.18981481481481474     ,
400:   5.27777777777777568E-002,
401:   1.38888888888889124E-002,
402:   4.90740740740740269E-002,
403:   4.72222222222221932E-002,
404:  -5.27777777777777637E-002,
405:  -8.24074074074073876E-002,
406:   1.38888888888888725E-002,
407:  -1.38888888888889037E-002,
408:  -5.64814814814814728E-002,
409:  -1.38888888888888985E-002,
410:   4.72222222222221932E-002,
411:   4.90740740740740478E-002,
412:   1.38888888888888968E-002,
413:  -1.38888888888888673E-002,
414:   1.20370370370370058E-002,
415:  -1.38888888888888742E-002,
416:  -4.72222222222221932E-002,
417:  -7.87037037037036785E-002,
418:  -4.72222222222222002E-002,
419:   1.38888888888888742E-002,
420:  -8.24074074074073598E-002,
421:  -5.27777777777777568E-002,
422:   5.27777777777777568E-002,
423:   5.27777777777777568E-002,
424:   0.18981481481481474     ,
425:   1.38888888888889055E-002,
426:   4.72222222222222002E-002,
427:   4.90740740740740269E-002,
428:  -1.38888888888888829E-002,
429:  -1.38888888888888829E-002,
430:   1.20370370370370180E-002,
431:   4.72222222222222002E-002,
432:   1.38888888888888985E-002,
433:   4.90740740740740339E-002,
434:  -1.38888888888888985E-002,
435:  -1.38888888888888968E-002,
436:  -5.64814814814814520E-002,
437:  -5.27777777777777568E-002,
438:   1.38888888888888777E-002,
439:  -8.24074074074073876E-002,
440:  -4.72222222222222002E-002,
441:  -4.72222222222221932E-002,
442:  -7.87037037037036646E-002,
443:   1.38888888888888794E-002,
444:  -5.27777777777777568E-002,
445:  -8.24074074074073598E-002,
446:  -5.64814814814814659E-002,
447:   1.38888888888889124E-002,
448:   1.38888888888889055E-002,
449:   0.18981481481481474     ,
450:  -5.27777777777777568E-002,
451:  -5.27777777777777499E-002,
452:   4.90740740740740269E-002,
453:  -1.38888888888889072E-002,
454:  -4.72222222222221932E-002,
455:  -8.24074074074073876E-002,
456:   5.27777777777777568E-002,
457:  -1.38888888888888812E-002,
458:  -8.24074074074073876E-002,
459:  -1.38888888888888742E-002,
460:   5.27777777777777499E-002,
461:   4.90740740740740269E-002,
462:  -4.72222222222221863E-002,
463:  -1.38888888888889089E-002,
464:   1.20370370370370162E-002,
465:   1.38888888888888673E-002,
466:   1.38888888888888742E-002,
467:  -7.87037037037036785E-002,
468:   4.72222222222222002E-002,
469:   4.72222222222222071E-002,
470:  -1.38888888888889072E-002,
471:   4.90740740740740269E-002,
472:   4.72222222222222002E-002,
473:  -5.27777777777777568E-002,
474:   0.18981481481481480     ,
475:   5.27777777777777568E-002,
476:   1.38888888888889020E-002,
477:  -5.64814814814814728E-002,
478:  -1.38888888888888951E-002,
479:   5.27777777777777637E-002,
480:  -8.24074074074073876E-002,
481:   1.38888888888888881E-002,
482:   1.38888888888888742E-002,
483:   1.20370370370370232E-002,
484:  -1.38888888888888812E-002,
485:  -4.72222222222221863E-002,
486:   4.90740740740740339E-002,
487:   1.38888888888888933E-002,
488:  -1.38888888888888812E-002,
489:  -8.24074074074073876E-002,
490:  -5.27777777777777568E-002,
491:   4.72222222222222071E-002,
492:  -7.87037037037036924E-002,
493:  -4.72222222222222140E-002,
494:  -1.38888888888889089E-002,
495:   4.72222222222221932E-002,
496:   4.90740740740740269E-002,
497:  -5.27777777777777499E-002,
498:   5.27777777777777568E-002,
499:   0.18981481481481477     ,
500:  -4.72222222222222071E-002,
501:   1.38888888888888968E-002,
502:   4.90740740740740131E-002,
503:   1.38888888888888812E-002,
504:  -1.38888888888888708E-002,
505:   1.20370370370370267E-002,
506:   5.27777777777777568E-002,
507:   1.38888888888888812E-002,
508:  -8.24074074074073876E-002,
509:   1.38888888888889124E-002,
510:  -1.38888888888889055E-002,
511:  -5.64814814814814589E-002,
512:  -1.38888888888888812E-002,
513:  -5.27777777777777568E-002,
514:  -8.24074074074073737E-002,
515:   4.72222222222222002E-002,
516:  -4.72222222222222002E-002,
517:  -7.87037037037036924E-002,
518:  -8.24074074074073876E-002,
519:  -5.27777777777777637E-002,
520:  -1.38888888888888829E-002,
521:   4.90740740740740269E-002,
522:   1.38888888888889020E-002,
523:  -4.72222222222222071E-002,
524:   0.18981481481481480     ,
525:   5.27777777777777637E-002,
526:  -5.27777777777777637E-002,
527:  -5.64814814814814728E-002,
528:  -1.38888888888889037E-002,
529:   1.38888888888888951E-002,
530:  -7.87037037037036785E-002,
531:  -4.72222222222222002E-002,
532:   4.72222222222221932E-002,
533:   1.20370370370370128E-002,
534:  -1.38888888888888725E-002,
535:   1.38888888888888812E-002,
536:   4.90740740740740408E-002,
537:   4.72222222222222002E-002,
538:  -1.38888888888888951E-002,
539:  -8.24074074074073876E-002,
540:   1.38888888888888812E-002,
541:   5.27777777777777637E-002,
542:  -5.27777777777777429E-002,
543:  -8.24074074074073876E-002,
544:  -1.38888888888888829E-002,
545:  -1.38888888888889072E-002,
546:  -5.64814814814814728E-002,
547:   1.38888888888888968E-002,
548:   5.27777777777777637E-002,
549:   0.18981481481481480     ,
550:  -5.27777777777777568E-002,
551:   1.38888888888888916E-002,
552:   4.90740740740740339E-002,
553:  -4.72222222222222210E-002,
554:  -4.72222222222221932E-002,
555:  -7.87037037037036924E-002,
556:   4.72222222222222002E-002,
557:   1.38888888888888742E-002,
558:  -8.24074074074073876E-002,
559:   5.27777777777777429E-002,
560:   4.72222222222222002E-002,
561:   4.90740740740740269E-002,
562:  -1.38888888888888951E-002,
563:  -1.38888888888888846E-002,
564:   1.20370370370370267E-002,
565:   1.38888888888888916E-002,
566:   1.38888888888888725E-002,
567:   1.38888888888888725E-002,
568:   1.20370370370370180E-002,
569:  -4.72222222222221932E-002,
570:  -1.38888888888888951E-002,
571:   4.90740740740740131E-002,
572:  -5.27777777777777637E-002,
573:  -5.27777777777777568E-002,
574:   0.18981481481481480     ,
575:  -1.38888888888888968E-002,
576:  -4.72222222222221932E-002,
577:   4.90740740740740339E-002,
578:   4.72222222222221932E-002,
579:   4.72222222222222071E-002,
580:  -7.87037037037036646E-002,
581:  -1.38888888888888742E-002,
582:   5.27777777777777499E-002,
583:  -8.24074074074073737E-002,
584:   1.38888888888888933E-002,
585:   1.38888888888889020E-002,
586:  -5.64814814814814589E-002,
587:   5.27777777777777568E-002,
588:  -1.38888888888888794E-002,
589:  -8.24074074074073876E-002,
590:   4.90740740740740339E-002,
591:  -1.38888888888889037E-002,
592:   4.72222222222222002E-002,
593:  -8.24074074074073876E-002,
594:   5.27777777777777637E-002,
595:   1.38888888888888812E-002,
596:  -5.64814814814814728E-002,
597:   1.38888888888888916E-002,
598:  -1.38888888888888968E-002,
599:   0.18981481481481480     ,
600:  -5.27777777777777499E-002,
601:   5.27777777777777707E-002,
602:   1.20370370370370180E-002,
603:   1.38888888888888812E-002,
604:  -1.38888888888888812E-002,
605:  -7.87037037037036785E-002,
606:   4.72222222222222002E-002,
607:  -4.72222222222222071E-002,
608:  -8.24074074074073876E-002,
609:  -1.38888888888888742E-002,
610:  -5.27777777777777568E-002,
611:   4.90740740740740616E-002,
612:  -4.72222222222222002E-002,
613:   1.38888888888888846E-002,
614:   1.38888888888889124E-002,
615:  -5.64814814814814728E-002,
616:   1.38888888888888985E-002,
617:   5.27777777777777568E-002,
618:  -8.24074074074073876E-002,
619:  -1.38888888888888708E-002,
620:  -1.38888888888889037E-002,
621:   4.90740740740740339E-002,
622:  -4.72222222222221932E-002,
623:  -5.27777777777777499E-002,
624:   0.18981481481481480     ,
625:  -5.27777777777777568E-002,
626:  -1.38888888888888673E-002,
627:  -8.24074074074073598E-002,
628:   5.27777777777777429E-002,
629:   4.72222222222222002E-002,
630:  -7.87037037037036785E-002,
631:   4.72222222222222002E-002,
632:   1.38888888888888708E-002,
633:   1.20370370370370128E-002,
634:   1.38888888888888760E-002,
635:  -4.72222222222222002E-002,
636:   4.90740740740740478E-002,
637:  -1.38888888888888951E-002,
638:   4.72222222222222071E-002,
639:  -1.38888888888888985E-002,
640:   4.90740740740740339E-002,
641:  -1.38888888888888812E-002,
642:   1.38888888888888881E-002,
643:   1.20370370370370267E-002,
644:   1.38888888888888951E-002,
645:  -4.72222222222222210E-002,
646:   4.90740740740740339E-002,
647:   5.27777777777777707E-002,
648:  -5.27777777777777568E-002,
649:   0.18981481481481477     ,
650:   1.38888888888888829E-002,
651:   5.27777777777777707E-002,
652:  -8.24074074074073598E-002,
653:  -4.72222222222222140E-002,
654:   4.72222222222222140E-002,
655:  -7.87037037037036646E-002,
656:  -5.27777777777777707E-002,
657:  -1.38888888888888829E-002,
658:  -8.24074074074073876E-002,
659:  -1.38888888888888881E-002,
660:   1.38888888888888881E-002,
661:  -5.64814814814814589E-002,
662:   4.90740740740740339E-002,
663:   4.72222222222221932E-002,
664:  -1.38888888888888985E-002,
665:  -8.24074074074073876E-002,
666:   1.38888888888888742E-002,
667:   5.27777777777777568E-002,
668:  -7.87037037037036785E-002,
669:  -4.72222222222221932E-002,
670:   4.72222222222221932E-002,
671:   1.20370370370370180E-002,
672:  -1.38888888888888673E-002,
673:   1.38888888888888829E-002,
674:   0.18981481481481469     ,
675:   5.27777777777777429E-002,
676:  -5.27777777777777429E-002,
677:  -5.64814814814814659E-002,
678:  -1.38888888888889055E-002,
679:   1.38888888888889055E-002,
680:  -8.24074074074074153E-002,
681:  -5.27777777777777429E-002,
682:  -1.38888888888888760E-002,
683:   4.90740740740740408E-002,
684:   1.38888888888888968E-002,
685:  -4.72222222222222071E-002,
686:   4.72222222222221932E-002,
687:   4.90740740740740478E-002,
688:  -1.38888888888888968E-002,
689:  -1.38888888888888742E-002,
690:   1.20370370370370232E-002,
691:   1.38888888888888812E-002,
692:  -4.72222222222222002E-002,
693:  -7.87037037037036924E-002,
694:   4.72222222222222071E-002,
695:   1.38888888888888812E-002,
696:  -8.24074074074073598E-002,
697:   5.27777777777777707E-002,
698:   5.27777777777777429E-002,
699:   0.18981481481481477     ,
700:  -5.27777777777777499E-002,
701:   1.38888888888889107E-002,
702:   4.90740740740740478E-002,
703:  -4.72222222222221932E-002,
704:  -5.27777777777777568E-002,
705:  -8.24074074074074153E-002,
706:  -1.38888888888888812E-002,
707:  -1.38888888888888846E-002,
708:  -5.64814814814814659E-002,
709:   1.38888888888888812E-002,
710:   1.38888888888888968E-002,
711:   1.38888888888888968E-002,
712:  -5.64814814814814520E-002,
713:   5.27777777777777499E-002,
714:  -1.38888888888888812E-002,
715:  -8.24074074074073876E-002,
716:   4.72222222222221932E-002,
717:   4.72222222222222002E-002,
718:  -7.87037037037036646E-002,
719:  -1.38888888888888812E-002,
720:   5.27777777777777429E-002,
721:  -8.24074074074073598E-002,
722:  -5.27777777777777429E-002,
723:  -5.27777777777777499E-002,
724:   0.18981481481481474     ,
725:  -1.38888888888888985E-002,
726:  -4.72222222222221863E-002,
727:   4.90740740740740339E-002,
728:   1.38888888888888829E-002,
729:   1.38888888888888777E-002,
730:   1.20370370370370249E-002,
731:  -4.72222222222222002E-002,
732:  -1.38888888888888933E-002,
733:   4.90740740740740339E-002,
734:  -8.24074074074073876E-002,
735:  -1.38888888888888673E-002,
736:  -5.27777777777777568E-002,
737:   4.90740740740740269E-002,
738:  -4.72222222222221863E-002,
739:   1.38888888888889124E-002,
740:   1.20370370370370128E-002,
741:   1.38888888888888742E-002,
742:  -1.38888888888888742E-002,
743:  -7.87037037037036785E-002,
744:   4.72222222222222002E-002,
745:  -4.72222222222222140E-002,
746:  -5.64814814814814659E-002,
747:   1.38888888888889107E-002,
748:  -1.38888888888888985E-002,
749:   0.18981481481481474     ,
750:  -5.27777777777777499E-002,
751:   5.27777777777777499E-002,
752:   4.90740740740740339E-002,
753:  -1.38888888888889055E-002,
754:   4.72222222222221932E-002,
755:  -8.24074074074074153E-002,
756:   5.27777777777777499E-002,
757:   1.38888888888888829E-002,
758:   1.38888888888888673E-002,
759:   1.20370370370370058E-002,
760:   1.38888888888888777E-002,
761:  -4.72222222222221863E-002,
762:   4.90740740740740339E-002,
763:  -1.38888888888889055E-002,
764:  -1.38888888888888725E-002,
765:  -8.24074074074073876E-002,
766:   5.27777777777777499E-002,
767:   4.72222222222222002E-002,
768:  -7.87037037037036785E-002,
769:   4.72222222222222140E-002,
770:  -1.38888888888889055E-002,
771:   4.90740740740740478E-002,
772:  -4.72222222222221863E-002,
773:  -5.27777777777777499E-002,
774:   0.18981481481481469     ,
775:  -5.27777777777777499E-002,
776:   1.38888888888889072E-002,
777:  -5.64814814814814659E-002,
778:   1.38888888888889003E-002,
779:   5.27777777777777429E-002,
780:  -8.24074074074074153E-002,
781:  -1.38888888888888812E-002,
782:  -5.27777777777777429E-002,
783:  -1.38888888888888742E-002,
784:  -8.24074074074073876E-002,
785:  -1.38888888888889089E-002,
786:   1.38888888888888933E-002,
787:  -5.64814814814814589E-002,
788:   1.38888888888888812E-002,
789:   5.27777777777777429E-002,
790:  -8.24074074074073737E-002,
791:  -4.72222222222222071E-002,
792:   4.72222222222222002E-002,
793:  -7.87037037037036646E-002,
794:   1.38888888888889055E-002,
795:  -4.72222222222221932E-002,
796:   4.90740740740740339E-002,
797:   5.27777777777777499E-002,
798:  -5.27777777777777499E-002,
799:   0.18981481481481474     ,
800:   4.72222222222222002E-002,
801:  -1.38888888888888985E-002,
802:   4.90740740740740339E-002,
803:  -1.38888888888888846E-002,
804:   1.38888888888888812E-002,
805:   1.20370370370370284E-002,
806:  -7.87037037037036785E-002,
807:  -4.72222222222221932E-002,
808:  -4.72222222222222002E-002,
809:   1.20370370370370162E-002,
810:  -1.38888888888888812E-002,
811:  -1.38888888888888812E-002,
812:   4.90740740740740408E-002,
813:   4.72222222222222002E-002,
814:   1.38888888888888933E-002,
815:  -8.24074074074073876E-002,
816:   1.38888888888888708E-002,
817:  -5.27777777777777707E-002,
818:  -8.24074074074074153E-002,
819:  -5.27777777777777568E-002,
820:   1.38888888888888829E-002,
821:   4.90740740740740339E-002,
822:   1.38888888888889072E-002,
823:   4.72222222222222002E-002,
824:   0.18981481481481477     ,
825:   5.27777777777777429E-002,
826:   5.27777777777777568E-002,
827:  -5.64814814814814659E-002,
828:  -1.38888888888888846E-002,
829:  -1.38888888888888881E-002,
830:  -4.72222222222221932E-002,
831:  -7.87037037037036785E-002,
832:  -4.72222222222221932E-002,
833:   1.38888888888888673E-002,
834:  -8.24074074074073876E-002,
835:  -5.27777777777777568E-002,
836:   4.72222222222222002E-002,
837:   4.90740740740740269E-002,
838:   1.38888888888889020E-002,
839:  -1.38888888888888742E-002,
840:   1.20370370370370128E-002,
841:  -1.38888888888888829E-002,
842:  -5.27777777777777429E-002,
843:  -8.24074074074074153E-002,
844:   1.38888888888888777E-002,
845:  -1.38888888888889055E-002,
846:  -5.64814814814814659E-002,
847:  -1.38888888888888985E-002,
848:   5.27777777777777429E-002,
849:   0.18981481481481469     ,
850:   5.27777777777777429E-002,
851:   1.38888888888888933E-002,
852:   4.90740740740740339E-002,
853:   4.72222222222222071E-002,
854:  -4.72222222222222071E-002,
855:  -4.72222222222222002E-002,
856:  -7.87037037037036646E-002,
857:   1.38888888888888742E-002,
858:  -5.27777777777777568E-002,
859:  -8.24074074074073737E-002,
860:  -1.38888888888888951E-002,
861:  -1.38888888888888951E-002,
862:  -5.64814814814814589E-002,
863:  -5.27777777777777568E-002,
864:   1.38888888888888760E-002,
865:  -8.24074074074073876E-002,
866:  -1.38888888888888760E-002,
867:  -1.38888888888888812E-002,
868:   1.20370370370370249E-002,
869:   4.72222222222221932E-002,
870:   1.38888888888889003E-002,
871:   4.90740740740740339E-002,
872:   5.27777777777777568E-002,
873:   5.27777777777777429E-002,
874:   0.18981481481481474     ,
875:   1.38888888888888933E-002,
876:   4.72222222222222071E-002,
877:   4.90740740740740339E-002,
878:   1.20370370370370180E-002,
879:   1.38888888888888742E-002,
880:   1.38888888888888794E-002,
881:  -7.87037037037036785E-002,
882:   4.72222222222222071E-002,
883:   4.72222222222222002E-002,
884:  -8.24074074074073876E-002,
885:  -1.38888888888888846E-002,
886:   5.27777777777777568E-002,
887:   4.90740740740740616E-002,
888:  -4.72222222222222002E-002,
889:  -1.38888888888888881E-002,
890:   4.90740740740740408E-002,
891:  -1.38888888888888846E-002,
892:  -4.72222222222222002E-002,
893:  -8.24074074074074153E-002,
894:   5.27777777777777429E-002,
895:  -1.38888888888888846E-002,
896:  -5.64814814814814659E-002,
897:   1.38888888888888933E-002,
898:   1.38888888888888933E-002,
899:   0.18981481481481477     ,
900:  -5.27777777777777568E-002,
901:  -5.27777777777777637E-002,
902:  -1.38888888888888742E-002,
903:  -8.24074074074073598E-002,
904:  -5.27777777777777568E-002,
905:   4.72222222222222002E-002,
906:  -7.87037037037036924E-002,
907:  -4.72222222222222002E-002,
908:   1.38888888888888812E-002,
909:   1.20370370370370267E-002,
910:  -1.38888888888888794E-002,
911:  -4.72222222222222002E-002,
912:   4.90740740740740478E-002,
913:   1.38888888888888881E-002,
914:   1.38888888888888968E-002,
915:  -5.64814814814814659E-002,
916:  -1.38888888888888933E-002,
917:   5.27777777777777499E-002,
918:  -8.24074074074074153E-002,
919:   1.38888888888888812E-002,
920:  -1.38888888888888846E-002,
921:   4.90740740740740339E-002,
922:   4.72222222222222071E-002,
923:  -5.27777777777777568E-002,
924:   0.18981481481481477     ,
925:   5.27777777777777637E-002,
926:  -1.38888888888888829E-002,
927:  -5.27777777777777568E-002,
928:  -8.24074074074073598E-002,
929:   4.72222222222222071E-002,
930:  -4.72222222222222140E-002,
931:  -7.87037037037036924E-002,
932:   5.27777777777777637E-002,
933:   1.38888888888888916E-002,
934:  -8.24074074074073876E-002,
935:   1.38888888888888846E-002,
936:  -1.38888888888888951E-002,
937:  -5.64814814814814589E-002,
938:  -4.72222222222222071E-002,
939:   1.38888888888888812E-002,
940:   4.90740740740740339E-002,
941:   1.38888888888888829E-002,
942:  -1.38888888888888812E-002,
943:   1.20370370370370284E-002,
944:  -1.38888888888888881E-002,
945:   4.72222222222222071E-002,
946:   4.90740740740740339E-002,
947:  -5.27777777777777637E-002,
948:   5.27777777777777637E-002,
949:   0.18981481481481477     ,
950:   };

953:   PetscMemcpy(dd,DD,sizeof(PetscScalar)*576);
954:   return(0);
955: }