Actual source code: ex35.c

  1: /*T
  2:    Concepts: KSP^solving a system of linear equations
  3:    Concepts: KSP^Laplacian, 2d
  4:    Processors: n
  5: T*/

  7: /*


 10: Subduction Zone Benchmark

 12: */

 14: static char help[] = "Subduction Zone Benchmark\n\n";

 16:  #include petscmesh.h
 17:  #include petscksp.h
 18:  #include petscmg.h
 19:  #include petscdmmg.h

 21: PetscErrorCode MeshView_Sieve_Newer(ALE::Obj<ALE::Two::Mesh>, PetscViewer);
 22: PetscErrorCode CreateZoneBoundary(ALE::Obj<ALE::Two::Mesh>);

 24: PetscErrorCode updateOperator(Mat, ALE::Obj<ALE::Two::Mesh::field_type>, const ALE::Two::Mesh::point_type&, PetscScalar [], InsertMode);


 29: typedef enum {DIRICHLET, NEUMANN} BCType;

 31: typedef struct {
 32:   PetscScalar   nu;
 33:   BCType        bcType;
 34:   VecScatter  injection;
 35: } UserContext;

 37: PetscInt debug;

 41: int main(int argc,char **argv)
 42: {
 43:   MPI_Comm       comm;
 44:   DMMG          *dmmg;
 45:   UserContext    user;
 46:   PetscViewer    viewer;
 47:   const char    *bcTypes[2] = {"dirichlet", "neumann"};
 48:   PetscReal      refinementLimit, norm;
 49:   PetscInt       dim, bc, l;

 52:   PetscInitialize(&argc,&argv,(char *)0,help);
 53:   comm = PETSC_COMM_WORLD;

 55:   PetscOptionsBegin(comm, "", "Options for the inhomogeneous Poisson equation", "DMMG");
 56:     debug = 0;
 57:     PetscOptionsInt("-debug", "The debugging flag", "ex33.c", 0, &debug, PETSC_NULL);
 58:     dim  = 2;
 59:     PetscOptionsInt("-dim", "The mesh dimension", "ex33.c", 2, &dim, PETSC_NULL);
 60:     refinementLimit = 0.0;
 61:     PetscOptionsReal("-refinement_limit", "The area of the largest triangle in the mesh", "ex33.c", 1.0, &refinementLimit, PETSC_NULL);
 62:     user.nu = 0.1;
 63:     PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex33.c", 0.1, &user.nu, PETSC_NULL);
 64:     bc = (PetscInt)DIRICHLET;
 65:     PetscOptionsEList("-bc_type","Type of boundary condition","ex33.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
 66:     user.bcType = (BCType) bc;
 67:   PetscOptionsEnd();


 70:   ALE::Obj<ALE::Two::Mesh> meshBoundary = ALE::Two::Mesh(comm, 1, debug);
 71:   ALE::Obj<ALE::Two::Mesh> mesh;

 73:   try {
 74:     ALE::LogStage stage = ALE::LogStageRegister("MeshCreation");
 75:     ALE::LogStagePush(stage);
 76:     PetscPrintf(comm, "Generating mesh\n");
 77:     CreateZoneBoundary(meshBoundary);
 78:     mesh = ALE::Two::Generator::generate(meshBoundary);
 79:     ALE::Obj<ALE::Two::Mesh::sieve_type> topology = mesh->getTopology();
 80:     PetscPrintf(comm, "  Read %d elements\n", topology->heightStratum(0)->size());
 81:     PetscPrintf(comm, "  Read %d vertices\n", topology->depthStratum(0)->size());
 82:     ALE::LogStagePop(stage);

 84:     stage = ALE::LogStageRegister("BndCreation");
 85:     ALE::LogStagePush(stage);
 86:     PetscPrintf(comm, "Creating boundary\n");
 87:     ALE::Obj<ALE::Two::Mesh::field_type> boundary = mesh->getBoundary();
 88:     ALE::Obj<ALE::Two::Mesh::sieve_type::traits::depthSequence> bdVertices = topology->depthStratum(0, 1);
 89:     ALE::Two::Mesh::field_type::patch_type patch;

 91:     boundary->setTopology(topology);
 92:     boundary->setPatch(topology->leaves(), patch);
 93:     //boundary->setFiberDimensionByDepth(patch, 0, 1);
 94:     for(ALE::Two::Mesh::sieve_type::traits::depthSequence::iterator v_iter = bdVertices->begin(); v_iter != bdVertices->end(); ++v_iter) {
 95:       boundary->setFiberDimension(patch, *v_iter, 1);
 96:     }
 97:     boundary->orderPatches();
 98:     for(ALE::Two::Mesh::sieve_type::traits::depthSequence::iterator v_iter = bdVertices->begin(); v_iter != bdVertices->end(); ++v_iter) {
 99:       //double *coords = mesh->getCoordinates()->restrict(patch, *v_iter);
100:       double values[1] = {0.0};

102:       boundary->update(patch, *v_iter, values);
103:     }
104:     boundary->view("Mesh Boundary");
105:     ALE::LogStagePop(stage);

107:     stage = ALE::LogStageRegister("MeshDistribution");
108:     ALE::LogStagePush(stage);
109:     PetscPrintf(comm, "Distributing mesh\n");
110:     mesh = mesh->distribute();
111:     ALE::LogStagePop(stage);
112:     mesh->getBoundary()->view("Mesh Boundary");

114:     /* Refine... ideally we will refine based on the distance from the singularity.
115:        This "function call" based refinement has not yet been implemented.  Matt's doing it... */
116:     if (refinementLimit > 0.0) {
117:       stage = ALE::LogStageRegister("MeshRefine");
118:       ALE::LogStagePush(stage);
119:       PetscPrintf(comm, "Refining mesh\n");
120:       mesh = ALE::Two::Generator::refine(mesh, refinementLimit);
121:       ALE::LogStagePop(stage);
122:     }
123:     topology = mesh->getTopology();

125:     ALE::Obj<ALE::Two::Mesh::field_type> coords = mesh->getCoordinates();
126:     ALE::Obj<ALE::Two::Mesh::field_type> u = mesh->getField("u");
127:     u->setPatch(topology->leaves(), ALE::Two::Mesh::field_type::patch_type());
128:     u->setFiberDimensionByDepth(patch, 0, 1);
129:     u->orderPatches();
130:     u->createGlobalOrder();
131:     ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
132:     ALE::Obj<ALE::Two::Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
133:     std::string orderName("element");

135:     for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_iter = elements->begin(); e_iter != elements->end(); e_iter++) {
136:       // setFiberDimensionByDepth() does not work here since we only want it to apply to the patch cone
137:       //   What we really need is the depthStratum relative to the patch
138:       ALE::Obj<ALE::Two::Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_iter);

140:       coords->setPatch(orderName, cone, *e_iter);
141:       u->setPatch(orderName, cone, *e_iter);
142:       for(ALE::Two::Mesh::bundle_type::order_type::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
143:         u->setFiberDimension(orderName, *e_iter, *c_iter, 1);
144:       }
145:     }
146:     u->orderPatches(orderName);

148:     Mesh petscMesh;
149:     MeshCreate(comm, &petscMesh);
150:     MeshSetMesh(petscMesh, mesh);
151:     DMMGCreate(comm,3,PETSC_NULL,&dmmg);
152:     DMMGSetDM(dmmg, (DM) petscMesh);
153:     MeshDestroy(petscMesh);
154:     for (l = 0; l < DMMGGetLevels(dmmg); l++) {
155:       DMMGSetUser(dmmg,l,&user);
156:     }

158:     DMMGSetKSP(dmmg,ComputeRHS,ComputeMatrix);
159:     if (user.bcType == NEUMANN) {
160:       DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
161:     }
162:     MeshGetGlobalScatter(mesh, "u", DMMGGetx(dmmg), &user.injection);

164:     DMMGSolve(dmmg);

166:     MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
167:     VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
168:     VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
169:     PetscPrintf(comm,"Residual norm %g\n",norm);
170:     VecAssemblyBegin(DMMGGetx(dmmg));
171:     VecAssemblyEnd(DMMGGetx(dmmg));

173:     stage = ALE::LogStageRegister("MeshOutput");
174:     ALE::LogStagePush(stage);
175:     PetscPrintf(comm, "Creating VTK mesh file\n");
176:     PetscViewerCreate(comm, &viewer);
177:     PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
178:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
179:     PetscViewerFileSetName(viewer, "poisson.vtk");
180:     MeshView_Sieve_Newer(mesh, viewer);
181:     //VecView(DMMGGetRHS(dmmg), viewer);
182:     VecView(DMMGGetx(dmmg), viewer);
183:     PetscViewerDestroy(viewer);
184:     ALE::LogStagePop(stage);

186:     DMMGDestroy(dmmg);
187:   } catch (ALE::Exception e) {
188:     std::cout << e << std::endl;
189:   }
190:   PetscFinalize();
191:   return 0;
192: }

