Actual source code: meshMovement.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: meshMovement.c,v 1.19 2000/10/17 13:48:57 knepley Exp $";
  3: #endif

  5: #include "src/gvec/meshMovementImpl.h" /*I "meshMovement.h" I*/
  6: #include "src/grid/gridimpl.h"    /*I "grid.h" I*//*I "gvec.h" I*/

  8: extern int MeshPreCopy_Private(Mesh);
  9: extern int MeshPostCopy_Private(Mesh, Mesh);

 11: /*-------------------------------------------------- Mesh Functions --------------------------------------------------*/
 12: #undef  __FUNCT__
 14: /*@
 15:   MeshSetMover - This function provides the MeshMover object for this mesh.

 17:   Not collective

 19:   Input Parameters:
 20: + mesh  - The mesh
 21: - mover - The mover, or PETSC_NULL

 23:   Level: intermediate

 25: .keywords: Mesh, movement, mover
 26: .seealso: MeshGetMover(), MeshSetMovement()
 27: @*/
 28: int MeshSetMover(Mesh mesh, MeshMover mover) {

 33:   PetscObjectCompose((PetscObject) mesh, "MeshMover", (PetscObject) mover);
 34:   if (mover != PETSC_NULL) {
 36:     PetscObjectDereference((PetscObject) mover);
 37:     MeshMoverSetMesh(mover, mesh);
 38:   }
 39:   return(0);
 40: }

 42: #undef  __FUNCT__
 44: /*@
 45:   MeshGetMover - This function returns the MeshMover object for this mesh.

 47:   Not collective

 49:   Input Parameter:
 50: . mesh  - The mesh

 52:   Output Parameters:
 53: . mover - The mover, or PETSC_NULL if it does not exist

 55:   Level: intermediate

 57: .keywords: Mesh, movement, mover
 58: .seealso: MeshSetMover(), MeshGetMovement()
 59: @*/
 60: int MeshGetMover(Mesh mesh, MeshMover *mover) {

 65:   PetscObjectQuery((PetscObject) mesh, "MeshMover", (PetscObject *) mover);
 66:   return(0);
 67: }

 69: #undef  __FUNCT__
 71: /*@
 72:   MeshMoveMesh - Recalculates the node coordinates based on the
 73:   current velocities and accelerations.

 75:   Collective on Mesh

 77:   Input Parameters:
 78: . mesh - The mesh
 79: . dt   - The timestep

 81:   Level: advanced

 83: .keywords: mesh, movement
 84: .seealso: MeshCalcNodeVelocitiess(), MeshCalcNodeAccelerations()
 85: @*/
 86: int MeshMoveMesh(Mesh mesh, double dt)
 87: {
 88:   MeshMover  mover;
 89:   PetscTruth isMoving;
 90:   int        ierr;


 95:   MeshGetMovement(mesh, &isMoving);
 96:   if (isMoving == PETSC_FALSE) return(0);
 97:   PetscLogEventBegin(MESH_MoveMesh, mesh, 0, 0, 0);
 98:   MeshGetMover(mesh, &mover);
 99:   (*mover->ops->movemesh)(mover, dt);
100:   PetscLogEventEnd(MESH_MoveMesh, mesh, 0, 0, 0);
101:   return(0);
102: }

104: #undef  __FUNCT__
106: /*@
107:   MeshUpdateNodeValues - Updates the mesh velocity, acceleration, or any other node-based quantity.

109:   Collective on Mesh

111:   Input Parameters:
112: + mesh     - The mesh
113: - sol      - The solution to an auxilliary problem

115:   Output Parameters:
116: + vec      - The vector which will hold he solution
117: - ghostVec - The corresponding ghost vector for vec

119:   Level: developer

121: .keywords: mesh, movement
122: .seealso: MeshMoverSetNodeVelocities()
123: @*/
124: int MeshUpdateNodeValues(Mesh mesh, GVec sol, Vec vec, Vec ghostVec)
125: {
126:   MeshMover mover;
127:   int       ierr;

134:   MeshGetMover(mesh, &mover);
135:   (*mover->ops->updatenodevalues)(mover, sol, vec, ghostVec);
136:   return(0);
137: }

139: EXTERN_C_BEGIN
140: #undef  __FUNCT__
142: /*
143:   MeshPostReform_Private - Create a suitable MeshMover after reforming a Mesh.

145:   Collective on Mesh

147:   Input Parameters:
148: + mesh    - The mesh
149: - newMesh - The reformed mesh

151:   Level: developer

153: .keywords: mesh, refine, reform
154: .seealso: MeshReform()
155: */
156: int MeshPostReform_Private(Mesh mesh, Mesh newMesh)
157: {
158:   MeshMover mover, newMover;
159:   int       ierr;

164:   MeshGetMover(mesh, &mover);
165:   if (mover != PETSC_NULL) {
166:     MeshMoverConvert(mover, newMesh, &newMover);
167:   }
168:   return(0);
169: }
170: EXTERN_C_END

172: #undef  __FUNCT__
174: /*@
175:   MeshReform - Reform a mesh using the original boundary.

177:   Collective on Mesh

179:   Input Parameters:
180: + mesh    - The previous mesh
181: . refine  - Flag indicating whether to refine or recreate the mesh
182: . area    - A function which gives an area constraint when evaluated inside an element
183: - newBd   - Flag to determine whether the boundary should be generated (PETSC_TRUE) or taken from storage

185:   Output Parameter:
186: . newMesh - The new mesh

188:   Level: beginner

190: .keywords: mesh, refinement, reform
191: .seealso: MeshReform(), MeshSetBoundary()
192: @*/
193: int MeshReform(Mesh mesh, PetscTruth refine, PointFunction area, PetscTruth newBd, Mesh *newMesh)
194: {
195:   PetscObject obj = (PetscObject) mesh;
196:   Mesh        refinedMesh;
197:   MeshMover   mover;
198:   PetscReal   maxArea, startX, endX, startY, endY;
199:   void       *tempCtx;
200:   int         ierr;

205:   PetscObjectComposeFunction(obj, "PostCopy", "MeshPostReform_Private", (void (*)(void)) MeshPostReform_Private);
206:   if (refine == PETSC_TRUE) {
207:     /* Refine grid with tolerance equal to the entire size
208:        This just "flips" the current mesh and does the necessary point insertion
209:     */
210:     MeshGetBoundingBox(mesh, &startX, &startY, PETSC_NULL, &endX, &endY, PETSC_NULL);
211:     maxArea = (endX - startX)*(endY - startY);
212:     MeshGetUserContext(mesh, &tempCtx);
213:     MeshSetUserContext(mesh, (void *) &maxArea);
214:     MeshRefine(mesh, PointFunctionConstant, newMesh);
215:     MeshSetUserContext(mesh, tempCtx);
216:     MeshSetUserContext(*newMesh, tempCtx);
217:   } else {
218:     MeshGetMover(mesh, &mover);
219:     if (mover->ops->reform == PETSC_NULL) SETERRQ(PETSC_ERR_SUP, " ");
220:     PetscLogEventBegin(MESH_Reform, mesh, 0, 0, 0);
221:     MeshPreCopy_Private(mesh);
222:     (*mover->ops->reform)(mesh, refine, area, newBd, newMesh);
223:     MeshPostCopy_Private(mesh, *newMesh);
224:     PetscLogEventEnd(MESH_Reform, mesh, 0, 0, 0);
225:   }

227:   /* Refine new mesh if necessary */
228:   if (area != PETSC_NULL) {
230:     PetscObjectComposeFunction((PetscObject) *newMesh, "PostCopy", "MeshPostReform_Private", (void (*)(void)) MeshPostReform_Private);
231:     MeshRefine(*newMesh, area, &refinedMesh);
232:     MeshDestroy(*newMesh);
233:     *newMesh = refinedMesh;
234:     PetscObjectComposeFunction((PetscObject) *newMesh, "PostCopy", "", PETSC_NULL);
235:   }
236:   PetscObjectComposeFunction(obj, "PostCopy", "", PETSC_NULL);

238:   return(0);
239: }

241: /*----------------------------------------------- MeshMover Functions ------------------------------------------------*/
242: #undef  __FUNCT__
244: /*@ 
245:   MeshMoverCreate - This function creates an empty mesh mover.

247:   Collective on Mesh

249:   Input Parameter:
250: . mesh  - The Mesh to be moved

252:   Output Parameter:
253: . mover - The mover

255:   Level: beginner

257: .keywords: Mesh, create, mover
258: .seealso: MeshSetUp(), MeshDestroy(), GridCreate()
259: @*/
260: int MeshMoverCreate(Mesh mesh, MeshMover *mover) {
261:   MeshMover m;
262:   MPI_Comm  comm;
263:   int       ierr;

267:   *mover = PETSC_NULL;

269:   PetscObjectGetComm((PetscObject) mesh, &comm);
270:   PetscHeaderCreate(m, _MeshMover, struct _MeshMoverOps, MESH_COOKIE, -1, "MeshMover", comm, MeshMoverDestroy, MeshMoverView);
271:   PetscLogObjectCreate(m);
272:   PetscLogObjectMemory(m, sizeof(struct _MeshMover));
273:   PetscMemzero(m->ops, sizeof(struct _MeshMoverOps));
274:   m->bops->publish  = PETSC_NULL /* MeshMoverPublish_Petsc */;
275:   m->type_name      = PETSC_NULL;
276:   m->serialize_name = PETSC_NULL;

278:   m->meshVelocityMethod      = MESH_LAPLACIAN;
279:   m->meshVelocityGrid        = PETSC_NULL;
280:   m->nodeVelocities          = PETSC_NULL;
281:   m->nodeVelocitiesGhost     = PETSC_NULL;
282:   m->meshVelocityMat         = PETSC_NULL;
283:   m->meshVelocityRhs         = PETSC_NULL;
284:   m->meshVelocitySol         = PETSC_NULL;
285:   m->meshVelocitySles        = PETSC_NULL;
286:   m->meshAccelerationMethod  = MESH_LAPLACIAN;
287:   m->meshAccelerationGrid    = PETSC_NULL;
288:   m->nodeAccelerations       = PETSC_NULL;
289:   m->nodeAccelerationsGhost  = PETSC_NULL;
290:   m->meshAccelerationMat     = PETSC_NULL;
291:   m->meshAccelerationRhs     = PETSC_NULL;
292:   m->meshAccelerationSol     = PETSC_NULL;
293:   m->meshAccelerationSles    = PETSC_NULL;
294:   m->meshVelocityScatter     = PETSC_NULL;
295:   m->movingMeshViewer        = PETSC_NULL;
296:   m->movingMeshViewerCaption = "Moving Mesh";
297:   m->movingMeshCtx           = PETSC_NULL;

299:   MeshSetMover(mesh, m);
300:   *mover = m;
301:   return(0);
302: }

