Actual source code: ex33.c

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

  7: /*
  8: Laplacian in 2D. Modeled by the partial differential equation

 10:    -div grad u = f,  0 < x,y < 2,

 12: with forcing function

 14:    f = 2
 15:    f = 2y (1 - y) + 2x (1 - x)

 17: with Dirichlet boundary conditions

 19:    u = 0 for x = 0, x = 2 and u = x (2 - x) for y = 0, y = 2
 20:    u = 0 for x = 0, x = 2, y = 0, y = 2

 22: or pure Neumman boundary conditions.

 24: This has exact solution

 26:    u = x (2 - x)
 27:    u = x (2 - x) y (2 - y)

 29: This uses multigrid to solve the linear system

 31: The 2D test mesh

 33:          13
 34:   14--29----31---12
 35:     |\    |\    |
 36:     2 2 5 2 3 7 3
 37:     7  6  8  0  2
 38:     | 4 \ | 6 \ |
 39:     |    \|    \|
 40:   15--20-16-24---11
 41:     |\    |\    |
 42:     1 1 1 2 2 3 2
 43:     8  7  1  2  5
 44:     | 0 \ | 2 \ |
 45:     |    \|    \|
 46:    8--19----23---10
 47:           9

 49: */

 51: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";

 53:  #include petscmesh.h
 54:  #include petscksp.h
 55:  #include petscmg.h
 56:  #include petscdmmg.h

 58: PetscErrorCode MeshView_Sieve_Newer(ALE::Obj<ALE::Two::Mesh>, PetscViewer);
 59: PetscErrorCode CreateMeshBoundary(ALE::Obj<ALE::Two::Mesh>);
 60: PetscErrorCode updateOperator(Mat, ALE::Obj<ALE::Two::Mesh::field_type>, const ALE::Two::Mesh::point_type&, PetscScalar [], InsertMode);


 67: typedef enum {DIRICHLET, NEUMANN} BCType;

 69: typedef struct {
 70:   BCType      bcType;
 71:   VecScatter  injection;
 72: } UserContext;

 74: PetscInt debug;

 78: int main(int argc,char **argv)
 79: {
 80:   MPI_Comm       comm;
 81:   DMMG          *dmmg;
 82:   UserContext    user;
 83:   PetscViewer    viewer;
 84:   const char    *bcTypes[2] = {"dirichlet", "neumann"};
 85:   PetscReal      refinementLimit, norm, error;
 86:   PetscInt       dim, bc, l, meshDebug;

 89:   PetscInitialize(&argc,&argv,(char *)0,help);
 90:   comm = PETSC_COMM_WORLD;
 91:   PetscOptionsBegin(comm, "", "Options for the inhomogeneous Poisson equation", "DMMG");
 92:     debug = 0;
 93:     PetscOptionsInt("-debug", "The debugging flag", "ex33.c", 0, &debug, PETSC_NULL);
 94:     meshDebug = 0;
 95:     PetscOptionsInt("-mesh_debug", "The mesh debugging flag", "ex33.c", 0, &meshDebug, PETSC_NULL);
 96:     dim  = 2;
 97:     PetscOptionsInt("-dim", "The mesh dimension", "ex33.c", 2, &dim, PETSC_NULL);
 98:     refinementLimit = 0.0;
 99:     PetscOptionsReal("-refinement_limit", "The area of the largest triangle in the mesh", "ex33.c", 1.0, &refinementLimit, PETSC_NULL);
100:     bc = (PetscInt)DIRICHLET;
101:     PetscOptionsEList("-bc_type","Type of boundary condition","ex33.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
102:     user.bcType = (BCType) bc;
103:   PetscOptionsEnd();

105:   ALE::Obj<ALE::Two::Mesh> meshBoundary = ALE::Two::Mesh(comm, dim-1, meshDebug);
106:   ALE::Obj<ALE::Two::Mesh> mesh;

108:   try {
109:     ALE::LogStage stage = ALE::LogStageRegister("MeshCreation");
110:     ALE::LogStagePush(stage);
111:     PetscPrintf(comm, "Generating mesh\n");
112:     CreateMeshBoundary(meshBoundary);
113:     mesh = ALE::Two::Generator::generate(meshBoundary);
114:     ALE::Obj<ALE::Two::Mesh::sieve_type> topology = mesh->getTopology();
115:     PetscPrintf(comm, "  Generated %d elements\n", topology->heightStratum(0)->size());
116:     PetscPrintf(comm, "  Generated %d vertices\n", topology->depthStratum(0)->size());
117:     ALE::LogStagePop(stage);

119:     stage = ALE::LogStageRegister("MeshDistribution");
120:     ALE::LogStagePush(stage);
121:     PetscPrintf(comm, "Distributing mesh\n");
122:     mesh = mesh->distribute();
123:     ALE::LogStagePop(stage);

125:     if (refinementLimit > 0.0) {
126:       stage = ALE::LogStageRegister("MeshRefine");
127:       ALE::LogStagePush(stage);
128:       PetscPrintf(comm, "Refining mesh\n");
129:       mesh = ALE::Two::Generator::refine(mesh, refinementLimit);
130:       ALE::LogStagePop(stage);
131:     }
132:     topology = mesh->getTopology();

134:     stage = ALE::LogStageRegister("BndValues");
135:     ALE::LogStagePush(stage);
136:     PetscPrintf(comm, "Calculating boundary values\n");
137:     ALE::Obj<ALE::Two::Mesh::field_type> boundary = mesh->getBoundary();
138:     ALE::Obj<ALE::Two::Mesh::sieve_type::traits::depthSequence> vertices = topology->depthStratum(0);
139:     ALE::Two::Mesh::field_type::patch_type bdPatch(0, 1);
140:     ALE::Two::Mesh::field_type::patch_type patch;

142:     for(ALE::Two::Mesh::sieve_type::traits::depthSequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
143:       if (boundary->getIndex(bdPatch, *v_iter).index > 0) {
144:         const double *coords = mesh->getCoordinates()->restrict(patch, *v_iter);
145:         double        values[1];

147:         values[0] = coords[0]*(2.0 - coords[0]);
148:         boundary->update(bdPatch, *v_iter, values);
149:       }
150:     }
151:     if (debug) {boundary->view("Mesh Boundary");}
152:     ALE::LogStagePop(stage);

154:     ALE::Obj<ALE::Two::Mesh::field_type> u = mesh->getField("u");
155:     ALE::Obj<ALE::Two::Mesh::field_type> b = mesh->getField("b");
156:     u->setPatch(topology->leaves(), ALE::Two::Mesh::field_type::patch_type());
157:     u->setFiberDimensionByDepth(patch, 0, 1);
158:     u->orderPatches();
159:     if (debug) {u->view("u");}
160:     u->createGlobalOrder();
161:     b->setPatch(topology->leaves(), ALE::Two::Mesh::field_type::patch_type());
162:     b->setFiberDimensionByDepth(patch, 0, 1);
163:     b->orderPatches();
164:     if (debug) {b->view("b");}
165:     b->createGlobalOrder();
166:     ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
167:     ALE::Obj<ALE::Two::Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
168:     std::string orderName("element");

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

175:       u->setPatch(orderName, cone, *e_iter);
176:       b->setPatch(orderName, cone, *e_iter);
177:       for(ALE::Two::Mesh::bundle_type::order_type::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
178:         u->setFiberDimension(orderName, *e_iter, *c_iter, 1);
179:         b->setFiberDimension(orderName, *e_iter, *c_iter, 1);
180:       }
181:     }
182:     u->orderPatches(orderName);
183:     b->orderPatches(orderName);
184:     CheckElementGeometry(mesh);

186:     Mesh petscMesh;
187:     MeshCreate(comm, &petscMesh);
188:     MeshSetMesh(petscMesh, mesh);
189:     DMMGCreate(comm,3,PETSC_NULL,&dmmg);
190:     DMMGSetDM(dmmg, (DM) petscMesh);
191:     MeshDestroy(petscMesh);
192:     for (l = 0; l < DMMGGetLevels(dmmg); l++) {
193:       DMMGSetUser(dmmg,l,&user);
194:     }

196:     DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);
197:     if (user.bcType == NEUMANN) {
198:       DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
199:     }
200:     MeshGetGlobalScatter(mesh, "u", DMMGGetx(dmmg), &user.injection);

202:     DMMGSolve(dmmg);

204:     if (debug) {
205:       PetscPrintf(mesh->comm(), "Solution vector:");
206:       VecView(DMMGGetx(dmmg), PETSC_VIEWER_STDOUT_WORLD);
207:     }
208:     ComputeError(dmmg[DMMGGetLevels(dmmg)-1], DMMGGetx(dmmg), &error);
209:     PetscPrintf(comm,"Error norm %g\n",error);

211:     MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
212:     VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
213:     VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
214:     PetscPrintf(comm,"Residual norm %g\n",norm);
215:     VecAssemblyBegin(DMMGGetx(dmmg));
216:     VecAssemblyEnd(DMMGGetx(dmmg));

218:     stage = ALE::LogStageRegister("MeshOutput");
219:     ALE::LogStagePush(stage);
220:     PetscPrintf(comm, "Creating VTK mesh file\n");
221:     PetscViewerCreate(comm, &viewer);
222:     PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
223:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
224:     PetscViewerFileSetName(viewer, "poisson.vtk");
225:     MeshView_Sieve_Newer(mesh, viewer);
226:     //VecView(DMMGGetRHS(dmmg), viewer);
227:     VecView(DMMGGetx(dmmg), viewer);
228:     PetscViewerDestroy(viewer);
229:     ALE::LogStagePop(stage);

231:     DMMGDestroy(dmmg);
232:   } catch (ALE::Exception e) {
233:     std::cout << e << std::endl;
234:   }
235:   PetscFinalize();
236:   return 0;
237: }