196: /*
197:   Subduction zone benchmark boundary

199:   Vertices       Edges

201:   0---------1     -----7-----
202:   |\        |    | 12        14
203:   | 2-------3    |  \---8----|
204:   |  \      |    11  \       |
205:   |   \     |    |    13     15
206:   |    \    |    |     \     |
207:   4-----5---6     --9-----10-


210: */

212: PetscErrorCode CreateZoneBoundary(ALE::Obj<ALE::Two::Mesh> mesh)
213: {
214:   MPI_Comm        comm = mesh->comm();
215:   ALE::Obj<ALE::Two::Mesh::sieve_type> topology = mesh->getTopology();
216:   PetscScalar   coords[14] = {          0.,  0.,
217:                                 660.,  0.,
218:                                  50., 50.,
219:                                 660., 50.,
220:                                   0.,600.,
221:                                 600.,600.,
222:                                 660.,600.};
223:   PetscInt        conns[18] = {        0, 1,
224:                                 1, 3,
225:                                 3, 6,
226:                                 6, 5,
227:                                 5, 4,
228:                                 4, 0,
229:                                 0, 2,
230:                                 2, 3,
231:                                 2, 5};


234:   PetscInt          order = 0;
235:   PetscMPIInt       rank;
236:   PetscErrorCode    ierr;

239:   MPI_Comm_rank(comm, &rank);
240:   mesh->populate(9, conns, 7, coords, false);

242:   /* Create boundary condition markers
243:         1 - surface
244:         2 - left inlet boundary
245:         3 - right corner flow boundary
246:         etc
247:   */
248:   if (rank == 0) {
249:       topology->setMarker(ALE::Two::Mesh::point_type(0, 0), 1);
250:       topology->setMarker(ALE::Two::Mesh::point_type(0, 1), 1);
251:       topology->setMarker(ALE::Two::Mesh::point_type(0, 3), 1);
252:       topology->setMarker(ALE::Two::Mesh::point_type(0, 4), 1);
253:       topology->setMarker(ALE::Two::Mesh::point_type(0, 5), 1);
254:       topology->setMarker(ALE::Two::Mesh::point_type(0, 6), 1);

256:       topology->setMarker(ALE::Two::Mesh::point_type(0, 7), 1);
257:       topology->setMarker(ALE::Two::Mesh::point_type(0, 8), 1);
258:       topology->setMarker(ALE::Two::Mesh::point_type(0, 9), 1);
259:       topology->setMarker(ALE::Two::Mesh::point_type(0, 10), 1);
260:       topology->setMarker(ALE::Two::Mesh::point_type(0, 11), 1);
261:       topology->setMarker(ALE::Two::Mesh::point_type(0, 12), 1);
262:   }
263:   return(0);
264: }