304: #undef  __FUNCT__
306: /*@
307:   MeshMoverSetup - Creates all the necessary structures for mesh movement.

309:   Collective on MeshMover

311:   Input Parameter:
312: . mover - The mover

314:   Note:
315:   This function should be called directly before beginning the
316:   calculation, and after all options have been set. It is
317:   normally called automatically by the solver.

319:   Level: advanced

321: .keywords: mesh, movement
322: .seealso: MeshMoverGetMovement()
323: @*/
324: int MeshMoverSetup(MeshMover mover)
325: {
326:   Mesh       mesh;
327:   Partition  part;
328:   IS         globalIS, localIS;
329:   int        dim;
330:   int       *indices;
331:   int        numNodes, numOverlapNodes;
332:   int        node, gNode, j;
333:   PetscTruth reduceVel, reduceAcc;
334:   PetscTruth opt;
335:   int        ierr;

339:   /* Destroy old structures */
340:   if (mover->nodeVelocities != PETSC_NULL) {
341:     VecDestroy(mover->nodeVelocities);
342:     VecDestroy(mover->nodeVelocitiesGhost);
343:   }
344:   if (mover->meshVelocityGrid != PETSC_NULL) {
345:     GridDestroy(mover->meshVelocityGrid);
346:     mover->meshVelocityGrid = PETSC_NULL;
347:     GMatDestroy(mover->meshVelocityMat);
348:     mover->meshVelocityMat  = PETSC_NULL;
349:     GVecDestroy(mover->meshVelocityRhs);
350:     mover->meshVelocityRhs  = PETSC_NULL;
351:     GVecDestroy(mover->meshVelocitySol);
352:     mover->meshVelocitySol  = PETSC_NULL;
353:     SLESDestroy(mover->meshVelocitySles);
354:     mover->meshVelocitySles = PETSC_NULL;
355:   }
356:   if (mover->nodeAccelerations != PETSC_NULL) {
357:     VecDestroy(mover->nodeAccelerations);
358:     VecDestroy(mover->nodeAccelerationsGhost);
359:   }
360:   if (mover->meshAccelerationGrid != PETSC_NULL) {
361:     GridDestroy(mover->meshAccelerationGrid);
362:     mover->meshAccelerationGrid = PETSC_NULL;
363:     GMatDestroy(mover->meshAccelerationMat);
364:     mover->meshAccelerationMat  = PETSC_NULL;
365:     GVecDestroy(mover->meshAccelerationRhs);
366:     mover->meshAccelerationRhs  = PETSC_NULL;
367:     GVecDestroy(mover->meshAccelerationSol);
368:     mover->meshAccelerationSol  = PETSC_NULL;
369:     SLESDestroy(mover->meshAccelerationSles);
370:     mover->meshAccelerationSles = PETSC_NULL;
371:   }
372:   if (mover->meshVelocityScatter != PETSC_NULL) {
373:     VecScatterDestroy(mover->meshVelocityScatter);
374:   }
375:   /* Determine solution method */
376:   ierr  = PetscOptionsHasName(mover->prefix, "-mesh_move_lap", &opt);
377:   if (opt == PETSC_TRUE) {
378:     mover->meshVelocityMethod     = MESH_LAPLACIAN;
379:     mover->meshAccelerationMethod = MESH_LAPLACIAN;
380:   }
381:   ierr  = PetscOptionsHasName(mover->prefix, "-mesh_move_weighted_lap", &opt);
382:   if (opt == PETSC_TRUE) {
383:     mover->meshVelocityMethod     = MESH_WEIGHTED_LAPLACIAN;
384:     mover->meshAccelerationMethod = MESH_WEIGHTED_LAPLACIAN;
385:   }
386:   /* Setup mesh node velocities vector */
387:   MeshMoverGetMesh(mover, &mesh);
388:   MeshGetPartition(mesh, &part);
389:   MeshGetDimension(mesh, &dim);
390:   MeshGetInfo(mesh, PETSC_NULL, &numNodes, PETSC_NULL, PETSC_NULL);
391:   VecCreateMPI(mover->comm, numNodes*dim, PETSC_DECIDE, &mover->nodeVelocities);
392:   PartitionGetNumOverlapNodes(part, &numOverlapNodes);
393:   VecCreateSeq(PETSC_COMM_SELF, numOverlapNodes*dim, &mover->nodeVelocitiesGhost);
394:   /* Setup mesh node accelerations vector */
395:   VecCreateMPI(mover->comm, numNodes*dim, PETSC_DECIDE, &mover->nodeAccelerations);
396:   VecCreateSeq(PETSC_COMM_SELF, numOverlapNodes*dim, &mover->nodeAccelerationsGhost);
397:   /* Setup grids and other implementation specific stuff */
398:   (*mover->ops->setup)(mover);
399:   /* Setup scatter from node values vector to ghost vector */
400:   ISCreateStride(mover->comm, numOverlapNodes*dim, 0, 1, &localIS);
401:   PetscMalloc(numOverlapNodes*dim * sizeof(int), &indices);
402:   for(node = 0; node < numOverlapNodes; node++) {
403:     PartitionLocalToGlobalNodeIndex(part, node, &gNode);
404:     for(j = 0; j < dim; j++)
405:       indices[node*dim+j] = gNode*dim+j;
406:   }
407:   ISCreateGeneral(mover->comm, numOverlapNodes*dim, indices, &globalIS);
408:   VecScatterCreate(mover->nodeVelocities, globalIS, mover->nodeVelocitiesGhost, localIS, &mover->meshVelocityScatter);
409: 
410:   ISDestroy(localIS);
411:   ISDestroy(globalIS);
412:   PetscFree(indices);
413:   /* Setup viewer */
414:   PetscOptionsHasName(mover->prefix, "-mesh_moving_view", &opt);
415:   if ((opt == PETSC_TRUE) && (mover->movingMeshViewer == PETSC_NULL)) {
416:     mover->movingMeshViewer = PETSC_VIEWER_STDOUT_WORLD;
417:   }
418:   PetscOptionsHasName(mover->prefix, "-mesh_moving_view_draw", &opt);
419:   if ((opt == PETSC_TRUE) && (mover->movingMeshViewer == PETSC_NULL)) {
420:     PetscReal startX, startY, endX, endY, sizeX, sizeY;
421:     int       xsize, ysize;

423:     ierr  = MeshGetBoundingBox(mesh, &startX, &startY, 0, &endX, &endY, 0);
424:     sizeX = endX - startX;
425:     sizeY = endY - startY;
426:     if (sizeX > sizeY) {
427:       ysize  = 300;
428:       xsize  = ysize * (int) (sizeX/sizeY);
429:     } else {
430:       xsize  = 300;
431:       ysize  = xsize * (int) (sizeY/sizeX);
432:     }
433:     PetscViewerDrawOpen(mover->comm, 0, 0, 0, 315, xsize, ysize, &mover->movingMeshViewer);
434:   }
435:   /* Check grids */
436:   GridGetReduceSystem(mover->meshVelocityGrid,     &reduceVel);
437:   GridGetReduceSystem(mover->meshAccelerationGrid, &reduceAcc);
438:   if (reduceVel != reduceAcc) {
439:     SETERRQ(PETSC_ERR_ARG_WRONG, "Velocity and acceleration grids must have same BC reduction");
440:   }
441:   return(0);
442: }

444: #undef __FUNCT__  
446: /*
447:   MeshMoverSetTypeFromOptions - Sets the type of mesh movement from user options.

449:   Collective on MeshMover

451:   Input Parameter:
452: . mover - The mover

454:   Level: intermediate

456: .keywords: MeshMover, set, options, database, type
457: .seealso: MeshMoverSetFromOptions(), MeshMoverSetType()
458: */
459: static int MeshMoverSetTypeFromOptions(MeshMover mover)
460: {
461:   Mesh       mesh;
462:   PetscTruth opt;
463:   char      *defaultType;
464:   char       typeName[256];
465:   int        dim;
466:   int        ierr;

469:   MeshMoverGetMesh(mover, &mesh);
470:   MeshGetDimension(mesh, &dim);
471:   if (mover->type_name != PETSC_NULL) {
472:     defaultType = mover->type_name;
473:   } else {
474:     switch(dim)
475:     {
476:     case 2:
477:       defaultType = MESH_TRIANGULAR_2D;
478:       break;
479:     default:
480:       SETERRQ1(PETSC_ERR_SUP, "MeshMover dimension %d is not supported", dim);
481:     }
482:   }

484:   if (!MeshMoverRegisterAllCalled) {
485:     MeshMoverRegisterAll(PETSC_NULL);
486:   }
487:   PetscOptionsList("-mesh_mover_type", "Mesh movement method", "MeshMoverSetType", MeshMoverList, defaultType, typeName, 256, &opt);
488: 
489:   if (opt == PETSC_TRUE) {
490:     MeshMoverSetType(mover, typeName);
491:   } else {
492:     MeshMoverSetType(mover, defaultType);
493:   }
494:   return(0);
495: }

497: #undef __FUNCT__  
499: /*@
500:   MeshMoverSetFromOptions - Sets various MeshMover parameters from user options.

502:   Collective on MeshMover

504:   Input Parameter:
505: . mover - The mover

507:   Notes:  To see all options, run your program with the -help option, or consult the users manual.
508:           Must be called after MeshMoverCreate() but before the vector is used.

510:   Level: intermediate

512: .keywords: MeshMover, set, options, database
513: .seealso: MeshMoverCreate(), MeshMoverPrintHelp(), MeshMoverSetOptionsPrefix()
514: @*/
515: int MeshMoverSetFromOptions(MeshMover mover)
516: {
517:   PetscTruth opt;
518:   int        ierr;

522:   PetscOptionsBegin(mover->comm, mover->prefix, "MeshMover options", "MeshMover");

524:   /* Handle generic mesh options */
525:   PetscOptionsHasName(PETSC_NULL, "-help", &opt);
526:   if (opt == PETSC_TRUE) {
527:     MeshMoverPrintHelp(mover);
528:   }

530:   /* Handle mesh type options */
531:   MeshMoverSetTypeFromOptions(mover);

533:   /* Handle specific mesh options */
534:   if (mover->ops->setfromoptions != PETSC_NULL) {
535:     (*mover->ops->setfromoptions)(mover);
536:   }

538:   PetscOptionsEnd();

540:   /* Handle subobject options */
541:   return(0);
542: }