241: /*
242:   Simple square boundary:

244:   6--14-5--13-4
245:   |     |     |
246:   15   19    12
247:   |     |     |
248:   7--20-8--18-3
249:   |     |     |
250:   16   17    11
251:   |     |     |
252:   0--9--1--10-2
253: */
254: PetscErrorCode CreateSquareBoundary(ALE::Obj<ALE::Two::Mesh> mesh)
255: {
256:   MPI_Comm          comm = mesh->comm();
257:   ALE::Obj<ALE::Two::Mesh::sieve_type> topology = mesh->getTopology();
258:   PetscScalar       coords[18] = {0.0, 0.0,
259:                                   1.0, 0.0,
260:                                   2.0, 0.0,
261:                                   2.0, 1.0,
262:                                   2.0, 2.0,
263:                                   1.0, 2.0,
264:                                   0.0, 2.0,
265:                                   0.0, 1.0,
266:                                   1.0, 1.0};
267:   PetscInt    connectivity[40] = {0, 1,
268:                                   1, 2,
269:                                   2, 3,
270:                                   3, 4,
271:                                   4, 5,
272:                                   5, 6,
273:                                   6, 7,
274:                                   7, 0,
275:                                   1, 8,
276:                                   3, 8,
277:                                   5, 8,
278:                                   7, 8};
279:   ALE::Two::Mesh::point_type vertices[9];
280:   PetscInt          order = 0;
281:   PetscMPIInt       rank;
282:   PetscErrorCode    ierr;

285:   MPI_Comm_rank(comm, &rank);
286:   if (rank == 0) {
287:     ALE::Two::Mesh::point_type edge;

289:     /* Create topology and ordering */
290:     for(int v = 0; v < 9; v++) {
291:       vertices[v] = ALE::Two::Mesh::point_type(0, v);
292:     }
293:     for(int e = 9; e < 17; e++) {
294:       edge = ALE::Two::Mesh::point_type(0, e);
295:       topology->addArrow(vertices[e-9],     edge, order++);
296:       topology->addArrow(vertices[(e-8)%8], edge, order++);
297:     }
298:     edge = ALE::Two::Mesh::point_type(0, 17);
299:     topology->addArrow(vertices[1], edge, order++);
300:     topology->addArrow(vertices[8], edge, order++);
301:     edge = ALE::Two::Mesh::point_type(0, 18);
302:     topology->addArrow(vertices[3], edge, order++);
303:     topology->addArrow(vertices[8], edge, order++);
304:     edge = ALE::Two::Mesh::point_type(0, 19);
305:     topology->addArrow(vertices[5], edge, order++);
306:     topology->addArrow(vertices[8], edge, order++);
307:     edge = ALE::Two::Mesh::point_type(0, 20);
308:     topology->addArrow(vertices[7], edge, order++);
309:     topology->addArrow(vertices[8], edge, order++);
310:   }
311:   topology->stratify();
312:   mesh->createVertexBundle(20, connectivity);
313:   mesh->createSerialCoordinates(2, 0, coords);
314:   /* Create boundary conditions */
315:   if (rank == 0) {
316:     for(int v = 0; v < 8; v++) {
317:       topology->setMarker(ALE::Two::Mesh::point_type(0, v), 1);
318:     }
319:     for(int e = 9; e < 17; e++) {
320:       topology->setMarker(ALE::Two::Mesh::point_type(0, e), 1);
321:     }
322:   }
323:   return(0);
324: }