267: #ifndef MESH_3D

269: #define NUM_QUADRATURE_POINTS 9

271: /* Quadrature points */
272: static double points[18] = {
273:   -0.794564690381,
274:   -0.822824080975,
275:   -0.866891864322,
276:   -0.181066271119,
277:   -0.952137735426,
278:   0.575318923522,
279:   -0.0885879595127,
280:   -0.822824080975,
281:   -0.409466864441,
282:   -0.181066271119,
283:   -0.787659461761,
284:   0.575318923522,
285:   0.617388771355,
286:   -0.822824080975,
287:   0.0479581354402,
288:   -0.181066271119,
289:   -0.623181188096,
290:   0.575318923522};

292: /* Quadrature weights */
293: static double weights[9] = {
294:   0.223257681932,
295:   0.2547123404,
296:   0.0775855332238,
297:   0.357212291091,
298:   0.407539744639,
299:   0.124136853158,
300:   0.223257681932,
301:   0.2547123404,
302:   0.0775855332238};

304: #define NUM_BASIS_FUNCTIONS 3

306: /* Nodal basis function evaluations */
307: static double Basis[27] = {
308:   0.808694385678,
309:   0.10271765481,
310:   0.0885879595127,
311:   0.52397906772,
312:   0.0665540678392,
313:   0.409466864441,
314:   0.188409405952,
315:   0.0239311322871,
316:   0.787659461761,
317:   0.455706020244,
318:   0.455706020244,
319:   0.0885879595127,
320:   0.29526656778,
321:   0.29526656778,
322:   0.409466864441,
323:   0.10617026912,
324:   0.10617026912,
325:   0.787659461761,
326:   0.10271765481,
327:   0.808694385678,
328:   0.0885879595127,
329:   0.0665540678392,
330:   0.52397906772,
331:   0.409466864441,
332:   0.0239311322871,
333:   0.188409405952,
334:   0.787659461761};

336: /* Nodal basis function derivative evaluations */
337: static double BasisDerivatives[54] = {
338:   -0.5,
339:   -0.5,
340:   0.5,
341:   4.74937635818e-17,
342:   0.0,
343:   0.5,
344:   -0.5,
345:   -0.5,
346:   0.5,
347:   4.74937635818e-17,
348:   0.0,
349:   0.5,
350:   -0.5,
351:   -0.5,
352:   0.5,
353:   4.74937635818e-17,
354:   0.0,
355:   0.5,
356:   -0.5,
357:   -0.5,
358:   0.5,
359:   4.74937635818e-17,
360:   0.0,
361:   0.5,
362:   -0.5,
363:   -0.5,
364:   0.5,
365:   4.74937635818e-17,
366:   0.0,
367:   0.5,
368:   -0.5,
369:   -0.5,
370:   0.5,
371:   4.74937635818e-17,
372:   0.0,
373:   0.5,
374:   -0.5,
375:   -0.5,
376:   0.5,
377:   4.74937635818e-17,
378:   0.0,
379:   0.5,
380:   -0.5,
381:   -0.5,
382:   0.5,
383:   4.74937635818e-17,
384:   0.0,
385:   0.5,
386:   -0.5,
387:   -0.5,
388:   0.5,
389:   4.74937635818e-17,
390:   0.0,
391:   0.5};

393: #else

395: #define NUM_QUADRATURE_POINTS 27