544: #undef __FUNCT__  
546: /*@
547:   MeshMoverPrintHelp - Prints all options for the MeshMover.

549:   Input Parameter:
550: . mover - The mover

552:   Options Database Keys:
553: $  -help, -h

555:   Level: intermediate

557: .keywords: MeshMover, help
558: .seealso: MeshMoverSetFromOptions()
559: @*/
560: int MeshMoverPrintHelp(MeshMover mover)
561: {
562:   char p[64];
563:   int  ierr;


568:   PetscStrcpy(p, "-");
569:   if (mover->prefix != PETSC_NULL) {
570:     PetscStrcat(p, mover->prefix);
571:   }

573:   (*PetscHelpPrintf)(mover->comm, "MeshMover options ------------------------------------------------n");
574:   (*PetscHelpPrintf)(mover->comm,"   %smesh_mover_type <typename> : Sets the mesh movement typen", p);
575:   return(0);
576: }

578: #undef  __FUNCT__
580: /*@
581:   MeshMoverDuplicate - This function returns an identical MeshMover.

583:   Not collective

585:   Input Parameter:
586: . mover    - The mover

588:   Output Parameter:
589: . newMover - The new mover

591:   Level: intermediate

593: .keywords: Mesh, movement, mover
594: .seealso: MeshMoverCopy()
595: @*/
596: int MeshMoverDuplicate(MeshMover mover, MeshMover *newMover) {
597:   Mesh mesh;
598:   int  ierr;

601:   MeshMoverGetMesh(mover, &mesh);
602:   MeshMoverConvert(mover, mesh, newMover);
603:   return(0);
604: }

606: #undef  __FUNCT__
608: /*@
609:   MeshMoverConvert - This function returns an identical MeshMover on the new Mesh.

611:   Not collective

613:   Input Parameters:
614: + mover   - The mover
615: - newMesh - The new Mesh

617:   Output Parameter:
618: . newMover - The new mover

620:   Level: intermediate

622: .keywords: Mesh, movement, mover
623: .seealso: MeshMoverCopy()
624: @*/
625: int MeshMoverConvert(MeshMover mover, Mesh newMesh, MeshMover *newMover) {
626:   MeshMover m;
627:   Mesh      mesh;
628:   int       ierr;

631:   MeshMoverGetMesh(mover, &mesh);
632:   if (mesh == newMesh) {
633:     /* Need to play reference games since MeshMoverCreate automatically sets it in the Mesh */
634:     PetscObjectReference((PetscObject) mover);
635:     MeshMoverCreate(newMesh, &m);
636:     PetscObjectReference((PetscObject) m);
637:     MeshSetMover(mesh, mover);
638:   } else {
639:     MeshMoverCreate(newMesh, &m);
640:   }

642:   MeshMoverSetFromOptions(m);
643:   m->meshVelocityMethod      = mover->meshVelocityMethod;
644:   m->meshAccelerationMethod  = mover->meshAccelerationMethod;
645:   m->movingMeshViewer        = mover->movingMeshViewer;
646:   if (mover->movingMeshViewer != PETSC_NULL) PetscObjectReference((PetscObject)  mover->movingMeshViewer);
647:   m->movingMeshViewerCaption = mover->movingMeshViewerCaption;
648:   m->movingMeshCtx           = mover->movingMeshCtx;
649:   MeshMoverSetup(m);
650:   GridDuplicateBC(mover->meshVelocityGrid,     m->meshVelocityGrid);
651:   GridDuplicateBC(mover->meshAccelerationGrid, m->meshAccelerationGrid);
652:   *newMover = m;
653:   return(0);
654: }

656: #undef  __FUNCT__
658: /*@
659:   MeshMoverView - Views a mesh mover object.

661:   Collective on MeshMover

663:   Input Parameters:
664: + mover  - The mover
665: - viewer - The viewer with which to view the mesh mover

667:   Level: beginner

669: .keywords: mesh, view, mover
670: .seealso: MeshMoverDestroy(), PetscViewerDrawOpen()
671: @*/
672: int MeshMoverView(MeshMover mover, PetscViewer viewer)
673: {
676:   if (viewer == PETSC_NULL) {
677:     viewer = PETSC_VIEWER_STDOUT_SELF;
678:   } else {
680:   }
681:   return(0);
682: }

684: #undef  __FUNCT__
686: /*@
687:   MeshMoverDestroy - Destroys a mesh mover object.

689:   Collective on MeshMover

691:   Input Parameter:
692: . mover - The mover

694:   Level: beginner

696: .keywords: mesh, destroy, mover
697: .seealso MeshMoverCreate()
698: @*/
699: int MeshMoverDestroy(MeshMover mover) {

704:   if (--mover->refct > 0) return(0);

706:   if (mover->nodeVelocities != PETSC_NULL) {
707:     GVecDestroy(mover->nodeVelocities);
708:     GVecDestroy(mover->nodeVelocitiesGhost);
709:   }
710:   if (mover->meshVelocitySles != PETSC_NULL) {
711:     GMatDestroy(mover->meshVelocityMat);
712:     GVecDestroy(mover->meshVelocityRhs);
713:     GVecDestroy(mover->meshVelocitySol);
714:     SLESDestroy(mover->meshVelocitySles);
715:   }
716:   if (mover->meshVelocityGrid != PETSC_NULL) {
717:     GridFinalizeBC(mover->meshVelocityGrid);
718:     GridDestroy(mover->meshVelocityGrid);
719:   }
720:   if (mover->nodeAccelerations != PETSC_NULL) {
721:     GVecDestroy(mover->nodeAccelerations);
722:     GVecDestroy(mover->nodeAccelerationsGhost);
723:   }
724:   if (mover->meshAccelerationSles != PETSC_NULL) {
725:     GMatDestroy(mover->meshAccelerationMat);
726:     GVecDestroy(mover->meshAccelerationRhs);
727:     GVecDestroy(mover->meshAccelerationSol);
728:     SLESDestroy(mover->meshAccelerationSles);
729:   }
730:   if (mover->meshAccelerationGrid != PETSC_NULL) {
731:     GridFinalizeBC(mover->meshAccelerationGrid);
732:     GridDestroy(mover->meshAccelerationGrid);
733:   }
734:   if (mover->meshVelocityScatter != PETSC_NULL) {
735:     VecScatterDestroy(mover->meshVelocityScatter);
736:   }
737:   if (mover->movingMeshViewer != PETSC_NULL) {
738:     PetscViewerDestroy(mover->movingMeshViewer);
739:   }

741:   (*mover->ops->destroy)(mover);

743:   PetscLogObjectDestroy(mover);
744:   PetscHeaderDestroy(mover);
745:   return(0);
746: }

748: #undef  __FUNCT__
750: /*@
751:   MeshMoverSetMesh - This function returns the associated Mesh.

753:   Not collective

755:   Input Parameter:
756: . mover - The mover

758:   Output Parameter:
759: . mesh  - The mesh, or PETSC_NULL if it does not exist

761:   Level: intermediate

763: .keywords: Mesh, movement, mover
764: .seealso: MeshSetMover()
765: @*/
766: int MeshMoverSetMesh(MeshMover mover, Mesh mesh) {
770:   mover->mesh = mesh;
771:   return(0);
772: }

774: #undef  __FUNCT__
776: /*@
777:   MeshMoverGetMesh - This function returns the associated Mesh.

779:   Not collective

781:   Input Parameter:
782: . mover - The mover

784:   Output Parameter:
785: . mesh  - The mesh, or PETSC_NULL if it does not exist

787:   Level: intermediate

789: .keywords: Mesh, movement, mover
790: .seealso: MeshGetMover()
791: @*/
792: int MeshMoverGetMesh(MeshMover mover, Mesh *mesh) {
796:   *mesh = mover->mesh;
797:   return(0);
798: }

800: #undef  __FUNCT__
802: /*@
803:   MeshMoverGetMovement - Returns the type of mesh movement calculation.

805:   Not collective

807:   Input Parameter:
808: . mover - The mover

810:   Output Parameters:
811: + vtype - The type of mesh velocity calculation, like MESH_LAPLACIAN
812: . atype - The type of mesh acceleration calculation, like MESH_LAPLACIAN
813: - ctx   - A user context

815:   Level: intermediate

817: .keywords: mesh, movement
818: .seealso: MeshMoverSetMovement()
819: @*/
820: int MeshMoverGetMovement(MeshMover mover, MeshSolveMethod *vtype, MeshSolveMethod *atype, PetscObject *ctx)
821: {
824:   if (vtype != PETSC_NULL) {
826:     *vtype = mover->meshVelocityMethod;
827:   }
828:   if (atype != PETSC_NULL) {
830:     *atype = mover->meshAccelerationMethod;
831:   }
832:   if (ctx   != PETSC_NULL) {
834:     *ctx   = mover->movingMeshCtx;
835:   }
836:   return(0);
837: }

839: #undef  __FUNCT__
841: /*@
842:   MeshMoverSetMovement - Determines the type of mesh mvoement calculation

844:   Collective on MeshMover

846:   Input Parameters:
847: + mover  - The mover
848: . vtype  - The type of mesh velocity calculation, like MESH_LAPLACIAN
849: . atype  - The type of mesh acceleration calculation, like MESH_LAPLACIAN
850: - ctx    - A user context

852:   Level: intermediate

854: .keywords: mesh, movement
855: .seealso: MeshMoverGetMovement()
856: @*/
857: int MeshMoverSetMovement(MeshMover mover, MeshSolveMethod vtype, MeshSolveMethod atype, PetscObject ctx)
858: {

863:   if (vtype != PETSC_DECIDE) mover->meshVelocityMethod     = vtype;
864:   if (atype != PETSC_DECIDE) mover->meshAccelerationMethod = atype;
865:   if (ctx   != PETSC_NULL)   mover->movingMeshCtx          = ctx;
866:   if (mover->meshVelocityGrid == PETSC_NULL) {
867:     MeshMoverSetup(mover);
868:   }
869:   return(0);
870: }

872: #undef  __FUNCT__
874: /*@
875:   MeshMoverGetVelocityGrid - Returns the grid which defines the mesh velocity.

877:   Not collective

879:   Input Parameter:
880: . mover    - The mover

882:   Output Parameter:
883: . velGrid - The mesh velocity grid

885:   Level: intermediate

887: .keywords: grid, mesh, mesh velocity, movement
888: .seealso: MeshMoverGetAccelerationGrid(), GridGetMeshMovement()
889: @*/
890: int MeshMoverGetVelocityGrid(MeshMover mover, Grid *velGrid)
891: {
894:   *velGrid = mover->meshVelocityGrid;
895:   return(0);
896: }