328: /*
329:   Simple cube boundary:

331:       7-----6
332:      /|    /|
333:     3-----2 |
334:     | |   | |
335:     | 4---|-5
336:     |/    |/
337:     0-----1
338: */
339: PetscErrorCode CreateCubeBoundary(ALE::Obj<ALE::Two::Mesh> mesh)
340: {
341:   MPI_Comm          comm = mesh->comm();
342:   ALE::Obj<ALE::Two::Mesh::sieve_type> topology = mesh->getTopology();
343:   PetscScalar       coords[24] = {0.0, 0.0, 0.0,
344:                                   1.0, 0.0, 0.0,
345:                                   1.0, 1.0, 0.0,
346:                                   0.0, 1.0, 0.0,
347:                                   0.0, 0.0, 1.0,
348:                                   1.0, 0.0, 1.0,
349:                                   1.0, 1.0, 1.0,
350:                                   0.0, 1.0, 1.0};
351:   PetscInt    connectivity[24] = {0, 1, 2, 3,
352:                                   7, 6, 5, 4,
353:                                   0, 4, 5, 1,
354:                                   1, 5, 6, 2,
355:                                   2, 6, 7, 3,
356:                                   3, 7, 4, 0};
357:   ALE::Obj<std::set<ALE::Two::Mesh::point_type> > cone = std::set<ALE::Two::Mesh::point_type>();
358:   ALE::Two::Mesh::point_type            vertices[8];
359:   ALE::Two::Mesh::point_type            edges[12];
360:   ALE::Two::Mesh::point_type            edge;
361:   PetscInt                              embedDim = 3;
362:   PetscInt                              order = 0;
363:   PetscMPIInt                           rank;
364:   PetscErrorCode                        ierr;

367:   MPI_Comm_rank(comm, &rank);
368:   if (rank == 0) {
369:     ALE::Two::Mesh::point_type face;

371:     /* Create topology and ordering */
372:     /* Vertices: 0 .. 3 on the bottom of the cube, 4 .. 7 on the top */
373:     for(int v = 0; v < 8; v++) {
374:       vertices[v] = ALE::Two::Mesh::point_type(0, v);
375:     }

377:     /* Edges on the bottom: Sieve element numbers e = 8 .. 11, edge numbers e - 8 = 0 .. 3 */
378:     for(int e = 8; e < 12; e++) {
379:       edge = ALE::Two::Mesh::point_type(0, e);
380:       edges[e-8] = edge;
381:       topology->addArrow(vertices[e-8],     edge, order++);
382:       topology->addArrow(vertices[(e-7)%4], edge, order++);
383:     }
384:     /* Edges on the top: Sieve element numbers e = 12 .. 15, edge numbers e - 8 = 4 .. 7 */
385:     for(int e = 12; e < 16; e++) {
386:       edge = ALE::Two::Mesh::point_type(0, e);
387:       edges[e-8] = edge;
388:       topology->addArrow(vertices[e-8],        edge, order++);
389:       topology->addArrow(vertices[(e-11)%4+4], edge, order++);
390:     }
391:     /* Edges from bottom to top: Sieve element numbers e = 16 .. 19, edge numbers e - 8 = 8 .. 11 */
392:     for(int e = 16; e < 20; e++) {
393:       edge = ALE::Two::Mesh::point_type(0, e);
394:       edges[e-8] = edge;
395:       topology->addArrow(vertices[e-16],   edge, order++);
396:       topology->addArrow(vertices[e-16+4], edge, order++);
397:     }

399:     /* Bottom face */
400:     face = ALE::Two::Mesh::point_type(0, 20);
401:     topology->addArrow(edges[0], face, order++);
402:     topology->addArrow(edges[1], face, order++);
403:     topology->addArrow(edges[2], face, order++);
404:     topology->addArrow(edges[3], face, order++);
405:     /* Top face */
406:     face = ALE::Two::Mesh::point_type(0, 21);
407:     topology->addArrow(edges[4], face, order++);
408:     topology->addArrow(edges[5], face, order++);
409:     topology->addArrow(edges[6], face, order++);
410:     topology->addArrow(edges[7], face, order++);
411:     /* Side faces: f = 22 .. 25 */
412:     for(int f = 22; f < 26; f++) {
413:       face = ALE::Two::Mesh::point_type(0, f);
414:       int v = f - 22;
415:       /* Covered by edges f - 22, f - 22 + 4, f - 22 + 8, (f - 21)%4 + 8 */
416:       topology->addArrow(edges[v],         face, order++);
417:       topology->addArrow(edges[(v+1)%4+8], face, order++);
418:       topology->addArrow(edges[v+4],       face, order++);
419:       topology->addArrow(edges[v+8],       face, order++);
420:     }
421:   }/* if(rank == 0) */
422:   topology->stratify();
423:   if (rank == 0) {
424:     ALE::Obj<ALE::Two::Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
425:     ALE::Obj<ALE::Two::Mesh::bundle_type::PointArray> points = ALE::Two::Mesh::bundle_type::PointArray();
426:     const std::string orderName("element");
427:     /* Bottom face */
428:     ALE::Two::Mesh::point_type face = ALE::Two::Mesh::point_type(0, 20);
429:     points->clear();
430:     points->push_back(vertices[0]);
431:     points->push_back(vertices[1]);
432:     points->push_back(vertices[2]);
433:     points->push_back(vertices[3]);
434:     vertexBundle->setPatch(orderName, points, face);
435:     /* Top face */
436:     face = ALE::Two::Mesh::point_type(0, 21);
437:     points->clear();
438:     points->push_back(vertices[4]);
439:     points->push_back(vertices[5]);
440:     points->push_back(vertices[6]);
441:     points->push_back(vertices[7]);
442:     vertexBundle->setPatch(orderName, points, face);
443:     /* Side faces: f = 22 .. 25 */
444:     for(int f = 22; f < 26; f++) {
445:       face = ALE::Two::Mesh::point_type(0, f);
446:       int v = f - 22;
447:       /* Covered by edges f - 22, f - 22 + 4, f - 22 + 8, (f - 21)%4 + 8 */
448:       points->clear();
449:       points->push_back(vertices[v]);
450:       points->push_back(vertices[(v+1)%4]);
451:       points->push_back(vertices[(v+1)%4+4]);
452:       points->push_back(vertices[v+4]);
453:       vertexBundle->setPatch(orderName, points, face);
454:     }
455:   }
456:   mesh->createVertexBundle(6, connectivity);
457:   mesh->createSerialCoordinates(embedDim, 0, coords);

459:   /* Create boundary conditions: set marker 1 to all of the sieve elements, 
460:      since everything is on the boundary (no internal faces, edges or vertices)  */
461:   if (rank == 0) {
462:     /* set marker to the base of the topology sieve -- the faces and the edges */
463:     topology->setMarker(topology->base(), 1);
464:     /* set marker to the vertices -- the 0-depth stratum */
465:     topology->setMarker(topology->depthStratum(0), 1);
466:   }

468:   return(0);
469: }

473: /*
474:   Simple square boundary:

476:   6--14-5--13-4
477:   |     |     |
478:   15   19    12
479:   |     |     |
480:   7--20-8--18-3
481:   |     |     |
482:   16   17    11
483:   |     |     |
484:   0--9--1--10-2
485: */
486: PetscErrorCode CreateMeshBoundary(ALE::Obj<ALE::Two::Mesh> mesh)
487: {
488:   int            dim = mesh->getDimension();

492:   if (dim == 1) {
493:     CreateSquareBoundary(mesh);
494:   } else if (dim == 2) {
495:     CreateCubeBoundary(mesh);
496:   } else {
497:     SETERRQ1(PETSC_ERR_SUP, "Cannot construct a boundary of dimension %d", dim);
498:   }
499:   return(0);
500: }

502: #ifndef MESH_3D

504: #define NUM_QUADRATURE_POINTS 9

506: /* Quadrature points */
507: static double points[18] = {
508:   -0.794564690381,
509:   -0.822824080975,
510:   -0.866891864322,
511:   -0.181066271119,
512:   -0.952137735426,
513:   0.575318923522,
514:   -0.0885879595127,
515:   -0.822824080975,
516:   -0.409466864441,
517:   -0.181066271119,
518:   -0.787659461761,
519:   0.575318923522,
520:   0.617388771355,
521:   -0.822824080975,
522:   0.0479581354402,
523:   -0.181066271119,
524:   -0.623181188096,
525:   0.575318923522};