397: /* Quadrature points */
398: static double points[81] = {
399:   -0.809560240317,
400:   -0.835756864273,
401:   -0.854011951854,
402:   -0.865851516496,
403:   -0.884304792128,
404:   -0.305992467923,
405:   -0.939397037651,
406:   -0.947733495427,
407:   0.410004419777,
408:   -0.876607962782,
409:   -0.240843539439,
410:   -0.854011951854,
411:   -0.913080888692,
412:   -0.465239359176,
413:   -0.305992467923,
414:   -0.960733394129,
415:   -0.758416359732,
416:   0.410004419777,
417:   -0.955631394718,
418:   0.460330056095,
419:   -0.854011951854,
420:   -0.968746121484,
421:   0.0286773243482,
422:   -0.305992467923,
423:   -0.985880737721,
424:   -0.53528439884,
425:   0.410004419777,
426:   -0.155115591937,
427:   -0.835756864273,
428:   -0.854011951854,
429:   -0.404851369974,
430:   -0.884304792128,
431:   -0.305992467923,
432:   -0.731135462175,
433:   -0.947733495427,
434:   0.410004419777,
435:   -0.452572254354,
436:   -0.240843539439,
437:   -0.854011951854,
438:   -0.61438408645,
439:   -0.465239359176,
440:   -0.305992467923,
441:   -0.825794030022,
442:   -0.758416359732,
443:   0.410004419777,
444:   -0.803159052121,
445:   0.460330056095,
446:   -0.854011951854,
447:   -0.861342428212,
448:   0.0286773243482,
449:   -0.305992467923,
450:   -0.937360010468,
451:   -0.53528439884,
452:   0.410004419777,
453:   0.499329056443,
454:   -0.835756864273,
455:   -0.854011951854,
456:   0.0561487765469,
457:   -0.884304792128,
458:   -0.305992467923,
459:   -0.522873886699,
460:   -0.947733495427,
461:   0.410004419777,
462:   -0.0285365459258,
463:   -0.240843539439,
464:   -0.854011951854,
465:   -0.315687284208,
466:   -0.465239359176,
467:   -0.305992467923,
468:   -0.690854665916,
469:   -0.758416359732,
470:   0.410004419777,
471:   -0.650686709523,
472:   0.460330056095,
473:   -0.854011951854,
474:   -0.753938734941,
475:   0.0286773243482,
476:   -0.305992467923,
477:   -0.888839283216,
478:   -0.53528439884,
479:   0.410004419777};

481: /* Quadrature weights */
482: static double weights[27] = {
483:   0.0701637994372,
484:   0.0653012061324,
485:   0.0133734490519,
486:   0.0800491405774,
487:   0.0745014590358,
488:   0.0152576273199,
489:   0.0243830167241,
490:   0.022693189565,
491:   0.0046474825267,
492:   0.1122620791,
493:   0.104481929812,
494:   0.021397518483,
495:   0.128078624924,
496:   0.119202334457,
497:   0.0244122037118,
498:   0.0390128267586,
499:   0.0363091033041,
500:   0.00743597204272,
501:   0.0701637994372,
502:   0.0653012061324,
503:   0.0133734490519,
504:   0.0800491405774,
505:   0.0745014590358,
506:   0.0152576273199,
507:   0.0243830167241,
508:   0.022693189565,
509:   0.0046474825267};

511: #define NUM_BASIS_FUNCTIONS 4

513: /* Nodal basis function evaluations */
514: static double Basis[108] = {
515:   0.749664528222,
516:   0.0952198798417,
517:   0.0821215678634,
518:   0.0729940240731,
519:   0.528074388273,
520:   0.0670742417521,
521:   0.0578476039361,
522:   0.347003766038,
523:   0.23856305665,
524:   0.0303014811743,
525:   0.0261332522867,
526:   0.705002209888,
527:   0.485731727037,
528:   0.0616960186091,
529:   0.379578230281,
530:   0.0729940240731,
531:   0.342156357896,
532:   0.0434595556538,
533:   0.267380320412,
534:   0.347003766038,
535:   0.154572667042,
536:   0.0196333029355,
537:   0.120791820134,
538:   0.705002209888,
539:   0.174656645238,
540:   0.0221843026408,
541:   0.730165028048,
542:   0.0729940240731,
543:   0.12303063253,
544:   0.0156269392579,
545:   0.514338662174,
546:   0.347003766038,
547:   0.0555803583921,
548:   0.00705963113955,
549:   0.23235780058,
550:   0.705002209888,
551:   0.422442204032,
552:   0.422442204032,
553:   0.0821215678634,
554:   0.0729940240731,
555:   0.297574315013,
556:   0.297574315013,
557:   0.0578476039361,
558:   0.347003766038,
559:   0.134432268912,
560:   0.134432268912,
561:   0.0261332522867,
562:   0.705002209888,
563:   0.273713872823,
564:   0.273713872823,
565:   0.379578230281,
566:   0.0729940240731,
567:   0.192807956775,
568:   0.192807956775,
569:   0.267380320412,
570:   0.347003766038,
571:   0.0871029849888,
572:   0.0871029849888,
573:   0.120791820134,
574:   0.705002209888,
575:   0.0984204739396,
576:   0.0984204739396,
577:   0.730165028048,
578:   0.0729940240731,
579:   0.0693287858938,
580:   0.0693287858938,
581:   0.514338662174,
582:   0.347003766038,
583:   0.0313199947658,
584:   0.0313199947658,
585:   0.23235780058,
586:   0.705002209888,
587:   0.0952198798417,
588:   0.749664528222,
589:   0.0821215678634,
590:   0.0729940240731,
591:   0.0670742417521,
592:   0.528074388273,
593:   0.0578476039361,
594:   0.347003766038,
595:   0.0303014811743,
596:   0.23856305665,
597:   0.0261332522867,
598:   0.705002209888,
599:   0.0616960186091,
600:   0.485731727037,
601:   0.379578230281,
602:   0.0729940240731,
603:   0.0434595556538,
604:   0.342156357896,
605:   0.267380320412,
606:   0.347003766038,
607:   0.0196333029355,
608:   0.154572667042,
609:   0.120791820134,
610:   0.705002209888,
611:   0.0221843026408,
612:   0.174656645238,
613:   0.730165028048,
614:   0.0729940240731,
615:   0.0156269392579,
616:   0.12303063253,
617:   0.514338662174,
618:   0.347003766038,
619:   0.00705963113955,
620:   0.0555803583921,
621:   0.23235780058,
622:   0.705002209888};