898: #undef  __FUNCT__
900: /*@
901:   MeshMoverGetAccelerationGrid - Returns the grid which defines the mesh acceleration.

903:   Not collective

905:   Input Parameter:
906: . mover    - The mover

908:   Output Parameter:
909: . accGrid - The mesh acceleration grid

911:   Level: intermediate

913: .keywords: grid, mesh, mesh acceleration, movement
914: .seealso: MeshMoverGetVelocityGrid(), GridGetMeshMovement()
915: @*/
916: int MeshMoverGetAccelerationGrid(MeshMover mover, Grid *accGrid)
917: {
920:   *accGrid = mover->meshAccelerationGrid;
921:   return(0);
922: }

924: #undef  __FUNCT__
926: /*@
927:   MeshMoverCalcVelocityBCValues - This function calculates the boundary values for the
928:   mesh velocity. It is normally called once a timestep when using time dependent boundary
929:   conditions.

931:   Collective on MeshMover

933:   Input Parameters:
934: + mover - The MeshMover
935: . save - A flag used to store old values, usually for timestepping
936: - ctx  - The context

938:   Level: advanced

940: .keywords: MeshMover, reduction, boundary conditions
941: .seealso: MeshMoverSetVelocityBCContext(), MeshMoverSetVelocityBC()
942: @*/
943: int MeshMoverCalcVelocityBCValues(MeshMover mover)
944: {
945:   Grid       g;
946:   PetscTruth reduce;
947:   void      *reduceCtx;
948:   int        ierr;

952:   MeshMoverGetMovement(mover, PETSC_NULL, PETSC_NULL, PETSC_NULL);
953:   MeshMoverGetVelocityGrid(mover, &g);
955:   GridGetReduceSystem(g, &reduce);
956:   if (reduce == PETSC_FALSE) return(0);
957:   GridSetUp(g);
958:   GridGetBCContext(g, &reduceCtx);
959:   GridCalcBCValues(g, PETSC_FALSE, reduceCtx);
960:   return(0);
961: }

963: #undef  __FUNCT__
965: /*@
966:   MeshMoverCalcAccelerationBCValues - This function calculates the boundary values for the
967:   mesh acceleration. It is normally called once a timestep when using time dependent boundary
968:   conditions.

970:   Collective on MeshMover

972:   Input Parameters:
973: + mover - The MeshMover
974: . save - A flag used to store old values, usually for timestepping
975: - ctx  - The context

977:   Level: advanced

979: .keywords: MeshMover, reduction, boundary conditions
980: .seealso: MeshMoverSetAccelerationBCContext(), MeshMoverSetAccelerationBC()
981: @*/
982: int MeshMoverCalcAccelerationBCValues(MeshMover mover)
983: {
984:   Grid       g;
985:   PetscTruth reduce;
986:   void      *reduceCtx;
987:   int        ierr;

991:   MeshMoverGetMovement(mover, PETSC_NULL, PETSC_NULL, PETSC_NULL);
992:   MeshMoverGetAccelerationGrid(mover, &g);
994:   GridGetReduceSystem(g, &reduce);
995:   if (reduce == PETSC_FALSE) return(0);
996:   GridSetUp(g);
997:   GridGetBCContext(g, &reduceCtx);
998:   GridCalcBCValues(g, PETSC_FALSE, reduceCtx);
999:   return(0);
1000: }

1002: #undef  __FUNCT__
1004: /*@
1005:   MeshMoverCalcNodeVelocities - Recalculates the node velocities.

1007:   Collective on MeshMover

1009:   Input Parameters:
1010: + mover - The mover
1011: - flag - Flag for reuse of old operator in calculation

1013:   Level: advanced

1015: .keywords: mesh, movement
1016: .seealso: MeshMoverCalcNodeAccelerations()
1017: @*/
1018: int MeshMoverCalcNodeVelocities(MeshMover mover, PetscTruth flag)
1019: {
1020:   Grid         grid          = mover->meshVelocityGrid;
1021:   MatStructure matStructFlag = DIFFERENT_NONZERO_PATTERN;
1022:   Mesh         mesh;
1023:   SLES         sles;
1024:   GMat         mat;
1025:   GVec         rhs, sol;
1026:   PetscScalar  zero = 0.0;
1027:   int          its;
1028:   int          ierr;


1033:   PetscLogEventBegin(MESH_CalcNodeVelocities, mover, 0, 0, 0);

1035:   /* Allocation */
1036:   if (mover->meshVelocitySles == PETSC_NULL) {
1037:     GVecCreate(grid, &mover->meshVelocityRhs);
1038:     PetscLogObjectParent(mover, mover->meshVelocityRhs);
1039:     GVecCreate(grid, &mover->meshVelocitySol);
1040:     PetscLogObjectParent(mover, mover->meshVelocitySol);
1041:     GMatCreate(grid, &mover->meshVelocityMat);
1042:     PetscLogObjectParent(mover, mover->meshVelocityMat);
1043:     SLESCreate(mover->comm, &mover->meshVelocitySles);
1044:     SLESAppendOptionsPrefix(mover->meshVelocitySles, "mesh_vel_");
1045:     SLESSetFromOptions(mover->meshVelocitySles);
1046:     PetscLogObjectParent(mover, mover->meshVelocitySles);
1047:   }
1048:   sles = mover->meshVelocitySles;
1049:   mat  = mover->meshVelocityMat;
1050:   rhs  = mover->meshVelocityRhs;
1051:   sol  = mover->meshVelocitySol;

1053:   /* Form System Matrix */
1054:   if (flag == PETSC_TRUE) {
1055:     MatZeroEntries(mat);
1056:     GridEvaluateSystemMatrix(grid, PETSC_NULL, &mat, &mat, &matStructFlag, mover->movingMeshCtx);
1057:     GMatSetBoundary(mat, 1.0, mover->movingMeshCtx);
1058:     MatSetOption(mat, MAT_NO_NEW_NONZERO_LOCATIONS);
1059:   }
1060:   /* Form Rhs */
1061:   VecSet(&zero, rhs);
1062:   GridEvaluateRhs(grid, PETSC_NULL, rhs, mover->movingMeshCtx);
1063:   GVecSetBoundary(rhs, mover->movingMeshCtx);
1064:   /* Solve system */
1065:   SLESSetOperators(sles, mat, mat, SAME_NONZERO_PATTERN);
1066:   SLESSolve(sles, rhs, sol, &its);
1067:   /* Copy values into the mesh vector */
1068:   MeshMoverGetMesh(mover, &mesh);
1069:   MeshUpdateNodeValues(mesh, sol, mover->nodeVelocities, mover->nodeVelocitiesGhost);
1070:   PetscLogEventEnd(MESH_CalcNodeVelocities, mover, 0, 0, 0);
1071:   return(0);
1072: }

1074: #undef  __FUNCT__
1076: /*@
1077:   MeshMoverCalcNodeAccelerations - Recalculates the node accelerations.

1079:   Collective on MeshMover

1081:   Input Parameters:
1082: + mover - The mover
1083: - flag - Flag for reuse of old operator in calculation

1085:   Level: advanced

1087: .keywords: mesh, movement
1088: .seealso: MeshMoverCalcNodeVelocities()
1089: @*/
1090: int MeshMoverCalcNodeAccelerations(MeshMover mover, PetscTruth flag)
1091: {
1092:   Grid         grid          = mover->meshAccelerationGrid;
1093:   MatStructure matStructFlag = DIFFERENT_NONZERO_PATTERN;
1094:   Mesh         mesh;
1095:   SLES         sles;
1096:   GMat         mat;
1097:   GVec         rhs, sol;
1098:   PetscScalar  zero = 0.0;
1099:   int          its;
1100:   int          ierr;


1105:   PetscLogEventBegin(MESH_CalcNodeAccelerations, mover, 0, 0, 0);

1107:   /* Allocation */
1108:   if (mover->meshAccelerationSles == PETSC_NULL) {
1109:     GVecCreate(grid, &mover->meshAccelerationRhs);
1110:     PetscLogObjectParent(mover, mover->meshAccelerationRhs);
1111:     GVecCreate(grid, &mover->meshAccelerationSol);
1112:     PetscLogObjectParent(mover, mover->meshAccelerationSol);
1113:     GMatCreate(grid, &mover->meshAccelerationMat);
1114:     PetscLogObjectParent(mover, mover->meshAccelerationMat);
1115:     SLESCreate(mover->comm, &mover->meshAccelerationSles);
1116:     SLESAppendOptionsPrefix(mover->meshAccelerationSles, "mesh_accel_");
1117:     SLESSetFromOptions(mover->meshAccelerationSles);
1118:     PetscLogObjectParent(mover, mover->meshAccelerationSles);
1119:   }
1120:   sles = mover->meshAccelerationSles;
1121:   mat  = mover->meshAccelerationMat;
1122:   rhs  = mover->meshAccelerationRhs;
1123:   sol  = mover->meshAccelerationSol;

1125:   /* Form System Matrix */
1126:   if (flag == PETSC_TRUE) {
1127:     MatZeroEntries(mat);
1128:     GridEvaluateSystemMatrix(grid, PETSC_NULL, &mat, &mat, &matStructFlag, mover->movingMeshCtx);
1129:     GMatSetBoundary(mat, 1.0, mover->movingMeshCtx);
1130:     MatSetOption(mat, MAT_NO_NEW_NONZERO_LOCATIONS);
1131:   }
1132:   /* Form Rhs */
1133:   VecSet(&zero, rhs);
1134:   GridEvaluateRhs(grid, PETSC_NULL, rhs, mover->movingMeshCtx);
1135:   GVecSetBoundary(rhs, mover->movingMeshCtx);
1136:   /* Solve system */
1137:   SLESSetOperators(sles, mat, mat, SAME_NONZERO_PATTERN);
1138:   SLESSolve(sles, rhs, sol, &its);
1139:   /* Copy values into the mesh vector */
1140:   MeshMoverGetMesh(mover, &mesh);
1141:   MeshUpdateNodeValues(mesh, sol, mover->nodeAccelerations, mover->nodeAccelerationsGhost);
1142:   PetscLogEventEnd(MESH_CalcNodeAccelerations, mover, 0, 0, 0);
1143:   return(0);
1144: }