527: /* Quadrature weights */
528: static double weights[9] = {
529:   0.223257681932,
530:   0.2547123404,
531:   0.0775855332238,
532:   0.357212291091,
533:   0.407539744639,
534:   0.124136853158,
535:   0.223257681932,
536:   0.2547123404,
537:   0.0775855332238};

539: #define NUM_BASIS_FUNCTIONS 3

541: /* Nodal basis function evaluations */
542: static double Basis[27] = {
543:   0.808694385678,
544:   0.10271765481,
545:   0.0885879595127,
546:   0.52397906772,
547:   0.0665540678392,
548:   0.409466864441,
549:   0.188409405952,
550:   0.0239311322871,
551:   0.787659461761,
552:   0.455706020244,
553:   0.455706020244,
554:   0.0885879595127,
555:   0.29526656778,
556:   0.29526656778,
557:   0.409466864441,
558:   0.10617026912,
559:   0.10617026912,
560:   0.787659461761,
561:   0.10271765481,
562:   0.808694385678,
563:   0.0885879595127,
564:   0.0665540678392,
565:   0.52397906772,
566:   0.409466864441,
567:   0.0239311322871,
568:   0.188409405952,
569:   0.787659461761};

571: /* Nodal basis function derivative evaluations */
572: static double BasisDerivatives[54] = {
573:   -0.5,
574:   -0.5,
575:   0.5,
576:   4.74937635818e-17,
577:   0.0,
578:   0.5,
579:   -0.5,
580:   -0.5,
581:   0.5,
582:   4.74937635818e-17,
583:   0.0,
584:   0.5,
585:   -0.5,
586:   -0.5,
587:   0.5,
588:   4.74937635818e-17,
589:   0.0,
590:   0.5,
591:   -0.5,
592:   -0.5,
593:   0.5,
594:   4.74937635818e-17,
595:   0.0,
596:   0.5,
597:   -0.5,
598:   -0.5,
599:   0.5,
600:   4.74937635818e-17,
601:   0.0,
602:   0.5,
603:   -0.5,
604:   -0.5,
605:   0.5,
606:   4.74937635818e-17,
607:   0.0,
608:   0.5,
609:   -0.5,
610:   -0.5,
611:   0.5,
612:   4.74937635818e-17,
613:   0.0,
614:   0.5,
615:   -0.5,
616:   -0.5,
617:   0.5,
618:   4.74937635818e-17,
619:   0.0,
620:   0.5,
621:   -0.5,
622:   -0.5,
623:   0.5,
624:   4.74937635818e-17,
625:   0.0,
626:   0.5};

628: #else

630: #define NUM_QUADRATURE_POINTS 27

632: /* Quadrature points */
633: static double points[81] = {
634:   -0.809560240317,
635:   -0.835756864273,
636:   -0.854011951854,
637:   -0.865851516496,
638:   -0.884304792128,
639:   -0.305992467923,
640:   -0.939397037651,
641:   -0.947733495427,
642:   0.410004419777,
643:   -0.876607962782,
644:   -0.240843539439,
645:   -0.854011951854,
646:   -0.913080888692,
647:   -0.465239359176,
648:   -0.305992467923,
649:   -0.960733394129,
650:   -0.758416359732,
651:   0.410004419777,
652:   -0.955631394718,
653:   0.460330056095,
654:   -0.854011951854,
655:   -0.968746121484,
656:   0.0286773243482,
657:   -0.305992467923,
658:   -0.985880737721,
659:   -0.53528439884,
660:   0.410004419777,
661:   -0.155115591937,
662:   -0.835756864273,
663:   -0.854011951854,
664:   -0.404851369974,
665:   -0.884304792128,
666:   -0.305992467923,
667:   -0.731135462175,
668:   -0.947733495427,
669:   0.410004419777,
670:   -0.452572254354,
671:   -0.240843539439,
672:   -0.854011951854,
673:   -0.61438408645,
674:   -0.465239359176,
675:   -0.305992467923,
676:   -0.825794030022,
677:   -0.758416359732,
678:   0.410004419777,
679:   -0.803159052121,
680:   0.460330056095,
681:   -0.854011951854,
682:   -0.861342428212,
683:   0.0286773243482,
684:   -0.305992467923,
685:   -0.937360010468,
686:   -0.53528439884,
687:   0.410004419777,
688:   0.499329056443,
689:   -0.835756864273,
690:   -0.854011951854,
691:   0.0561487765469,
692:   -0.884304792128,
693:   -0.305992467923,
694:   -0.522873886699,
695:   -0.947733495427,
696:   0.410004419777,
697:   -0.0285365459258,
698:   -0.240843539439,
699:   -0.854011951854,
700:   -0.315687284208,
701:   -0.465239359176,
702:   -0.305992467923,
703:   -0.690854665916,
704:   -0.758416359732,
705:   0.410004419777,
706:   -0.650686709523,
707:   0.460330056095,
708:   -0.854011951854,
709:   -0.753938734941,
710:   0.0286773243482,
711:   -0.305992467923,
712:   -0.888839283216,
713:   -0.53528439884,
714:   0.410004419777};

716: /* Quadrature weights */
717: static double weights[27] = {
718:   0.0701637994372,
719:   0.0653012061324,
720:   0.0133734490519,
721:   0.0800491405774,
722:   0.0745014590358,
723:   0.0152576273199,
724:   0.0243830167241,
725:   0.022693189565,
726:   0.0046474825267,
727:   0.1122620791,
728:   0.104481929812,
729:   0.021397518483,
730:   0.128078624924,
731:   0.119202334457,
732:   0.0244122037118,
733:   0.0390128267586,
734:   0.0363091033041,
735:   0.00743597204272,
736:   0.0701637994372,
737:   0.0653012061324,
738:   0.0133734490519,
739:   0.0800491405774,
740:   0.0745014590358,
741:   0.0152576273199,
742:   0.0243830167241,
743:   0.022693189565,
744:   0.0046474825267};

746: #define NUM_BASIS_FUNCTIONS 4