624: /* Nodal basis function derivative evaluations */
625: static double BasisDerivatives[324] = {
626:   -0.5,
627:   -0.5,
628:   -0.5,
629:   0.5,
630:   2.6964953901e-17,
631:   8.15881875835e-17,
632:   0.0,
633:   0.5,
634:   1.52600313755e-17,
635:   0.0,
636:   0.0,
637:   0.5,
638:   -0.5,
639:   -0.5,
640:   -0.5,
641:   0.5,
642:   2.6938349868e-17,
643:   1.08228622783e-16,
644:   0.0,
645:   0.5,
646:   -1.68114298577e-17,
647:   0.0,
648:   0.0,
649:   0.5,
650:   -0.5,
651:   -0.5,
652:   -0.5,
653:   0.5,
654:   2.69035912409e-17,
655:   1.43034809879e-16,
656:   0.0,
657:   0.5,
658:   -5.87133459397e-17,
659:   0.0,
660:   0.0,
661:   0.5,
662:   -0.5,
663:   -0.5,
664:   -0.5,
665:   0.5,
666:   2.6964953901e-17,
667:   2.8079494593e-17,
668:   0.0,
669:   0.5,
670:   1.52600313755e-17,
671:   0.0,
672:   0.0,
673:   0.5,
674:   -0.5,
675:   -0.5,
676:   -0.5,
677:   0.5,
678:   2.6938349868e-17,
679:   7.0536336094e-17,
680:   0.0,
681:   0.5,
682:   -1.68114298577e-17,
683:   0.0,
684:   0.0,
685:   0.5,
686:   -0.5,
687:   -0.5,
688:   -0.5,
689:   0.5,
690:   2.69035912409e-17,
691:   1.26006930238e-16,
692:   0.0,
693:   0.5,
694:   -5.87133459397e-17,
695:   0.0,
696:   0.0,
697:   0.5,
698:   -0.5,
699:   -0.5,
700:   -0.5,
701:   0.5,
702:   2.6964953901e-17,
703:   -3.49866380524e-17,
704:   0.0,
705:   0.5,
706:   1.52600313755e-17,
707:   0.0,
708:   0.0,
709:   0.5,
710:   -0.5,
711:   -0.5,
712:   -0.5,
713:   0.5,
714:   2.6938349868e-17,
715:   2.61116525673e-17,
716:   0.0,
717:   0.5,
718:   -1.68114298577e-17,
719:   0.0,
720:   0.0,
721:   0.5,
722:   -0.5,
723:   -0.5,
724:   -0.5,
725:   0.5,
726:   2.69035912409e-17,
727:   1.05937620823e-16,
728:   0.0,
729:   0.5,
730:   -5.87133459397e-17,
731:   0.0,
732:   0.0,
733:   0.5,
734:   -0.5,
735:   -0.5,
736:   -0.5,
737:   0.5,
738:   2.6964953901e-17,
739:   1.52807111565e-17,
740:   0.0,
741:   0.5,
742:   1.52600313755e-17,
743:   0.0,
744:   0.0,
745:   0.5,
746:   -0.5,
747:   -0.5,
748:   -0.5,
749:   0.5,
750:   2.6938349868e-17,
751:   6.1520690456e-17,
752:   0.0,
753:   0.5,
754:   -1.68114298577e-17,
755:   0.0,
756:   0.0,
757:   0.5,
758:   -0.5,
759:   -0.5,
760:   -0.5,
761:   0.5,
762:   2.69035912409e-17,
763:   1.21934019246e-16,
764:   0.0,
765:   0.5,
766:   -5.87133459397e-17,
767:   0.0,
768:   0.0,
769:   0.5,
770:   -0.5,
771:   -0.5,
772:   -0.5,
773:   0.5,
774:   2.6964953901e-17,
775:   -1.48832491782e-17,
776:   0.0,
777:   0.5,
778:   1.52600313755e-17,
779:   0.0,
780:   0.0,
781:   0.5,
782:   -0.5,
783:   -0.5,
784:   -0.5,
785:   0.5,
786:   2.6938349868e-17,
787:   4.0272766482e-17,
788:   0.0,
789:   0.5,
790:   -1.68114298577e-17,
791:   0.0,
792:   0.0,
793:   0.5,
794:   -0.5,
795:   -0.5,
796:   -0.5,
797:   0.5,
798:   2.69035912409e-17,
799:   1.1233505023e-16,
800:   0.0,
801:   0.5,
802:   -5.87133459397e-17,
803:   0.0,
804:   0.0,
805:   0.5,
806:   -0.5,
807:   -0.5,
808:   -0.5,
809:   0.5,
810:   2.6964953901e-17,
811:   -5.04349365259e-17,
812:   0.0,
813:   0.5,
814:   1.52600313755e-17,
815:   0.0,
816:   0.0,
817:   0.5,
818:   -0.5,
819:   -0.5,
820:   -0.5,
821:   0.5,
822:   2.6938349868e-17,
823:   1.52296507396e-17,
824:   0.0,
825:   0.5,
826:   -1.68114298577e-17,
827:   0.0,
828:   0.0,
829:   0.5,
830:   -0.5,
831:   -0.5,
832:   -0.5,
833:   0.5,
834:   2.69035912409e-17,
835:   1.01021564153e-16,
836:   0.0,
837:   0.5,
838:   -5.87133459397e-17,
839:   0.0,
840:   0.0,
841:   0.5,
842:   -0.5,
843:   -0.5,
844:   -0.5,
845:   0.5,
846:   2.6964953901e-17,
847:   -5.10267652705e-17,
848:   0.0,
849:   0.5,
850:   1.52600313755e-17,
851:   0.0,
852:   0.0,
853:   0.5,
854:   -0.5,
855:   -0.5,
856:   -0.5,
857:   0.5,
858:   2.6938349868e-17,
859:   1.4812758129e-17,
860:   0.0,
861:   0.5,
862:   -1.68114298577e-17,
863:   0.0,
864:   0.0,
865:   0.5,
866:   -0.5,
867:   -0.5,
868:   -0.5,
869:   0.5,
870:   2.69035912409e-17,
871:   1.00833228612e-16,
872:   0.0,
873:   0.5,
874:   -5.87133459397e-17,
875:   0.0,
876:   0.0,
877:   0.5,
878:   -0.5,
879:   -0.5,
880:   -0.5,
881:   0.5,
882:   2.6964953901e-17,
883:   -5.78459929494e-17,
884:   0.0,
885:   0.5,
886:   1.52600313755e-17,
887:   0.0,
888:   0.0,
889:   0.5,
890:   -0.5,
891:   -0.5,
892:   -0.5,
893:   0.5,
894:   2.6938349868e-17,
895:   1.00091968699e-17,
896:   0.0,
897:   0.5,
898:   -1.68114298577e-17,
899:   0.0,
900:   0.0,
901:   0.5,
902:   -0.5,
903:   -0.5,
904:   -0.5,
905:   0.5,
906:   2.69035912409e-17,
907:   9.86631702223e-17,
908:   0.0,
909:   0.5,
910:   -5.87133459397e-17,
911:   0.0,
912:   0.0,
913:   0.5,
914:   -0.5,
915:   -0.5,
916:   -0.5,
917:   0.5,
918:   2.6964953901e-17,
919:   -6.58832349994e-17,
920:   0.0,
921:   0.5,
922:   1.52600313755e-17,
923:   0.0,
924:   0.0,
925:   0.5,
926:   -0.5,
927:   -0.5,
928:   -0.5,
929:   0.5,
930:   2.6938349868e-17,
931:   4.34764891191e-18,
932:   0.0,
933:   0.5,
934:   -1.68114298577e-17,
935:   0.0,
936:   0.0,
937:   0.5,
938:   -0.5,
939:   -0.5,
940:   -0.5,
941:   0.5,
942:   2.69035912409e-17,
943:   9.61055074835e-17,
944:   0.0,
945:   0.5,
946:   -5.87133459397e-17,
947:   0.0,
948:   0.0,
949:   0.5};