1146: #undef  __FUNCT__
1148: /*@
1149:   MeshMoverSetNodeVelocities - Allows the user to specify node velocities.

1151:   Collective on MeshMover

1153:   Input Parameters:
1154: + mover  - The mover
1155: . func  - A PointFunction defining the node velocities
1156: . alpha - A scalar multiplier
1157: - ctx   - A user context

1159:   Level: advanced

1161: .keywords: mesh, movement
1162: .seealso: MeshMoverSetNodeAccelerations()
1163: @*/
1164: int MeshMoverSetNodeVelocities(MeshMover mover, PointFunction func, PetscScalar alpha, void *ctx)
1165: {
1166: #if 0
1167:   int field = 0;

1171:   GVecEvaluateFunction(mover->nodeVelocities, 1, &field, func, alpha, ctx);
1172:   /* Now need a way to get mesh velocities back to the system */
1173:   return(0);
1174: #else
1175:   SETERRQ(PETSC_ERR_SUP, "Not yet a grid vector");
1176: #endif
1177: }

1179: #undef  __FUNCT__
1181: /*@
1182:   MeshMoverSetNodeAccelerations - Allows the user to specify node accelerations.

1184:   Collective on MeshMover

1186:   Input Parameters:
1187: + mover  - The mover
1188: . func  - A function defining the node accelerations
1189: . alpha - A scalar multiplier
1190: - ctx   - A user context

1192:   Level: intermediate

1194: .keywords: mesh, movement
1195: .seealso: MeshMoverSetNodeVelocities()
1196: @*/
1197: int MeshMoverSetNodeAccelerations(MeshMover mover, PointFunction func, PetscScalar alpha, void *ctx)
1198: {
1199: #if 0
1200:   int field = 0;

1204:   GVecEvaluateFunction(mover->nodeAccelerations, 1, &field, func, alpha, ctx);
1205:   /* Now need a way to get mesh accelerations back to the system */
1206:   return(0);
1207: #else
1208:   SETERRQ(PETSC_ERR_SUP, "Not yet a grid vector");
1209: #endif
1210: }