748: /* Nodal basis function evaluations */
749: static double Basis[108] = {
750:   0.749664528222,
751:   0.0952198798417,
752:   0.0821215678634,
753:   0.0729940240731,
754:   0.528074388273,
755:   0.0670742417521,
756:   0.0578476039361,
757:   0.347003766038,
758:   0.23856305665,
759:   0.0303014811743,
760:   0.0261332522867,
761:   0.705002209888,
762:   0.485731727037,
763:   0.0616960186091,
764:   0.379578230281,
765:   0.0729940240731,
766:   0.342156357896,
767:   0.0434595556538,
768:   0.267380320412,
769:   0.347003766038,
770:   0.154572667042,
771:   0.0196333029355,
772:   0.120791820134,
773:   0.705002209888,
774:   0.174656645238,
775:   0.0221843026408,
776:   0.730165028048,
777:   0.0729940240731,
778:   0.12303063253,
779:   0.0156269392579,
780:   0.514338662174,
781:   0.347003766038,
782:   0.0555803583921,
783:   0.00705963113955,
784:   0.23235780058,
785:   0.705002209888,
786:   0.422442204032,
787:   0.422442204032,
788:   0.0821215678634,
789:   0.0729940240731,
790:   0.297574315013,
791:   0.297574315013,
792:   0.0578476039361,
793:   0.347003766038,
794:   0.134432268912,
795:   0.134432268912,
796:   0.0261332522867,
797:   0.705002209888,
798:   0.273713872823,
799:   0.273713872823,
800:   0.379578230281,
801:   0.0729940240731,
802:   0.192807956775,
803:   0.192807956775,
804:   0.267380320412,
805:   0.347003766038,
806:   0.0871029849888,
807:   0.0871029849888,
808:   0.120791820134,
809:   0.705002209888,
810:   0.0984204739396,
811:   0.0984204739396,
812:   0.730165028048,
813:   0.0729940240731,
814:   0.0693287858938,
815:   0.0693287858938,
816:   0.514338662174,
817:   0.347003766038,
818:   0.0313199947658,
819:   0.0313199947658,
820:   0.23235780058,
821:   0.705002209888,
822:   0.0952198798417,
823:   0.749664528222,
824:   0.0821215678634,
825:   0.0729940240731,
826:   0.0670742417521,
827:   0.528074388273,
828:   0.0578476039361,
829:   0.347003766038,
830:   0.0303014811743,
831:   0.23856305665,
832:   0.0261332522867,
833:   0.705002209888,
834:   0.0616960186091,
835:   0.485731727037,
836:   0.379578230281,
837:   0.0729940240731,
838:   0.0434595556538,
839:   0.342156357896,
840:   0.267380320412,
841:   0.347003766038,
842:   0.0196333029355,
843:   0.154572667042,
844:   0.120791820134,
845:   0.705002209888,
846:   0.0221843026408,
847:   0.174656645238,
848:   0.730165028048,
849:   0.0729940240731,
850:   0.0156269392579,
851:   0.12303063253,
852:   0.514338662174,
853:   0.347003766038,
854:   0.00705963113955,
855:   0.0555803583921,
856:   0.23235780058,
857:   0.705002209888};

859: /* Nodal basis function derivative evaluations */
860: static double BasisDerivatives[324] = {
861:   -0.5,
862:   -0.5,
863:   -0.5,
864:   0.5,
865:   2.6964953901e-17,
866:   8.15881875835e-17,
867:   0.0,
868:   0.5,
869:   1.52600313755e-17,
870:   0.0,
871:   0.0,
872:   0.5,
873:   -0.5,
874:   -0.5,
875:   -0.5,
876:   0.5,
877:   2.6938349868e-17,
878:   1.08228622783e-16,
879:   0.0,
880:   0.5,
881:   -1.68114298577e-17,
882:   0.0,
883:   0.0,
884:   0.5,
885:   -0.5,
886:   -0.5,
887:   -0.5,
888:   0.5,
889:   2.69035912409e-17,
890:   1.43034809879e-16,
891:   0.0,
892:   0.5,
893:   -5.87133459397e-17,
894:   0.0,
895:   0.0,
896:   0.5,
897:   -0.5,
898:   -0.5,
899:   -0.5,
900:   0.5,
901:   2.6964953901e-17,
902:   2.8079494593e-17,
903:   0.0,
904:   0.5,
905:   1.52600313755e-17,
906:   0.0,
907:   0.0,
908:   0.5,
909:   -0.5,
910:   -0.5,
911:   -0.5,
912:   0.5,
913:   2.6938349868e-17,
914:   7.0536336094e-17,
915:   0.0,
916:   0.5,
917:   -1.68114298577e-17,
918:   0.0,
919:   0.0,
920:   0.5,
921:   -0.5,
922:   -0.5,
923:   -0.5,
924:   0.5,
925:   2.69035912409e-17,
926:   1.26006930238e-16,
927:   0.0,
928:   0.5,
929:   -5.87133459397e-17,
930:   0.0,
931:   0.0,
932:   0.5,
933:   -0.5,
934:   -0.5,
935:   -0.5,
936:   0.5,
937:   2.6964953901e-17,
938:   -3.49866380524e-17,
939:   0.0,
940:   0.5,
941:   1.52600313755e-17,
942:   0.0,
943:   0.0,
944:   0.5,
945:   -0.5,
946:   -0.5,
947:   -0.5,
948:   0.5,
949:   2.6938349868e-17,
950:   2.61116525673e-17,
951:   0.0,
952:   0.5,
953:   -1.68114298577e-17,
954:   0.0,
955:   0.0,
956:   0.5,
957:   -0.5,
958:   -0.5,
959:   -0.5,
960:   0.5,
961:   2.69035912409e-17,
962:   1.05937620823e-16,
963:   0.0,
964:   0.5,
965:   -5.87133459397e-17,
966:   0.0,
967:   0.0,
968:   0.5,
969:   -0.5,
970:   -0.5,
971:   -0.5,
972:   0.5,
973:   2.6964953901e-17,
974:   1.52807111565e-17,
975:   0.0,
976:   0.5,
977:   1.52600313755e-17,
978:   0.0,
979:   0.0,
980:   0.5,
981:   -0.5,
982:   -0.5,
983:   -0.5,
984:   0.5,
985:   2.6938349868e-17,
986:   6.1520690456e-17,
987:   0.0,
988:   0.5,
989:   -1.68114298577e-17,
990:   0.0,
991:   0.0,
992:   0.5,
993:   -0.5,
994:   -0.5,
995:   -0.5,
996:   0.5,
997:   2.69035912409e-17,
998:   1.21934019246e-16,
999:   0.0,
1000:   0.5,
1001:   -5.87133459397e-17,
1002:   0.0,
1003:   0.0,
1004:   0.5,
1005:   -0.5,
1006:   -0.5,
1007:   -0.5,
1008:   0.5,
1009:   2.6964953901e-17,
1010:   -1.48832491782e-17,
1011:   0.0,
1012:   0.5,
1013:   1.52600313755e-17,
1014:   0.0,
1015:   0.0,
1016:   0.5,
1017:   -0.5,
1018:   -0.5,
1019:   -0.5,
1020:   0.5,
1021:   2.6938349868e-17,
1022:   4.0272766482e-17,
1023:   0.0,
1024:   0.5,
1025:   -1.68114298577e-17,
1026:   0.0,
1027:   0.0,
1028:   0.5,
1029:   -0.5,
1030:   -0.5,
1031:   -0.5,
1032:   0.5,
1033:   2.69035912409e-17,
1034:   1.1233505023e-16,
1035:   0.0,
1036:   0.5,
1037:   -5.87133459397e-17,
1038:   0.0,
1039:   0.0,
1040:   0.5,
1041:   -0.5,
1042:   -0.5,
1043:   -0.5,
1044:   0.5,
1045:   2.6964953901e-17,
1046:   -5.04349365259e-17,
1047:   0.0,
1048:   0.5,
1049:   1.52600313755e-17,
1050:   0.0,
1051:   0.0,
1052:   0.5,
1053:   -0.5,
1054:   -0.5,
1055:   -0.5,
1056:   0.5,
1057:   2.6938349868e-17,
1058:   1.52296507396e-17,
1059:   0.0,
1060:   0.5,
1061:   -1.68114298577e-17,
1062:   0.0,
1063:   0.0,
1064:   0.5,
1065:   -0.5,
1066:   -0.5,
1067:   -0.5,
1068:   0.5,
1069:   2.69035912409e-17,
1070:   1.01021564153e-16,
1071:   0.0,
1072:   0.5,
1073:   -5.87133459397e-17,
1074:   0.0,
1075:   0.0,
1076:   0.5,
1077:   -0.5,
1078:   -0.5,
1079:   -0.5,
1080:   0.5,
1081:   2.6964953901e-17,
1082:   -5.10267652705e-17,
1083:   0.0,
1084:   0.5,
1085:   1.52600313755e-17,
1086:   0.0,
1087:   0.0,
1088:   0.5,
1089:   -0.5,
1090:   -0.5,
1091:   -0.5,
1092:   0.5,
1093:   2.6938349868e-17,
1094:   1.4812758129e-17,
1095:   0.0,
1096:   0.5,
1097:   -1.68114298577e-17,
1098:   0.0,
1099:   0.0,
1100:   0.5,
1101:   -0.5,
1102:   -0.5,
1103:   -0.5,
1104:   0.5,
1105:   2.69035912409e-17,
1106:   1.00833228612e-16,
1107:   0.0,
1108:   0.5,
1109:   -5.87133459397e-17,
1110:   0.0,
1111:   0.0,
1112:   0.5,
1113:   -0.5,
1114:   -0.5,
1115:   -0.5,
1116:   0.5,
1117:   2.6964953901e-17,
1118:   -5.78459929494e-17,
1119:   0.0,
1120:   0.5,
1121:   1.52600313755e-17,
1122:   0.0,
1123:   0.0,
1124:   0.5,
1125:   -0.5,
1126:   -0.5,
1127:   -0.5,
1128:   0.5,
1129:   2.6938349868e-17,
1130:   1.00091968699e-17,
1131:   0.0,
1132:   0.5,
1133:   -1.68114298577e-17,
1134:   0.0,
1135:   0.0,
1136:   0.5,
1137:   -0.5,
1138:   -0.5,
1139:   -0.5,
1140:   0.5,
1141:   2.69035912409e-17,
1142:   9.86631702223e-17,
1143:   0.0,
1144:   0.5,
1145:   -5.87133459397e-17,
1146:   0.0,
1147:   0.0,
1148:   0.5,
1149:   -0.5,
1150:   -0.5,
1151:   -0.5,
1152:   0.5,
1153:   2.6964953901e-17,
1154:   -6.58832349994e-17,
1155:   0.0,
1156:   0.5,
1157:   1.52600313755e-17,
1158:   0.0,
1159:   0.0,
1160:   0.5,
1161:   -0.5,
1162:   -0.5,
1163:   -0.5,
1164:   0.5,
1165:   2.6938349868e-17,
1166:   4.34764891191e-18,
1167:   0.0,
1168:   0.5,
1169:   -1.68114298577e-17,
1170:   0.0,
1171:   0.0,
1172:   0.5,
1173:   -0.5,
1174:   -0.5,
1175:   -0.5,
1176:   0.5,
1177:   2.69035912409e-17,
1178:   9.61055074835e-17,
1179:   0.0,
1180:   0.5,
1181:   -5.87133459397e-17,
1182:   0.0,
1183:   0.0,
1184:   0.5};