951: #endif

955: PetscErrorCode ElementGeometry(ALE::Obj<ALE::Two::Mesh> mesh, const ALE::Two::Mesh::point_type& e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
956: {
957:   const double  *coords = mesh->getCoordinates()->restrict(std::string("element"), e);
958:   int            dim = mesh->getDimension();
959:   PetscReal      det, invDet;

963:   if (v0) {
964:     for(int d = 0; d < dim; d++) {
965:       v0[d] = coords[d];
966:     }
967:   }
968:   if (J) {
969:     for(int d = 0; d < dim; d++) {
970:       for(int f = 0; f < dim; f++) {
971:         J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
972:       }
973:     }
974:     if (debug) {
975:       MPI_Comm    comm = mesh->comm();
976:       PetscMPIInt rank;

978:       MPI_Comm_rank(comm, &rank);
979:       for(int d = 0; d < dim; d++) {
980:         if (d == 0) {
981:           PetscSynchronizedPrintf(comm, "[%d]J = /", rank);
982:         } else if (d == dim-1) {
983:           PetscSynchronizedPrintf(comm, "[%d]    \\", rank);
984:         } else {
985:           PetscSynchronizedPrintf(comm, "[%d]    |", rank);
986:         }
987:         for(int e = 0; e < dim; e++) {
988:           PetscSynchronizedPrintf(comm, " %g", J[d*dim+e]);
989:         }
990:         if (d == 0) {
991:           PetscSynchronizedPrintf(comm, " \\\n");
992:         } else if (d == dim-1) {
993:           PetscSynchronizedPrintf(comm, " /\n");
994:         } else {
995:           PetscSynchronizedPrintf(comm, " |\n");
996:         }
997:       }
998:     }
999:     if (dim == 2) {
1000:       det = J[0]*J[3] - J[1]*J[2];
1001:     } else if (dim == 3) {
1002:       det = J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
1003:             J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
1004:             J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1005:     }
1006:     invDet = 1.0/det;
1007:     if (detJ) {
1008:       if (det < 0) {SETERRQ(PETSC_ERR_ARG_WRONG, "Negative Matrix determinant");}
1009:       *detJ = det;
1010:     }
1011:     if (invJ) {
1012:       if (dim == 2) {
1013:         invJ[0] =  invDet*J[3];
1014:         invJ[1] = -invDet*J[1];
1015:         invJ[2] = -invDet*J[2];
1016:         invJ[3] =  invDet*J[0];
1017:       } else if (dim == 3) {
1018:         // FIX: This may be wrong
1019:         invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
1020:         invJ[0*3+1] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
1021:         invJ[0*3+2] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1022:         invJ[1*3+0] = invDet*(J[0*3+1]*J[2*3+2] - J[0*3+2]*J[2*3+1]);
1023:         invJ[1*3+1] = invDet*(J[0*3+2]*J[2*3+0] - J[0*3+0]*J[2*3+2]);
1024:         invJ[1*3+2] = invDet*(J[0*3+0]*J[2*3+1] - J[0*3+1]*J[2*3+0]);
1025:         invJ[2*3+0] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
1026:         invJ[2*3+1] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
1027:         invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
1028:       }
1029:       if (debug) {
1030:         MPI_Comm    comm = mesh->comm();
1031:         PetscMPIInt rank;

1033:         MPI_Comm_rank(comm, &rank);
1034:         for(int d = 0; d < dim; d++) {
1035:           if (d == 0) {
1036:             PetscSynchronizedPrintf(comm, "[%d]Jinv = /", rank);
1037:           } else if (d == dim-1) {
1038:             PetscSynchronizedPrintf(comm, "[%d]       \\", rank);
1039:           } else {
1040:             PetscSynchronizedPrintf(comm, "[%d]       |", rank);
1041:           }
1042:           for(int e = 0; e < dim; e++) {
1043:             PetscSynchronizedPrintf(comm, " %g", invJ[d*dim+e]);
1044:           }
1045:           if (d == 0) {
1046:             PetscSynchronizedPrintf(comm, " \\\n");
1047:           } else if (d == dim-1) {
1048:             PetscSynchronizedPrintf(comm, " /\n");
1049:           } else {
1050:             PetscSynchronizedPrintf(comm, " |\n");
1051:           }
1052:         }
1053:       }
1054:     }
1055:   }
1056:   return(0);
1057: }

1061: PetscErrorCode ComputeRho(PetscReal x, PetscReal y, PetscScalar *rho)
1062: {
1064:   if ((x > 1.0/3.0) && (x < 2.0/3.0) && (y > 1.0/3.0) && (y < 2.0/3.0)) {
1065:     //*rho = 100.0;
1066:     *rho = 1.0;
1067:   } else {
1068:     *rho = 1.0;
1069:   }
1070:   return(0);
1071: }

1075: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
1076: {
1077:   ALE::Obj<ALE::Two::Mesh> m;
1078:   Mesh                mesh = (Mesh) dmmg->dm;
1079:   UserContext        *user = (UserContext *) dmmg->user;
1080:   MPI_Comm            comm;
1081:   PetscReal           elementVec[NUM_BASIS_FUNCTIONS];
1082:   PetscReal           *v0, *Jac;
1083:   PetscReal           xi, eta, x_q, y_q, detJ, funcValue;
1084:   PetscInt            dim;
1085:   PetscInt            f, q;
1086:   PetscErrorCode      ierr;

1089:   PetscObjectGetComm((PetscObject) mesh, &comm);
1090:   MeshGetMesh(mesh, &m);
1091:   dim  = m->getDimension();
1092:   PetscMalloc(dim * sizeof(PetscReal), &v0);
1093:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1094:   ALE::Obj<ALE::Two::Mesh::field_type> field = m->getField("u");
1095:   ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
1096:   ALE::Two::Mesh::field_type::patch_type patch;
1097:   for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
1098:     ElementGeometry(m, *e_itor, v0, Jac, PETSC_NULL, &detJ);
1099:     /* Element integral */
1100:     PetscMemzero(elementVec, NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
1101:     for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1102:       xi = points[q*2+0] + 1.0;
1103:       eta = points[q*2+1] + 1.0;
1104:       x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1105:       y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1106:       funcValue = PetscExpScalar(-(x_q*x_q)/user->nu)*PetscExpScalar(-(y_q*y_q)/user->nu);
1107:       for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1108:         elementVec[f] += Basis[q*NUM_BASIS_FUNCTIONS+f]*funcValue*weights[q]*detJ;
1109:       }
1110:     }
1111:     if (debug) {PetscSynchronizedPrintf(comm, "elementVec = [%g %g %g]\n", elementVec[0], elementVec[1], elementVec[2]);}
1112:     /* Assembly */
1113:     field->updateAdd("element", *e_itor, elementVec);
1114:     if (debug) {PetscSynchronizedFlush(comm);}
1115:   }
1116:   PetscFree(v0);
1117:   PetscFree(Jac);