1212: #undef  __FUNCT__
1214: /*@C
1215:   MeshMoverSetVelocityBC -        This function sets the boundary conditions to use
1216:   for the mesh velocity calculation.

1218:   Collective on MeshMover

1220:   Input Parameters:
1221: + mover   - The mover
1222: . bd     - The marker for the boundary along which conditions are applied
1223: . f      - The function which defines the boundary condition
1224: - reduce - The flag for explicit reduction of the system

1226:   Note:
1227:   This function clears any previous conditions, whereas MeshMoverAddVelocityBC()
1228:   appends the condition.

1230:   Level: intermediate

1232: .keywords mesh, velocity, boundary condition
1233: .seealso MeshMoverAddVelocityBC(), MeshMoverSetAccelerationBC(), GridSetBC(), GridAddBC()
1234: @*/
1235: int MeshMoverSetVelocityBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce)
1236: {


1242:   if (mover->meshVelocityGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Velocity grid not created");
1243:   GridSetBC(mover->meshVelocityGrid, bd, 0, f, reduce);
1244:   return(0);
1245: }

1247: #undef  __FUNCT__
1249: /*@C
1250:   MeshMoverAddVelocityBC - This function adds a boundary condition to use
1251:   for the mesh velocity calculation.

1253:   Collective on MeshMover

1255:   Input Parameters:
1256: + mover   - The mover
1257: . bd     - The marker for the boundary along which conditions are applied
1258: . f      - The function which defines the boundary condition
1259: - reduce - The flag for explicit reduction of the system

1261:   Note:
1262:   This function appends the condition, whereas MeshMoverSetVelocityBC()
1263:   clears any previous conditions.

1265:   Level: intermediate

1267: .keywords mesh, velocity, boundary condition
1268: .seealso MeshMoverSetVelocityBC(), MeshMoverAddAccelerationBC(), GridSetBC(), GridAddBC()
1269: @*/
1270: int MeshMoverAddVelocityBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce)
1271: {


1277:   if (mover->meshVelocityGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Velocity grid not created");
1278:   GridAddBC(mover->meshVelocityGrid, bd, 0, f, reduce);
1279:   return(0);
1280: }

1282: #undef  __FUNCT__
1284: /*@C
1285:   MeshMoverSetVelocityBCContext - This function sets the boundary condition context
1286:   to use for the mesh velocity calculation.

1288:   Collective on MeshMover

1290:   Input Parameters:
1291: + mover - The mover
1292: - ctx  - The user context

1294:   Level: intermediate

1296: .keywords mesh, velocity, boundary condition
1297: .seealso MeshMoverAddVelocityBC(), MeshMoverSetVelocityBC()
1298: @*/
1299: int MeshMoverSetVelocityBCContext(MeshMover mover, void *ctx)
1300: {


1306:   if (mover->meshVelocityGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Velocity grid not created");
1307:   GridSetBCContext(mover->meshVelocityGrid, ctx);
1308:   return(0);
1309: }

1311: #undef  __FUNCT__
1313: /*@C
1314:   MeshMoverSetAccelerationBC - This function sets the boundary condition to use
1315:   for the mesh acceleration calculation.

1317:   Collective on MeshMover

1319:   Input Parameters:
1320: + mover   - The mover
1321: . bd     - The marker for the boundary along which conditions are applied
1322: . f      - The function which defines the boundary condition
1323: - reduce - The flag for explicit reduction of the system

1325:   Note:
1326:   This function clears any previous conditions, whereas MeshMoverAddAccelerationBC()
1327:   appends the condition.

1329:   Level: intermediate

1331: .keywords mesh, velocity, boundary condition
1332: .seealso MeshMoverAddAccelerationBC(), MeshMoverSetVelocityBC(), GridSetBC(), GridAddBC()
1333: @*/
1334: int MeshMoverSetAccelerationBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce)
1335: {


1341:   if (mover->meshAccelerationGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Acceleration grid not created");
1342:   GridSetBC(mover->meshAccelerationGrid, bd, 0, f, reduce);
1343:   return(0);
1344: }

1346: #undef  __FUNCT__
1348: /*@C
1349:   MeshMoverAddAccelerationBC - This function adds a boundary condition to use
1350:   for the mesh acceleration calculation.

1352:   Collective on MeshMover

1354:   Input Parameters:
1355: + mover   - The mover
1356: . bd     - The marker for the boundary along which conditions are applied
1357: . f      - The function which defines the boundary condition
1358: - reduce - The flag for explicit reduction of the system

1360:   Note:
1361:   This function appends the condition, whereas MeshMoverSetVelocityBC()
1362:   clears any previous conditions.

1364:   Level: intermediate

1366: .keywords mesh, velocity, boundary condition
1367: .seealso MeshMoverSetAccelerationBC(), MeshMoverAddVelocityBC(), GridSetBC(), GridAddBC()
1368: @*/
1369: int MeshMoverAddAccelerationBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce)
1370: {


1376:   if (mover->meshAccelerationGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Acceleration grid not created");
1377:   GridAddBC(mover->meshAccelerationGrid, bd, 0, f, reduce);
1378:   return(0);
1379: }

1381: #undef  __FUNCT__
1383: /*@C
1384:   MeshMoverSetAccelerationBCContext - This function sets the boundary condition context
1385:   to use for the mesh acceleration calculation.

1387:   Collective on MeshMover

1389:   Input Parameter:
1390: + mover - The mover
1391: - ctx  - The user context

1393:   Level: intermediate

1395: .keywords mesh, velocity, boundary condition
1396: .seealso MeshMoverAddAccelerationBC(), MeshMoverSetAccelerationBC()
1397: @*/
1398: int MeshMoverSetAccelerationBCContext(MeshMover mover, void *ctx)
1399: {


1405:   if (mover->meshAccelerationGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Acceleration grid not created");
1406:   GridSetBCContext(mover->meshAccelerationGrid, ctx);
1407:   return(0);
1408: }

1410: #undef  __FUNCT__
1412: /*@
1413:   MeshMoverSetMovementCaption - Sets the caption of the moving mesh
1414:   viewer. Can be used to provide auxilliary information, e.g. time.

1416:   Collective on MeshMover

1418:   Input Parameters:
1419: + mover    - The mover
1420: - caption - The caption

1422:   Level: intermediate

1424: .keywords: mesh, movement
1425: .seealso: MeshMoverMoveMesh()
1426: @*/
1427: int MeshMoverSetMovementCaption(MeshMover mover, char *caption)
1428: {
1432:   mover->movingMeshViewerCaption = caption;
1433:   return(0);
1434: }

1436: /*----------------------------------------- MeshMover Registration Functions -----------------------------------------*/
1437: PetscFList MeshMoverList              = 0;
1438: int        MeshMoverRegisterAllCalled = 0;

1440: #undef __FUNCT__  
1442: /*@C
1443:   MeshMoverSetType - Builds a MeshMover, for a particular Mesh implementation.

1445:   Collective on MeshMover

1447:   Input Parameters:
1448: + mover - The MeshMover
1449: - type  - The type name

1451:   Options Database Key:
1452: . -mesh_mover_type <type> - Sets the mesh mvoer type

1454:   Level: intermediate

1456: .keywords: vector, set, type
1457: .seealso: MeshMoverCreate()
1458: @*/
1459: int MeshMoverSetType(MeshMover mover, MeshMoverType method) {
1460:   int      (*r)(MeshMover);
1461:   PetscTruth match;
1462:   int        ierr;

1466:   PetscTypeCompare((PetscObject) mover, method, &match);
1467:   if (match == PETSC_TRUE) return(0);

1469:   /* Get the function pointers for the vector requested */
1470:   if (!MeshMoverRegisterAllCalled) {
1471:     MeshMoverRegisterAll(PETSC_NULL);
1472:   }
1473:   PetscFListFind(mover->comm, MeshMoverList, method,(void (**)(void)) &r);
1474:   if (!r) SETERRQ1(PETSC_ERR_ARG_WRONG, "Unknown mesh mover type: %s", method);

1476:   if (mover->ops->destroy != PETSC_NULL) {
1477:     (*mover->ops->destroy)(mover);
1478:   }
1479:   (*r)(mover);

1481:   PetscObjectChangeTypeName((PetscObject) mover, method);
1482:   return(0);
1483: }

1485: #undef __FUNCT__
1487: /*@C
1488:   MeshMoverGetType - Gets the MeshMover method type name (as a string).

1490:   Not collective

1492:   Input Parameter:
1493: . mover - The mover

1495:   Output Parameter:
1496: . type - The name of MeshMover method

1498:   Level: intermediate

1500: .keywords: mover, get, type, name
1501: .seealso MeshMoverSetType()
1502: @*/
1503: int MeshMoverGetType(MeshMover mover, MeshMoverType *type)
1504: {

1510:   if (!MeshMoverRegisterAllCalled) {
1511:     MeshMoverRegisterAll(PETSC_NULL);
1512:   }
1513:   *type = mover->type_name;
1514:   return(0);
1515: }

1517: /*@C
1518:   MeshMoverRegister - Adds a creation method to the mesh movement package.

1520:   Synopsis:

1522:   MeshMoverRegister(char *name, char *path, char *func_name, int (*create_func)(MeshMover))

1524:   Not Collective

1526:   Input Parameters:
1527: + name        - The name of a new user-defined creation routine
1528: . path        - The path (either absolute or relative) of the library containing this routine
1529: . func_name   - The name of the creation routine
1530: - create_func - The creation routine itself

1532:   Notes:
1533:   MeshMoverRegister() may be called multiple times to add several user-defined meshes.

1535:   If dynamic libraries are used, then the fourth input argument (create_func) is ignored.

1537:   Sample usage:
1538: .vb
1539:   MeshMoverRegisterDynamic("my_mesh", "/home/username/my_lib/lib/libO/solaris/libmy.a", "MyMeshMoverCreate", MyMeshMoverCreate);
1540: .ve

1542:   Then, your mesh type can be chosen with the procedural interface via
1543: .vb
1544:     MeshMoverCreate(MPI_Comm, MeshMover *);
1545:     MeshMoverSetType(vec, "my_mesh_mover")
1546: .ve
1547:   or at runtime via the option
1548: .vb
1549:     -mesh_mover_type my_mesh_mover
1550: .ve

1552:   Note: $PETSC_ARCH and $BOPT occuring in pathname will be replaced with appropriate values.

1554:   Level: advanced

1556: .keywords: MeshMover, register
1557: .seealso: MeshMoverRegisterAll(), MeshMoverRegisterDestroy()
1558: @*/
1559: #undef __FUNCT__  
1561: int MeshMoverRegister(const char sname[], const char path[], const char name[], int (*function)(MeshMover))
1562: {
1563:   char fullname[256];
1564:   int  ierr;

1567:   PetscStrcpy(fullname, path);
1568:   PetscStrcat(fullname, ":");
1569:   PetscStrcat(fullname, name);
1570:   PetscFListAdd(&MeshMoverList, sname, fullname, (void (*)(void)) function);
1571:   return(0);
1572: }

1574: #undef __FUNCT__  
1576: /*@C
1577:   MeshMoverRegisterDestroy - Frees the list of creation routines for
1578:   meshes that were registered by MeshMoverRegister().

1580:   Not collective

1582:   Level: advanced

1584: .keywords: mesh, register, destroy
1585: .seealso: MeshMoverRegister(), MeshMoverRegisterAll()
1586: @*/
1587: int MeshMoverRegisterDestroy() {

1591:   if (MeshMoverList != PETSC_NULL) {
1592:     PetscFListDestroy(&MeshMoverList);
1593:     MeshMoverList = PETSC_NULL;
1594:   }
1595:   MeshMoverRegisterAllCalled = 0;
1596:   return(0);
1597: }

1599: EXTERN_C_BEGIN
1600: extern int MeshMoverCreate_Triangular_2D(MeshMover);
1601: EXTERN_C_END

1603: #undef __FUNCT__
1605: /*@C
1606:   MeshMoverRegisterAll - Registers all of the generation routines in the MeshMover package. 

1608:   Not Collective

1610:   Input parameter:
1611: . path - The dynamic library path

1613:   Level: advanced

1615: .keywords: MeshMover, register, all
1616: .seealso: MeshMoverRegister(), MeshMoverRegisterDestroy()
1617: @*/
1618: int MeshMoverRegisterAll(const char path[])
1619: {

1623:   MeshMoverRegisterAllCalled = 1;

1625:   MeshMoverRegister(MESH_MOVER_TRIANGULAR_2D, path, "MeshMoverCreate_Triangular_2D", MeshMoverCreate_Triangular_2D);
1626: 
1627:   return(0);
1628: }

1630: /*-------------------------------------------------- Grid Functions --------------------------------------------------*/
1631: #undef  __FUNCT__
1633: /*@
1634:   GridReformMesh - Reconstruct the grid from the current description
1635:   and underlying mesh boundary. The old mesh is returned if requested.
1636:   This allows postprocessing of the mesh before GridReform() is called.

1638:   Collective on Grid

1640:   Input Parameters:
1641: + grid   - The grid
1642: . refine - Flag indicating whether to refine or recreate the mesh
1643: . area   - [Optional] A function which gives area constraints for refinement after reform, PETSC_NULL for no refinement
1644: - newBd  - Flag to determine whether the boundary should be generated (PETSC_TRUE) or taken from storage

1646:   Output Parameter:
1647: . oldMesh - The old mesh

1649:   Level: advanced

1651: .keywords: Grid, mesh, reform, refine
1652: .seealso: MeshReform(), GridReform(), GridRefineMesh()
1653: @*/
1654: int GridReformMesh(Grid grid, PetscTruth refine, PointFunction area, PetscTruth newBd, Mesh *oldMesh)
1655: {
1656:   Mesh mesh;
1657:   Mesh newMesh;
1658:   int  ierr;


1663:   /* Reform mesh */
1664:   GridGetMesh(grid, &mesh);
1665:   MeshReform(mesh, refine, area, newBd, &newMesh);

1667:   /* Save old information */
1668:   if (oldMesh != PETSC_NULL){
1670:     *oldMesh = mesh;
1671:   }
1672:   grid->mesh = newMesh;

1674:   /* Recalculate grid structures for mesh */
1675:   GridCalcPointBCNodes(grid);

1677:   /* Recreate class structure */
1678:   FieldClassMapDestroy(grid->cm);
1679:   FieldClassMapCreateTriangular2D(grid, grid->numActiveFields, grid->defaultFields, &grid->cm);

1681:   /* Specific setup */
1682:   if (grid->ops->gridreformmesh != PETSC_NULL) {
1683:     (*grid->ops->gridreformmesh)(grid);
1684:   }
1685:   return(0);
1686: }

1688: #undef  __FUNCT__
1690: /*@
1691:   GridReform - Reconstruct the variable structures, usually after a change
1692:   in mesh structure. The optional arguments allow projection of the
1693:   previous solution onto the new mesh.

1695:   Collective on Grid

1697:   Input Parameters:
1698: + grid    - The grid
1699: . oldMesh - [Optional] The old mesh
1700: - x       - [Optional] A grid vector to project onto the new mesh

1702:   Output Parameter:
1703: . y       - [Optional] The projected grid vector

1705:   Level: advanced

1707: .keywords: grid, mesh, reform
1708: .seealso: GridReformMesh(), MeshReform()
1709: @*/
1710: int GridReform(Grid grid, Mesh oldMesh, GVec x, GVec *y)
1711: {
1712:   MPI_Comm             comm;
1713:   VarOrdering          order;
1714:   FieldClassMap        map;
1715:   int                  numFields;
1716:   int                 *fields;
1717:   PetscTruth           isConstrained, explicitConstraints;
1718:   InterpolationContext iCtx;             /* Used to form the rhs */
1719:   /* L_2 projection variables */
1720: #if 0
1721:   GMat                 M;                /* The (local) mass matrix M */
1722:   GMat                 m;                /* The restriction of M to the boundary m */
1723:   GMat                 S;                /* The Schur complement S = m^T M^{-1} m */
1724:   GVec                 sol, rhs;         /* The solution x and rhs f */
1725:   GVec                 projRhs;          /* The projection rhs g */
1726:   GVec                 uzawaRhs;         /* The Uzawa rhs m^T M^{-1} f - g */
1727:   GVec                 uzawaSol;         /* The Uzawa solution lambda = (m^T M^{-1} m)^{-1} (m^T M^{-1} f - g) */
1728:   KSP                  ksp;              /* The Krylov solver for M */
1729:   PC                   pc;               /* The preconditioner for M */
1730:   SLES                 sles;             /* The solver for M */
1731:   SLES                 uzawaSles;        /* The solver for S */
1732:   PetscScalar          minusOne = -1.0;
1733:   PetscScalar          zero     = 0.0;
1734:   int                  its;
1735: #endif
1736:   /* Local Interpolation variables */
1737:   VarOrdering          reduceOrder,     reduceOrderOld;      /* TEMPORARY, until I fix the interpolation */
1738:   GVec                 bdReduceVec,     prevBdReduceVec;     /* TEMPORARY, until I fix the interpolation */
1739:   GVec                 bdReduceVecOld,  prevBdReduceVecOld;  /* TEMPORARY, until I fix the interpolation */
1740:   GVec                 bdReduceVecDiff, prevBdReduceVecDiff; /* TEMPORARY, until I fix the interpolation */
1741:   PetscReal            norm;
1742:   int                  maxComp;
1743:   int                  field, f;
1744:   int                  ierr;

1748:   PetscLogEventBegin(GRID_Reform, grid, 0, 0, 0);
1749:   /* Make sure the vector is saved in local storage */
1750:   GridGlobalToLocal(grid, INSERT_VALUES, x);
1751:   /* Need to setup grid again */
1752:   grid->setupcalled   = PETSC_FALSE;
1753:   grid->bdSetupCalled = PETSC_FALSE;
1754:   /* Check arguments */
1755:   if ((x == PETSC_NULL) || (y == PETSC_NULL)) {
1756:     GridSetupGhostScatter(grid);
1757:     return(0);
1758:   }

1762:   PetscObjectGetComm((PetscObject) oldMesh, &comm);
1763:   GridIsConstrained(grid, &isConstrained);
1764:   GridGetExplicitConstraints(grid, &explicitConstraints);
1765:   if (explicitConstraints == PETSC_FALSE) {
1766:     order = grid->order;
1767:   } else {
1768:     order = grid->constraintOrder;
1769:   }
1770:   VarOrderingGetClassMap(order, &map);
1771: #ifdef INTERFACE_UPDATE
1772:   FieldClassMapGetNumFields(map, &numFields);
1773:   FieldClassMapGetFields(map, &fields);
1774: #endif
1775:   numFields = map->numFields;
1776:   fields    = map->fields;

1778:   /* Project fields onto new mesh */
1779:   switch(grid->interpolationType)
1780:   {
1781: #if 0
1782:   case INTERPOLATION_L2:
1783:     /* Interpolation Procedure:
1784:          We must locally interpolate the value of the old field at each Gaussian intergration
1785:          point in an element of the new field. The function PointFunctionInterpolateField()
1786:          calls GVecInterpolateField(), which calls DiscretizationInterpolateField(). The vector, field,
1787:          and old mesh are stored in a context so that PointFunctionInterpolateField() can pass
1788:          them to GVecInterpolateField(). The indices and values for DiscretizationInterpolateField() are
1789:          calculated by GridCalcLocalElementVecIndices() and GridLocalToElement(). These
1790:          functions need the old mesh, classes, variable ordering, and local vector ghostVec in
1791:          order to correctly calculate the old vector indices and values. However, the
1792:          GVecEvaluateFunctionGalerkin() uses the new values to get the weak form. The old mesh
1793:          we have, and the classes are stored in classesOld[]. We will keep the variable
1794:          ordering and local vector in temporary veriables. The old copies must be set to NULL
1795:          or they will be deleted.
1796:     */

1798:     /* Save old structures */
1799:     iCtx.vec      = x;
1800:     iCtx.mesh     = oldMesh;
1801:     iCtx.order    = order;
1802:     iCtx.ghostVec = grid->ghostVec;

1804:     /* Setup context */
1805:     for(f = 0, maxComp = 1; f < numFields; f++) {
1806:       maxComp = PetscMax(maxComp, grid->comp[fields[f]]);
1807:     }
1808:     MPI_Comm_size(comm, &iCtx.numProcs);
1809:     MPI_Comm_rank(comm, &iCtx.rank);
1810:     iCtx.batchMode = PETSC_FALSE;
1811:     iCtx.curPoint  = 0;
1812:     iCtx.maxPoints = iCtx.numProcs;
1813:     PetscMalloc((iCtx.numProcs+1)     * sizeof(int),         &iCtx.numPoints);
1814:     PetscMalloc(iCtx.numProcs         * sizeof(int),         &iCtx.activeProcs);
1815:     PetscMalloc(iCtx.numProcs         * sizeof(int),         &iCtx.calcProcs);
1816:     PetscMalloc(iCtx.maxPoints*3      * sizeof(PetscScalar), &iCtx.coords);
1817:     PetscMalloc(iCtx.numProcs*maxComp * sizeof(PetscScalar), &iCtx.values);

1819:     /* Create new variable numbering and ghost structures */
1820:     if (explicitConstraints == PETSC_TRUE) {
1821:       grid->constraintOrder = PETSC_NULL;
1822:     } else {
1823:       grid->order         = PETSC_NULL;
1824:     }
1825:     grid->ghostVec        = PETSC_NULL;
1826:     reduceOrderOld        = grid->reduceOrder;
1827:     prevBdReduceVec       = grid->bdReduceVec;
1828:     prevBdReduceVecOld    = grid->bdReduceVecOld;
1829:     prevBdReduceVecDiff   = grid->bdReduceVecDiff;
1830:     grid->reduceOrder     = PETSC_NULL;
1831:     grid->bdReduceVec     = PETSC_NULL;
1832:     grid->bdReduceVecOld  = PETSC_NULL;
1833:     grid->bdReduceVecDiff = PETSC_NULL;
1834:     GridSetupGhostScatter(grid);
1835:     reduceOrder           = grid->reduceOrder;
1836:     bdReduceVec           = grid->bdReduceVec;
1837:     bdReduceVecOld        = grid->bdReduceVecOld;
1838:     bdReduceVecDiff       = grid->bdReduceVecDiff;
1839:     grid->reduceOrder     = reduceOrderOld;
1840:     grid->bdReduceVec     = prevBdReduceVec;
1841:     grid->bdReduceVecOld  = prevBdReduceVecOld;
1842:     grid->bdReduceVecDiff = prevBdReduceVecDiff;
1843:     if (grid->bdReduceVecCur == bdReduceVec) {
1844:       grid->bdReduceVecCur = grid->bdReduceVec;
1845:     } else if (grid->bdReduceVecCur == bdReduceVecOld) {
1846:       grid->bdReduceVecCur = grid->bdReduceVecOld;
1847:     } else if (grid->bdReduceVecCur == bdReduceVecDiff) {
1848:       grid->bdReduceVecCur = grid->bdReduceVecDiff;
1849:     } else {
1850:       SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector");
1851:     }

1853:     /* Allocate structures */
1854:     if (explicitConstraints) {
1855:       GVecCreateConstrained(grid, y);
1856:     } else {
1857:       GVecCreate(grid, y);
1858:     }
1859:     GVecDuplicate(*y, &rhs);
1860:     GVecDuplicate(*y, &sol);
1861:     GVecCreateBoundaryRestriction(grid, &uzawaSol);
1862:     GVecDuplicate(uzawaSol, &projRhs);
1863:     GVecDuplicate(uzawaSol, &uzawaRhs);

1865:     /* Create the rhs */
1866:     for(f = 0; f < numFields; f++) {
1867:       field = fields[f];
1868:       iCtx.field = field;
1869:       if (iCtx.batchMode == PETSC_FALSE) {
1870:         GVecEvaluateFunctionGalerkinCollective(rhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
1871: 
1872: #if 1
1873:         GVecEvaluateBoundaryFunctionGalerkinCollective(projRhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
1874: 
1875: #else
1876:         VecSet(&zero, projRhs);
1877: #endif
1878:       } else {
1879:         GVecEvaluateFunctionGalerkin(rhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
1880:         GVecInterpolateFieldBatch(rhs, field, &iCtx);
1881:         iCtx.curPoint = iCtx.numPoints[iCtx.rank];
1882:         GVecEvaluateFunctionGalerkin(rhs, 1, &field, PointFunctionInterpolateFieldBatch, 1.0, &iCtx);
1883: 
1884:         if (iCtx.curPoint != iCtx.numPoints[iCtx.rank+1]) SETERRQ(PETSC_ERR_PLIB, "Invalid batch point processing");
1885:         GVecEvaluateBoundaryFunctionGalerkin(projRhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
1886: 
1887:         GVecInterpolateFieldBatch(projRhs, field, &iCtx);
1888:         iCtx.curPoint = iCtx.numPoints[iCtx.rank];
1889:         GVecEvaluateFunctionGalerkin(projRhs, 1, &field, PointFunctionInterpolateFieldBatch, 1.0, &iCtx);
1890: 
1891:         if (iCtx.curPoint != iCtx.numPoints[iCtx.rank+1]) SETERRQ(PETSC_ERR_PLIB, "Invalid batch point processing");
1892:       }
1893:     }
1894:     if (grid->bdReduceVecCur == grid->bdReduceVec) {
1895:       grid->bdReduceVecCur = bdReduceVec;
1896:     } else if (grid->bdReduceVecCur == grid->bdReduceVecOld) {
1897:       grid->bdReduceVecCur = bdReduceVecOld;
1898:     } else if (grid->bdReduceVecCur == grid->bdReduceVecDiff) {
1899:       grid->bdReduceVecCur = bdReduceVecDiff;
1900:     } else {
1901:       SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector");
1902:     }
1903:     grid->reduceOrder     = reduceOrder;
1904:     grid->bdReduceVec     = bdReduceVec;
1905:     grid->bdReduceVecOld  = bdReduceVecOld;
1906:     grid->bdReduceVecDiff = bdReduceVecDiff;
1907:     if (grid->reduceSystem == PETSC_TRUE) {
1908:       VarOrderingDestroy(reduceOrderOld);
1909:       VecDestroy(prevBdReduceVec);
1910:       VecDestroy(prevBdReduceVecOld);
1911:       VecDestroy(prevBdReduceVecDiff);
1912:     }

1914:     /* Destroy old structures */
1915:     VarOrderingDestroy(iCtx.order);
1916:     VecDestroy(iCtx.ghostVec);

1918:     /* Cleanup */
1919:     PetscFree(iCtx.numPoints);
1920:     PetscFree(iCtx.activeProcs);
1921:     PetscFree(iCtx.calcProcs);
1922:     PetscFree(iCtx.values);
1923:     PetscFree(iCtx.coords);

1925:     /* Create the inner solver */
1926:     SLESCreate(PETSC_COMM_SELF, &sles);
1927:     SLESAppendOptionsPrefix(sles, "grid_proj_");
1928:     SLESGetKSP(sles, &ksp);
1929:     SLESGetPC(sles, &pc);
1930:     KSPSetType(ksp, KSPCG);
1931:     PCSetType(pc, PCJACOBI);
1932:     SLESSetFromOptions(sles);
1933:     /* Create the outer solver */
1934:     SLESCreate(PETSC_COMM_SELF, &uzawaSles);
1935:     SLESAppendOptionsPrefix(uzawaSles, "grid_proj_uzawa_");
1936:     SLESGetKSP(uzawaSles, &ksp);
1937:     SLESGetPC(uzawaSles, &pc);
1938:     KSPSetType(ksp, KSPGMRES);
1939:     PCSetType(pc, PCNONE);
1940:     SLESSetFromOptions(uzawaSles);

1942:     /* Create M, the mass matrix */
1943:     GMatCreate(grid, &M);
1944:     GMatEvaluateOperatorGalerkin(M, PETSC_NULL, numFields, fields, fields, IDENTITY, 1.0, MAT_FINAL_ASSEMBLY, PETSC_NULL);
1945: 
1946:     MatNorm(M, NORM_FROBENIUS, &norm);
1947:     PetscPrintf(PETSC_COMM_SELF, "Mass matrix norm %lfn", norm);

1949:     /* Create m, the restriction of M to the boundary */
1950:     GMatCreateBoundaryRestriction(grid, &m);
1951: #if 1
1952:     GMatEvaluateBoundaryOperatorGalerkin(m, numFields, grid->cm->fields, grid->bdOrder, grid->bdLocOrder,
1953:                                                 fields, grid->order, grid->locOrder, IDENTITY, 1.0, MAT_FINAL_ASSEMBLY,
1954:                                                 PETSC_NULL);
1955: 
1956: #else
1957:     MatAssemblyBegin(m, MAT_FINAL_ASSEMBLY);
1958:     MatAssemblyEnd(m, MAT_FINAL_ASSEMBLY);
1959:     MatZeroEntries(m);
1960: #endif
1961:     MatNorm(m, NORM_FROBENIUS, &norm);
1962:     PetscPrintf(PETSC_COMM_SELF, "Restricted mass matrix norm %lfn", norm);

1964:     VecNorm(rhs,     NORM_2, &norm);
1965:     PetscPrintf(PETSC_COMM_SELF, "Rhs norm %lfn", norm);
1966:     VecNorm(projRhs, NORM_2, &norm);
1967:     PetscPrintf(PETSC_COMM_SELF, "Restricted rhs norm %lfn", norm);

1969:     /* Create the Schur complement */
1970:     SLESSetOperators(sles, M, M, SAME_NONZERO_PATTERN);
1971:     GMatCreateUzawa(M, m, sles, &S);

1973:     /* Create the restriction of the rhs */
1974:     SLESSolve(sles, rhs, sol, &its);
1975:     MatMultTranspose(m, sol, uzawaRhs);
1976:     VecAXPY(&minusOne, projRhs, uzawaRhs);
1977:     VecNorm(uzawaRhs, NORM_2, &norm);
1978:     PetscPrintf(PETSC_COMM_SELF, "Uzawa rhs norm %lfn", norm);

1980:     /* Solve system */
1981:     SLESSetOperators(uzawaSles, S, S, SAME_NONZERO_PATTERN);
1982:     SLESSolve(uzawaSles, uzawaRhs, uzawaSol, &its);
1983:     PetscPrintf(PETSC_COMM_SELF, "Uzawa solution took %d iterationsn", its);
1984:     VecNorm(uzawaSol, NORM_2, &norm);
1985:     PetscPrintf(PETSC_COMM_SELF, "Uzawa solution norm %lfn", norm);

1987:     /* Get solution */
1988:     MatMult(m, uzawaSol, sol);
1989:     VecNorm(sol, NORM_2, &norm);
1990:     PetscPrintf(PETSC_COMM_SELF, "Boundary solution update norm %lfn", norm);
1991:     VecAYPX(&minusOne, rhs, sol);
1992:     VecNorm(sol, NORM_2, &norm);
1993:     PetscPrintf(PETSC_COMM_SELF, "Solution update norm %lfn", norm);
1994:     SLESSolve(sles, sol, *y, &its);
1995:     PetscPrintf(PETSC_COMM_SELF, "Projection solution took %d iterationsn", its);

1997:     VecNorm(*y, NORM_2, &norm);
1998:     PetscPrintf(PETSC_COMM_SELF, "Projected solution norm %lfn", norm);

2000:     /* Cleanup */
2001:     SLESDestroy(sles);
2002:     SLESDestroy(uzawaSles);
2003:     GMatDestroyUzawa(S);
2004:     GMatDestroy(M);
2005:     GMatDestroy(m);
2006:     GVecDestroy(rhs);
2007:     GVecDestroy(sol);
2008:     GVecDestroy(projRhs);
2009:     GVecDestroy(uzawaRhs);
2010:     GVecDestroy(uzawaSol);

2012:     break;
2013: #endif
2014:   case INTERPOLATION_LOCAL:
2015:     /* Interpolation Procedure:
2016:          We must locally interpolate the value of the old field at each node in an element of
2017:          the new field. The function PointFunctionInterpolateField()calls GVecInterpolateField(),
2018:          which calls DiscretizationInterpolateField(). The vector, field, and old mesh are stored in a
2019:          context so that PointFunctionInterpolateField() can pass them to GVecInterpolateField().
2020:          The indices and values for DiscretizationInterpolateField() are calculated by
2021:          GridCalcLocalElementVecIndices() and GridLocalToElement(). These functions need the old
2022:          mesh, classes, variable ordering, and local vector ghostVec in order to correctly calculate
2023:          the old vector indices and values. The old mesh we have, and the classes are stored in
2024:          classesOld[]. We will keep the variable ordering and local vector in temporary veriables.
2025:          The old copies must be set to NULL or they will be deleted.
2026:     */

2028:     /* Save old structures */
2029:     iCtx.vec      = x;
2030:     iCtx.mesh     = oldMesh;
2031:     iCtx.order    = order;
2032:     iCtx.ghostVec = grid->ghostVec;

2034:     /* Setup context */
2035:     for(f = 0, maxComp = 1; f < numFields; f++) {
2036:       maxComp = PetscMax(maxComp, grid->fields[fields[f]].numComp);
2037:     }
2038:     MPI_Comm_size(comm, &iCtx.numProcs);
2039:     MPI_Comm_rank(comm, &iCtx.rank);
2040:     iCtx.batchMode   = PETSC_FALSE;
2041:     iCtx.curPoint    = 0;
2042:     iCtx.maxPoints   = iCtx.numProcs;
2043:     PetscMalloc((iCtx.numProcs+1)     * sizeof(int),         &iCtx.numPoints);
2044:     PetscMalloc(iCtx.numProcs         * sizeof(int),         &iCtx.activeProcs);
2045:     PetscMalloc(iCtx.numProcs         * sizeof(int),         &iCtx.calcProcs);
2046:     PetscMalloc(iCtx.maxPoints*3      * sizeof(PetscScalar), &iCtx.coords);
2047:     PetscMalloc(iCtx.numProcs*maxComp * sizeof(PetscScalar), &iCtx.values);
2048:     PetscMemzero(iCtx.numPoints, (iCtx.numProcs+1) * sizeof(int));
2049:     PetscOptionsHasName(PETSC_NULL, "-grid_int_mode", &iCtx.batchMode);

2051:     /* Create new variable numbering and ghost structures */
2052:     if (explicitConstraints == PETSC_TRUE) {
2053:       grid->constraintOrder = PETSC_NULL;
2054:     } else {
2055:       grid->order         = PETSC_NULL;
2056:     }
2057:     grid->ghostVec        = PETSC_NULL;
2058:     reduceOrderOld        = grid->reduceOrder;
2059:     prevBdReduceVec       = grid->bdReduceVec;
2060:     prevBdReduceVecOld    = grid->bdReduceVecOld;
2061:     prevBdReduceVecDiff   = grid->bdReduceVecDiff;
2062:     grid->reduceOrder     = PETSC_NULL;
2063:     grid->bdReduceVec     = PETSC_NULL;
2064:     grid->bdReduceVecOld  = PETSC_NULL;
2065:     grid->bdReduceVecDiff = PETSC_NULL;
2066:     GridSetupGhostScatter(grid);
2067:     reduceOrder           = grid->reduceOrder;
2068:     bdReduceVec           = grid->bdReduceVec;
2069:     bdReduceVecOld        = grid->bdReduceVecOld;
2070:     bdReduceVecDiff       = grid->bdReduceVecDiff;
2071:     grid->reduceOrder     = reduceOrderOld;
2072:     grid->bdReduceVec     = prevBdReduceVec;
2073:     grid->bdReduceVecOld  = prevBdReduceVecOld;
2074:     grid->bdReduceVecDiff = prevBdReduceVecDiff;
2075:     if (grid->bdReduceVecCur == bdReduceVec) {
2076:       grid->bdReduceVecCur = grid->bdReduceVec;
2077:     } else if (grid->bdReduceVecCur == bdReduceVecOld) {
2078:       grid->bdReduceVecCur = grid->bdReduceVecOld;
2079:     } else if (grid->bdReduceVecCur == bdReduceVecDiff) {
2080:       grid->bdReduceVecCur = grid->bdReduceVecDiff;
2081:     } else {
2082:       SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector");
2083:     }

2085:     /* Allocate structures */
2086:     if (explicitConstraints) {
2087:       GVecCreateConstrained(grid, y);
2088:     } else {
2089:       GVecCreate(grid, y);
2090:     }

2092:     /* Locally interpolate the solution values */
2093:     for(f = 0; f < numFields; f++) {
2094:       field = fields[f];
2095:       iCtx.field = field;
2096:       if (iCtx.batchMode == PETSC_FALSE) {
2097:         GVecEvaluateFunctionCollective(*y, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
2098:       } else {
2099:         GVecEvaluateFunction(*y, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx);
2100:         GVecInterpolateFieldBatch(*y, field, &iCtx);
2101:         iCtx.curPoint = iCtx.numPoints[iCtx.rank];
2102:         GVecEvaluateFunction(*y, 1, &field, PointFunctionInterpolateFieldBatch, 1.0, &iCtx);
2103:         if (iCtx.curPoint != iCtx.numPoints[iCtx.rank+1])
2104:           SETERRQ(PETSC_ERR_PLIB, "Invalid batch point processing");
2105:       }
2106:     }
2107:     if (grid->bdReduceVecCur == grid->bdReduceVec) {
2108:       grid->bdReduceVecCur = bdReduceVec;
2109:     } else if (grid->bdReduceVecCur == grid->bdReduceVecOld) {
2110:       grid->bdReduceVecCur = bdReduceVecOld;
2111:     } else if (grid->bdReduceVecCur == grid->bdReduceVecDiff) {
2112:       grid->bdReduceVecCur = bdReduceVecDiff;
2113:     } else {
2114:       SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector");
2115:     }
2116:     grid->reduceOrder     = reduceOrder;
2117:     grid->bdReduceVec     = bdReduceVec;
2118:     grid->bdReduceVecOld  = bdReduceVecOld;
2119:     grid->bdReduceVecDiff = bdReduceVecDiff;
2120:     if (grid->reduceSystem == PETSC_TRUE) {
2121:       VarOrderingDestroy(reduceOrderOld);
2122:       VecDestroy(prevBdReduceVec);
2123:       VecDestroy(prevBdReduceVecOld);
2124:       VecDestroy(prevBdReduceVecDiff);
2125:     }

2127:     VecNorm(*y, NORM_2, &norm);
2128:     PetscPrintf(PETSC_COMM_SELF, "Projected solution norm %gn", norm);

2130:     /* Destroy old structures */
2131:     VarOrderingDestroy(iCtx.order);
2132:     VecDestroy(iCtx.ghostVec);

2134:     /* Cleanup */
2135:     PetscFree(iCtx.numPoints);
2136:     PetscFree(iCtx.activeProcs);
2137:     PetscFree(iCtx.calcProcs);
2138:     PetscFree(iCtx.values);
2139:     PetscFree(iCtx.coords);
2140:     break;
2141:   case INTERPOLATION_HALO:
2142:   case INTERPOLATION_L2_LOCAL:
2143:     SETERRQ1(PETSC_ERR_SUP, "Interpolation type %d not yet supported", grid->interpolationType);
2144:   default:
2145:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid interpolation type %d", grid->interpolationType);
2146:   }
2147:   /* Anything specific to the implementation */
2148:   if (grid->ops->gridreform != PETSC_NULL) {
2149:     (*grid->ops->gridreform)(grid, oldMesh, x, y);
2150:   }
2151:   PetscLogEventEnd(GRID_Reform, grid, 0, 0, 0);
2152:   return(0);
2153: }