1186: #endif

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

1197:   if (debug) {
1198:     MPI_Comm comm = mesh->comm();
1199:     int      rank = mesh->commRank();

1201:     PetscSynchronizedPrintf(comm, "[%d]Element (%d, %d)\n", rank, e.prefix, e.index);
1202:     PetscSynchronizedPrintf(comm, "[%d]Coordinates:\n[%d]  ", rank, rank);
1203:     for(int f = 0; f <= dim; f++) {
1204:       PetscSynchronizedPrintf(comm, " (");
1205:       for(int d = 0; d < dim; d++) {
1206:         if (d > 0) PetscSynchronizedPrintf(comm, ", ");
1207:         PetscSynchronizedPrintf(comm, "%g", coords[f*dim+d]);
1208:       }
1209:       PetscSynchronizedPrintf(comm, ")");
1210:     }
1211:     PetscSynchronizedPrintf(comm, "\n");
1212:   }
1213:   if (v0) {
1214:     for(int d = 0; d < dim; d++) {
1215:       v0[d] = coords[d];
1216:     }
1217:   }
1218:   if (J) {
1219:     for(int d = 0; d < dim; d++) {
1220:       for(int f = 0; f < dim; f++) {
1221:         J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
1222:       }
1223:     }
1224:     if (debug) {
1225:       MPI_Comm comm = mesh->comm();
1226:       int      rank = mesh->commRank();

1228:       for(int d = 0; d < dim; d++) {
1229:         if (d == 0) {
1230:           PetscSynchronizedPrintf(comm, "[%d]J = /", rank);
1231:         } else if (d == dim-1) {
1232:           PetscSynchronizedPrintf(comm, "[%d]    \\", rank);
1233:         } else {
1234:           PetscSynchronizedPrintf(comm, "[%d]    |", rank);
1235:         }
1236:         for(int e = 0; e < dim; e++) {
1237:           PetscSynchronizedPrintf(comm, " %g", J[d*dim+e]);
1238:         }
1239:         if (d == 0) {
1240:           PetscSynchronizedPrintf(comm, " \\\n");
1241:         } else if (d == dim-1) {
1242:           PetscSynchronizedPrintf(comm, " /\n");
1243:         } else {
1244:           PetscSynchronizedPrintf(comm, " |\n");
1245:         }
1246:       }
1247:     }
1248:     if (dim == 2) {
1249:       det = J[0]*J[3] - J[1]*J[2];
1250:     } else if (dim == 3) {
1251:       det = J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
1252:             J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
1253:             J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1254:     }
1255:     invDet = 1.0/det;
1256:     if (detJ) {
1257:       if (det < 0) {SETERRQ(PETSC_ERR_ARG_WRONG, "Negative Jacobian determinant");}
1258:       *detJ = det;
1259:     }
1260:     if (invJ) {
1261:       if (dim == 2) {
1262:         invJ[0] =  invDet*J[3];
1263:         invJ[1] = -invDet*J[1];
1264:         invJ[2] = -invDet*J[2];
1265:         invJ[3] =  invDet*J[0];
1266:       } else if (dim == 3) {
1267:         // FIX: This may be wrong
1268:         invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
1269:         invJ[0*3+1] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
1270:         invJ[0*3+2] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1271:         invJ[1*3+0] = invDet*(J[0*3+1]*J[2*3+2] - J[0*3+2]*J[2*3+1]);
1272:         invJ[1*3+1] = invDet*(J[0*3+2]*J[2*3+0] - J[0*3+0]*J[2*3+2]);
1273:         invJ[1*3+2] = invDet*(J[0*3+0]*J[2*3+1] - J[0*3+1]*J[2*3+0]);
1274:         invJ[2*3+0] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
1275:         invJ[2*3+1] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
1276:         invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
1277:       }
1278:       if (debug) {
1279:         MPI_Comm comm = mesh->comm();
1280:         int      rank = mesh->commRank();

1282:         for(int d = 0; d < dim; d++) {
1283:           if (d == 0) {
1284:             PetscSynchronizedPrintf(comm, "[%d]Jinv = /", rank);
1285:           } else if (d == dim-1) {
1286:             PetscSynchronizedPrintf(comm, "[%d]       \\", rank);
1287:           } else {
1288:             PetscSynchronizedPrintf(comm, "[%d]       |", rank);
1289:           }
1290:           for(int e = 0; e < dim; e++) {
1291:             PetscSynchronizedPrintf(comm, " %g", invJ[d*dim+e]);
1292:           }
1293:           if (d == 0) {
1294:             PetscSynchronizedPrintf(comm, " \\\n");
1295:           } else if (d == dim-1) {
1296:             PetscSynchronizedPrintf(comm, " /\n");
1297:           } else {
1298:             PetscSynchronizedPrintf(comm, " |\n");
1299:           }
1300:         }
1301:       }
1302:     }
1303:   }
1304:   return(0);
1305: }