1119:   Vec locB;
1120:   VecCreateSeqWithArray(PETSC_COMM_SELF, field->getSize(patch), field->restrict(patch), &locB);
1121:   VecScatterBegin(user->injection, locB, b, ADD_VALUES, SCATTER_FORWARD);
1122:   VecScatterEnd(user->injection, locB, b, ADD_VALUES, SCATTER_FORWARD);
1123:   VecDestroy(locB);

1125:   /* force right hand side to be consistent for singular matrix */
1126:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
1127:   if (user->bcType == NEUMANN) {
1128:     MatNullSpace nullspace;

1130:     KSPGetNullSpace(dmmg->ksp,&nullspace);
1131:     MatNullSpaceRemove(nullspace,b,PETSC_NULL);
1132:   }
1133:   return(0);
1134: }

1138: PetscErrorCode ComputeMatrix(DMMG dmmg, Mat J, Mat jac)
1139: {
1140:   ALE::Obj<ALE::Two::Mesh> m;
1141:   Mesh              mesh = (Mesh) dmmg->dm;
1142:   UserContext      *user = (UserContext *) dmmg->user;
1143:   MPI_Comm          comm;
1144:   PetscReal         elementMat[NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS];
1145:   PetscReal        *v0, *Jac, *Jinv, *t_der, *b_der;
1146:   PetscReal         xi, eta, x_q, y_q, detJ, rho;
1147:   PetscInt          dim;
1148:   PetscInt          f, g, q;
1149:   PetscMPIInt       rank;
1150:   PetscErrorCode    ierr;

1153:   PetscObjectGetComm((PetscObject) mesh, &comm);
1154:   MPI_Comm_rank(comm, &rank);
1155:   MeshGetMesh(mesh, &m);
1156:   dim  = m->getDimension();
1157:   PetscMalloc(dim * sizeof(PetscReal), &v0);
1158:   PetscMalloc(dim * sizeof(PetscReal), &t_der);
1159:   PetscMalloc(dim * sizeof(PetscReal), &b_der);
1160:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1161:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jinv);
1162:   ALE::Obj<ALE::Two::Mesh::field_type> field = m->getField("u");
1163:   ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
1164:   for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
1165:    CHKMEMQ;
1166:     ElementGeometry(m, *e_itor, v0, Jac, Jinv, &detJ);
1167:     /* Element integral */
1168:     PetscMemzero(elementMat, NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
1169:     for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1170:       xi = points[q*2+0] + 1.0;
1171:       eta = points[q*2+1] + 1.0;
1172:       x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1173:       y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1174:       ComputeRho(x_q, y_q, &rho);
1175:       for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1176:         t_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
1177:         t_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
1178:         for(g = 0; g < NUM_BASIS_FUNCTIONS; g++) {
1179:           b_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
1180:           b_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
1181:           elementMat[f*NUM_BASIS_FUNCTIONS+g] += rho*(t_der[0]*b_der[0] + t_der[1]*b_der[1])*weights[q]*detJ;
1182:         }
1183:       }
1184:     }
1185:     if (debug) {
1186:       PetscSynchronizedPrintf(comm, "[%d]elementMat = [%g %g %g]\n                [%g %g %g]\n                [%g %g %g]\n",
1187:                                      rank, elementMat[0], elementMat[1], elementMat[2], elementMat[3], elementMat[4],
1188:                                      elementMat[5], elementMat[6], elementMat[7], elementMat[8]);
1189:     }
1190:     /* Assembly */
1191:     updateOperator(jac, field, *e_itor, elementMat, ADD_VALUES);
1192:     if (debug) {PetscSynchronizedFlush(comm);}
1193:   }
1194:   PetscFree(v0);
1195:   PetscFree(t_der);
1196:   PetscFree(b_der);
1197:   PetscFree(Jac);
1198:   PetscFree(Jinv);
1199:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
1200:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);