1309: PetscErrorCode CheckElementGeometry(ALE::Obj<ALE::Two::Mesh> mesh)
1310: {
1311:   ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = mesh->getTopology()->heightStratum(0);
1312:   PetscInt       dim = mesh->getDimension();
1313:   PetscReal     *v0, *Jac;
1314:   PetscReal      detJ;

1318:   PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&Jac);
1319:   for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1320:     ElementGeometry(mesh, *e_iter, v0, Jac, PETSC_NULL, &detJ);
1321:   }
1322:   PetscSynchronizedFlush(mesh->comm());
1323:   PetscFree2(v0,Jac);
1324:   return(0);
1325: }

1329: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
1330: {
1331:   ALE::Obj<ALE::Two::Mesh> m;
1332:   Mesh                mesh = (Mesh) dmmg->dm;
1333:   UserContext        *user = (UserContext *) dmmg->user;
1334:   MPI_Comm            comm;
1335:   PetscReal           elementVec[NUM_BASIS_FUNCTIONS];
1336:   PetscReal           *v0, *Jac;
1337:   PetscReal           xi, eta, x_q, y_q, detJ, funcValue;
1338:   PetscInt            dim;
1339:   PetscInt            f, q;
1340:   PetscErrorCode      ierr;

1343:   PetscObjectGetComm((PetscObject) mesh, &comm);
1344:   MeshGetMesh(mesh, &m);
1345:   dim  = m->getDimension();
1346:   PetscMalloc(dim * sizeof(PetscReal), &v0);
1347:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1348:   ALE::Obj<ALE::Two::Mesh::field_type> field = m->getField("b");
1349:   ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
1350:   ALE::Two::Mesh::field_type::patch_type patch;
1351:   for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
1352:     ElementGeometry(m, *e_itor, v0, Jac, PETSC_NULL, &detJ);
1353:     /* Element integral */
1354:     PetscMemzero(elementVec, NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
1355:     for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1356:       xi = points[q*2+0] + 1.0;
1357:       eta = points[q*2+1] + 1.0;
1358:       x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1359:       y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1360:       funcValue = 2.0;
1361:       //funcValue = 2.0*y_q*(2.0 - y_q) + 2.0*x_q*(2.0 - x_q);
1362:       for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1363:         elementVec[f] += Basis[q*NUM_BASIS_FUNCTIONS+f]*funcValue*weights[q]*detJ;
1364:       }
1365:     }
1366:     if (debug) {PetscSynchronizedPrintf(comm, "elementVec = [%g %g %g]\n", elementVec[0], elementVec[1], elementVec[2]);}
1367:     /* Assembly */
1368:     field->updateAdd("element", *e_itor, elementVec);
1369:     if (debug) {PetscSynchronizedFlush(comm);}
1370:   }
1371:   PetscFree(v0);
1372:   PetscFree(Jac);

1374:   Vec locB;
1375:   VecCreateSeqWithArray(PETSC_COMM_SELF, field->getSize(patch), field->restrict(patch), &locB);
1376:   VecScatterBegin(locB, b, ADD_VALUES, SCATTER_FORWARD, user->injection);
1377:   VecScatterEnd(locB, b, ADD_VALUES, SCATTER_FORWARD, user->injection);
1378:   VecDestroy(locB);

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

1385:     KSPGetNullSpace(dmmg->ksp,&nullspace);
1386:     MatNullSpaceRemove(nullspace,b,PETSC_NULL);
1387:   }
1388:   if (user->bcType == DIRICHLET) {
1389:     /* Zero out BC rows */
1390:     ALE::Two::Mesh::field_type::patch_type patch;
1391:     ALE::Two::Mesh::field_type::patch_type bdPatch(0, 1);
1392:     ALE::Obj<ALE::Two::Mesh::field_type> boundary = m->getBoundary();
1393:     ALE::Obj<ALE::Two::Mesh::field_type::order_type::coneSequence> cone = boundary->getPatch(bdPatch);
1394:     PetscScalar *boundaryValues;
1395:     PetscInt    *boundaryIndices;
1396:     PetscInt     numBoundaryIndices = 0;
1397:     PetscInt     k = 0;

1399:     for(ALE::Two::Mesh::field_type::order_type::coneSequence::iterator p = cone->begin(); p != cone->end(); ++p) {
1400:       numBoundaryIndices += field->getGlobalOrder()->getIndex(patch, *p).index;
1401:     }
1402:     PetscMalloc2(numBoundaryIndices,PetscInt,&boundaryIndices,numBoundaryIndices,PetscScalar,&boundaryValues);
1403:     for(ALE::Two::Mesh::field_type::order_type::coneSequence::iterator p = cone->begin(); p != cone->end(); ++p) {
1404:       const ALE::Two::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);
1405:       const double *data = boundary->restrict(bdPatch, *p);

1407:       for(int i = 0; i < idx.index; i++) {
1408:         boundaryIndices[k] = idx.prefix + i;
1409:         boundaryValues[k] = data[i];
1410:         k++;
1411:       }
1412:     }
1413:     if (debug) {
1414:       boundary->view("Boundary for rhs conditions");
1415:       for(int i = 0; i < numBoundaryIndices; i++) {
1416:         PetscSynchronizedPrintf(comm, "[%d]boundaryIndices[%d] = %d\n", m->commRank(), i, boundaryIndices[i]);
1417:       }
1418:     }
1419:     PetscSynchronizedFlush(comm);
1420:     VecSetValues(b, numBoundaryIndices, boundaryIndices, boundaryValues, INSERT_VALUES);
1421:     PetscFree2(boundaryIndices, boundaryValues);
1422:   }
1423:   if (debug) {
1424:     PetscPrintf(comm, "Rhs vector:");
1425:     VecView(b, PETSC_VIEWER_STDOUT_WORLD);
1426:   }
1427:   return(0);
1428: }

1432: PetscErrorCode ComputeJacobian(DMMG dmmg, Mat J, Mat jac)
1433: {
1434:   ALE::Obj<ALE::Two::Mesh> m;
1435:   Mesh              mesh = (Mesh) dmmg->dm;
1436:   UserContext      *user = (UserContext *) dmmg->user;
1437:   MPI_Comm          comm;
1438:   PetscReal         elementMat[NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS];
1439:   PetscReal        *v0, *Jac, *Jinv, *t_der, *b_der;
1440:   PetscReal         xi, eta, x_q, y_q, detJ;
1441:   PetscInt          dim;
1442:   PetscInt          f, g, q;
1443:   PetscMPIInt       rank;
1444:   PetscErrorCode    ierr;

1447:   PetscObjectGetComm((PetscObject) mesh, &comm);
1448:   MPI_Comm_rank(comm, &rank);
1449:   MeshGetMesh(mesh, &m);
1450:   dim  = m->getDimension();
1451:   PetscMalloc(dim * sizeof(PetscReal), &v0);
1452:   PetscMalloc(dim * sizeof(PetscReal), &t_der);
1453:   PetscMalloc(dim * sizeof(PetscReal), &b_der);
1454:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1455:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jinv);
1456:   ALE::Obj<ALE::Two::Mesh::field_type> field = m->getField("u");
1457:   ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
1458:   for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
1459:     CHKMEMQ;
1460:     ElementGeometry(m, *e_itor, v0, Jac, Jinv, &detJ);
1461:     /* Element integral */
1462:     PetscMemzero(elementMat, NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
1463:     for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1464:       xi = points[q*2+0] + 1.0;
1465:       eta = points[q*2+1] + 1.0;
1466:       x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1467:       y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1468:       for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1469:         t_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
1470:         t_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
1471:         for(g = 0; g < NUM_BASIS_FUNCTIONS; g++) {
1472:           b_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
1473:           b_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
1474:           elementMat[f*NUM_BASIS_FUNCTIONS+g] += (t_der[0]*b_der[0] + t_der[1]*b_der[1])*weights[q]*detJ;
1475:         }
1476:       }
1477:     }
1478:     if (debug) {
1479:       PetscSynchronizedPrintf(comm, "[%d]elementMat = [%g %g %g]\n                [%g %g %g]\n                [%g %g %g]\n",
1480:                                      rank, elementMat[0], elementMat[1], elementMat[2], elementMat[3], elementMat[4],
1481:                                      elementMat[5], elementMat[6], elementMat[7], elementMat[8]);
1482:     }
1483:     /* Assembly */
1484:     updateOperator(jac, field, *e_itor, elementMat, ADD_VALUES);
1485:     if (debug) {PetscSynchronizedFlush(comm);}
1486:   }
1487:   PetscFree(v0);
1488:   PetscFree(t_der);
1489:   PetscFree(b_der);
1490:   PetscFree(Jac);
1491:   PetscFree(Jinv);
1492:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
1493:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);

1495:   if (user->bcType == DIRICHLET) {
1496:     /* Zero out BC rows */
1497:     ALE::Two::Mesh::field_type::patch_type patch;
1498:     ALE::Two::Mesh::field_type::patch_type bdPatch(0, 1);
1499:     ALE::Obj<ALE::Two::Mesh::field_type> boundary = m->getBoundary();
1500:     ALE::Obj<ALE::Two::Mesh::field_type::order_type::coneSequence> cone = boundary->getPatch(bdPatch);
1501:     PetscInt *boundaryIndices;
1502:     PetscInt  numBoundaryIndices = 0;
1503:     PetscInt  k = 0;

1505:     for(ALE::Two::Mesh::field_type::order_type::coneSequence::iterator p = cone->begin(); p != cone->end(); ++p) {
1506:       numBoundaryIndices += field->getGlobalOrder()->getIndex(patch, *p).index;
1507:     }
1508:     PetscMalloc(numBoundaryIndices * sizeof(PetscInt), &boundaryIndices);
1509:     for(ALE::Two::Mesh::field_type::order_type::coneSequence::iterator p = cone->begin(); p != cone->end(); ++p) {
1510:       const ALE::Two::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);

1512:       for(int i = 0; i < idx.index; i++) {
1513:         boundaryIndices[k++] = idx.prefix + i;
1514:       }
1515:     }
1516:     if (debug) {
1517:       for(int i = 0; i < numBoundaryIndices; i++) {
1518:         PetscSynchronizedPrintf(comm, "[%d]boundaryIndices[%d] = %d\n", rank, i, boundaryIndices[i]);
1519:       }
1520:     }
1521:     PetscSynchronizedFlush(comm);
1522:     MatZeroRows(jac, numBoundaryIndices, boundaryIndices, 1.0);
1523:     PetscFree(boundaryIndices);
1524:   }
1525:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
1526:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
1527:   return(0);
1528: }

1532: PetscErrorCode ComputeError(DMMG dmmg, Vec u, double *error)
1533: {
1534:   ALE::Obj<ALE::Two::Mesh> m;
1535:   Mesh                mesh = (Mesh) dmmg->dm;
1536:   UserContext        *user = (UserContext *) dmmg->user;
1537:   MPI_Comm            comm;
1538:   PetscReal           *v0, *Jac;
1539:   PetscReal           detJ, totalError;
1540:   PetscInt            dim;
1541:   PetscInt            f, q;
1542:   PetscErrorCode      ierr;

1545:   PetscObjectGetComm((PetscObject) mesh, &comm);
1546:   MeshGetMesh(mesh, &m);
1547:   dim  = m->getDimension();
1548:   PetscMalloc(dim * sizeof(PetscReal), &v0);
1549:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1550:   ALE::Obj<ALE::Two::Mesh::field_type> field = m->getField("u");
1551:   ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
1552:   ALE::Two::Mesh::field_type::patch_type patch;

1554:   Vec locU;
1555:   VecCreateSeqWithArray(PETSC_COMM_SELF, field->getSize(patch), field->restrict(patch), &locU);
1556:   VecScatterBegin(u, locU, INSERT_VALUES, SCATTER_REVERSE, user->injection);
1557:   VecScatterEnd(u, locU, INSERT_VALUES, SCATTER_REVERSE, user->injection);
1558:   VecDestroy(locU);

1560:   totalError = 0.0;
1561:   for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
1562:     const double *elementVec = field->restrict("element", *e_itor);
1563:     double e;

1565:     e = 0.0;
1566:     ElementGeometry(m, *e_itor, v0, Jac, PETSC_NULL, &detJ);
1567:     /* Element integral */
1568:     for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1569:       PetscReal xi, eta, x_q, y_q, uExact, uComputed;

1571:       xi = points[q*2+0] + 1.0;
1572:       eta = points[q*2+1] + 1.0;
1573:       x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1574:       y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1575:       uExact = x_q*(2.0 - x_q);
1576:       //uExact = x_q*(2.0 - x_q)*y_q*(2.0 - y_q);
1577:       uComputed = 0.0;
1578:       for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1579:         uComputed += elementVec[f]*Basis[q*NUM_BASIS_FUNCTIONS+f];
1580:       }
1581:       e += (uComputed - uExact)*(uComputed - uExact)*weights[q]*detJ;
1582:     }
1583:     totalError += e;
1584:     if (debug) {PetscSynchronizedPrintf(comm, "elementError = %g\n", e);}
1585:     if (debug) {PetscSynchronizedFlush(comm);}
1586:   }
1587:   PetscFree(v0);
1588:   PetscFree(Jac);
1589:   *error = sqrt(totalError);
1590:   return(0);
1591: }