1202:   if (user->bcType == DIRICHLET) {
1203:     /* Zero out BC rows */
1204:     ALE::Obj<ALE::Two::Mesh::field_type> boundary = m->getBoundary();
1205:     ALE::Obj<ALE::Two::Mesh::sieve_type::traits::depthSequence> vertices = m->getTopology()->depthStratum(0);
1206:     ALE::Two::Mesh::field_type::patch_type patch;
1207:     PetscInt *boundaryIndices;
1208:     PetscInt  numBoundaryIndices = 0;
1209:     PetscInt  k = 0;

1211:     for(ALE::Two::Mesh::sieve_type::traits::depthSequence::iterator p = vertices->begin(); p != vertices->end(); ++p) {
1212:       if (boundary->getIndex(patch, *p).index > 0) {
1213:         const ALE::Two::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);

1215:         if (idx.index > 0) {
1216:           numBoundaryIndices += idx.index;
1217:         }
1218:       }
1219:     }
1220:     PetscMalloc(numBoundaryIndices * sizeof(PetscInt), &boundaryIndices);
1221:     for(ALE::Two::Mesh::sieve_type::traits::depthSequence::iterator p = vertices->begin(); p != vertices->end(); ++p) {
1222:       if (boundary->getIndex(patch, *p).index > 0) {
1223:         const ALE::Two::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);

1225:         for(int i = 0; i < idx.index; i++) {
1226:           boundaryIndices[k++] = idx.prefix + i;
1227:         }
1228:       }
1229:     }
1230:     if (debug) {
1231:       for(int i = 0; i < numBoundaryIndices; i++) {
1232:         PetscSynchronizedPrintf(comm, "[%d]boundaryIndices[%d] = %d\n", rank, i, boundaryIndices[i]);
1233:       }
1234:     }
1235:     PetscSynchronizedFlush(comm);
1236:     MatZeroRows(jac, numBoundaryIndices, boundaryIndices, 1.0);
1237:     PetscFree(boundaryIndices);
1238:   }
1239:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
1240:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
1241:   return(0);
1242: }