#ifdef PETSC_RCS_HEADER static char vcid[] = "$Id: meshMovement.c,v 1.19 2000/10/17 13:48:57 knepley Exp $"; #endif #include "src/gvec/meshMovementImpl.h" /*I "meshMovement.h" I*/ #include "src/grid/gridimpl.h" /*I "grid.h" I*//*I "gvec.h" I*/ extern int MeshPreCopy_Private(Mesh); extern int MeshPostCopy_Private(Mesh, Mesh); /*-------------------------------------------------- Mesh Functions --------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "MeshSetMover" /*@ MeshSetMover - This function provides the MeshMover object for this mesh. Not collective Input Parameters: + mesh - The mesh - mover - The mover, or PETSC_NULL Level: intermediate .keywords: Mesh, movement, mover .seealso: MeshGetMover(), MeshSetMovement() @*/ int MeshSetMover(Mesh mesh, MeshMover mover) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE); ierr = PetscObjectCompose((PetscObject) mesh, "MeshMover", (PetscObject) mover); CHKERRQ(ierr); if (mover != PETSC_NULL) { PetscValidHeaderSpecific(mover, MESH_COOKIE); ierr = PetscObjectDereference((PetscObject) mover); CHKERRQ(ierr); ierr = MeshMoverSetMesh(mover, mesh); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetMover" /*@ MeshGetMover - This function returns the MeshMover object for this mesh. Not collective Input Parameter: . mesh - The mesh Output Parameters: . mover - The mover, or PETSC_NULL if it does not exist Level: intermediate .keywords: Mesh, movement, mover .seealso: MeshSetMover(), MeshGetMovement() @*/ int MeshGetMover(Mesh mesh, MeshMover *mover) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE); ierr = PetscObjectQuery((PetscObject) mesh, "MeshMover", (PetscObject *) mover); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoveMesh" /*@ MeshMoveMesh - Recalculates the node coordinates based on the current velocities and accelerations. Collective on Mesh Input Parameters: . mesh - The mesh . dt - The timestep Level: advanced .keywords: mesh, movement .seealso: MeshCalcNodeVelocitiess(), MeshCalcNodeAccelerations() @*/ int MeshMoveMesh(Mesh mesh, double dt) { MeshMover mover; PetscTruth isMoving; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE); ierr = MeshGetMovement(mesh, &isMoving); CHKERRQ(ierr); if (isMoving == PETSC_FALSE) PetscFunctionReturn(0); ierr = PetscLogEventBegin(MESH_MoveMesh, mesh, 0, 0, 0); CHKERRQ(ierr); ierr = MeshGetMover(mesh, &mover); CHKERRQ(ierr); ierr = (*mover->ops->movemesh)(mover, dt); CHKERRQ(ierr); ierr = PetscLogEventEnd(MESH_MoveMesh, mesh, 0, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshUpdateNodeValues" /*@ MeshUpdateNodeValues - Updates the mesh velocity, acceleration, or any other node-based quantity. Collective on Mesh Input Parameters: + mesh - The mesh - sol - The solution to an auxilliary problem Output Parameters: + vec - The vector which will hold he solution - ghostVec - The corresponding ghost vector for vec Level: developer .keywords: mesh, movement .seealso: MeshMoverSetNodeVelocities() @*/ int MeshUpdateNodeValues(Mesh mesh, GVec sol, Vec vec, Vec ghostVec) { MeshMover mover; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE); PetscValidHeaderSpecific(sol, VEC_COOKIE); PetscValidHeaderSpecific(vec, VEC_COOKIE); PetscValidHeaderSpecific(ghostVec, VEC_COOKIE); ierr = MeshGetMover(mesh, &mover); CHKERRQ(ierr); ierr = (*mover->ops->updatenodevalues)(mover, sol, vec, ghostVec); CHKERRQ(ierr); PetscFunctionReturn(0); } EXTERN_C_BEGIN #undef __FUNCT__ #define __FUNCT__ "MeshPostReform_Private" /* MeshPostReform_Private - Create a suitable MeshMover after reforming a Mesh. Collective on Mesh Input Parameters: + mesh - The mesh - newMesh - The reformed mesh Level: developer .keywords: mesh, refine, reform .seealso: MeshReform() */ int MeshPostReform_Private(Mesh mesh, Mesh newMesh) { MeshMover mover, newMover; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE); PetscValidHeaderSpecific(newMesh, MESH_COOKIE); ierr = MeshGetMover(mesh, &mover); CHKERRQ(ierr); if (mover != PETSC_NULL) { ierr = MeshMoverConvert(mover, newMesh, &newMover); CHKERRQ(ierr); } PetscFunctionReturn(0); } EXTERN_C_END #undef __FUNCT__ #define __FUNCT__ "MeshReform" /*@ MeshReform - Reform a mesh using the original boundary. Collective on Mesh Input Parameters: + mesh - The previous mesh . refine - Flag indicating whether to refine or recreate the mesh . area - A function which gives an area constraint when evaluated inside an element - newBd - Flag to determine whether the boundary should be generated (PETSC_TRUE) or taken from storage Output Parameter: . newMesh - The new mesh Level: beginner .keywords: mesh, refinement, reform .seealso: MeshReform(), MeshSetBoundary() @*/ int MeshReform(Mesh mesh, PetscTruth refine, PointFunction area, PetscTruth newBd, Mesh *newMesh) { PetscObject obj = (PetscObject) mesh; Mesh refinedMesh; MeshMover mover; PetscReal maxArea, startX, endX, startY, endY; void *tempCtx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE); PetscValidPointer(newMesh); ierr = PetscObjectComposeFunction(obj, "PostCopy", "MeshPostReform_Private", (void (*)(void)) MeshPostReform_Private);CHKERRQ(ierr); if (refine == PETSC_TRUE) { /* Refine grid with tolerance equal to the entire size This just "flips" the current mesh and does the necessary point insertion */ ierr = MeshGetBoundingBox(mesh, &startX, &startY, PETSC_NULL, &endX, &endY, PETSC_NULL); CHKERRQ(ierr); maxArea = (endX - startX)*(endY - startY); ierr = MeshGetUserContext(mesh, &tempCtx); CHKERRQ(ierr); ierr = MeshSetUserContext(mesh, (void *) &maxArea); CHKERRQ(ierr); ierr = MeshRefine(mesh, PointFunctionConstant, newMesh); CHKERRQ(ierr); ierr = MeshSetUserContext(mesh, tempCtx); CHKERRQ(ierr); ierr = MeshSetUserContext(*newMesh, tempCtx); CHKERRQ(ierr); } else { ierr = MeshGetMover(mesh, &mover); CHKERRQ(ierr); if (mover->ops->reform == PETSC_NULL) SETERRQ(PETSC_ERR_SUP, " "); ierr = PetscLogEventBegin(MESH_Reform, mesh, 0, 0, 0); CHKERRQ(ierr); ierr = MeshPreCopy_Private(mesh); CHKERRQ(ierr); ierr = (*mover->ops->reform)(mesh, refine, area, newBd, newMesh); CHKERRQ(ierr); ierr = MeshPostCopy_Private(mesh, *newMesh); CHKERRQ(ierr); ierr = PetscLogEventEnd(MESH_Reform, mesh, 0, 0, 0); CHKERRQ(ierr); } /* Refine new mesh if necessary */ if (area != PETSC_NULL) { PetscValidPointer(area); ierr = PetscObjectComposeFunction((PetscObject) *newMesh, "PostCopy", "MeshPostReform_Private", (void (*)(void)) MeshPostReform_Private);CHKERRQ(ierr); ierr = MeshRefine(*newMesh, area, &refinedMesh); CHKERRQ(ierr); ierr = MeshDestroy(*newMesh); CHKERRQ(ierr); *newMesh = refinedMesh; ierr = PetscObjectComposeFunction((PetscObject) *newMesh, "PostCopy", "", PETSC_NULL); CHKERRQ(ierr); } ierr = PetscObjectComposeFunction(obj, "PostCopy", "", PETSC_NULL); CHKERRQ(ierr); PetscFunctionReturn(0); } /*----------------------------------------------- MeshMover Functions ------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "MeshMoverCreate" /*@ MeshMoverCreate - This function creates an empty mesh mover. Collective on Mesh Input Parameter: . mesh - The Mesh to be moved Output Parameter: . mover - The mover Level: beginner .keywords: Mesh, create, mover .seealso: MeshSetUp(), MeshDestroy(), GridCreate() @*/ int MeshMoverCreate(Mesh mesh, MeshMover *mover) { MeshMover m; MPI_Comm comm; int ierr; PetscFunctionBegin; PetscValidPointer(mover); *mover = PETSC_NULL; ierr = PetscObjectGetComm((PetscObject) mesh, &comm); CHKERRQ(ierr); PetscHeaderCreate(m, _MeshMover, struct _MeshMoverOps, MESH_COOKIE, -1, "MeshMover", comm, MeshMoverDestroy, MeshMoverView); PetscLogObjectCreate(m); PetscLogObjectMemory(m, sizeof(struct _MeshMover)); ierr = PetscMemzero(m->ops, sizeof(struct _MeshMoverOps)); CHKERRQ(ierr); m->bops->publish = PETSC_NULL /* MeshMoverPublish_Petsc */; m->type_name = PETSC_NULL; m->serialize_name = PETSC_NULL; m->meshVelocityMethod = MESH_LAPLACIAN; m->meshVelocityGrid = PETSC_NULL; m->nodeVelocities = PETSC_NULL; m->nodeVelocitiesGhost = PETSC_NULL; m->meshVelocityMat = PETSC_NULL; m->meshVelocityRhs = PETSC_NULL; m->meshVelocitySol = PETSC_NULL; m->meshVelocitySles = PETSC_NULL; m->meshAccelerationMethod = MESH_LAPLACIAN; m->meshAccelerationGrid = PETSC_NULL; m->nodeAccelerations = PETSC_NULL; m->nodeAccelerationsGhost = PETSC_NULL; m->meshAccelerationMat = PETSC_NULL; m->meshAccelerationRhs = PETSC_NULL; m->meshAccelerationSol = PETSC_NULL; m->meshAccelerationSles = PETSC_NULL; m->meshVelocityScatter = PETSC_NULL; m->movingMeshViewer = PETSC_NULL; m->movingMeshViewerCaption = "Moving Mesh"; m->movingMeshCtx = PETSC_NULL; ierr = MeshSetMover(mesh, m); CHKERRQ(ierr); *mover = m; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetup" /*@ MeshMoverSetup - Creates all the necessary structures for mesh movement. Collective on MeshMover Input Parameter: . mover - The mover Note: This function should be called directly before beginning the calculation, and after all options have been set. It is normally called automatically by the solver. Level: advanced .keywords: mesh, movement .seealso: MeshMoverGetMovement() @*/ int MeshMoverSetup(MeshMover mover) { Mesh mesh; Partition part; IS globalIS, localIS; int dim; int *indices; int numNodes, numOverlapNodes; int node, gNode, j; PetscTruth reduceVel, reduceAcc; PetscTruth opt; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); /* Destroy old structures */ if (mover->nodeVelocities != PETSC_NULL) { ierr = VecDestroy(mover->nodeVelocities); CHKERRQ(ierr); ierr = VecDestroy(mover->nodeVelocitiesGhost); CHKERRQ(ierr); } if (mover->meshVelocityGrid != PETSC_NULL) { ierr = GridDestroy(mover->meshVelocityGrid); CHKERRQ(ierr); mover->meshVelocityGrid = PETSC_NULL; ierr = GMatDestroy(mover->meshVelocityMat); CHKERRQ(ierr); mover->meshVelocityMat = PETSC_NULL; ierr = GVecDestroy(mover->meshVelocityRhs); CHKERRQ(ierr); mover->meshVelocityRhs = PETSC_NULL; ierr = GVecDestroy(mover->meshVelocitySol); CHKERRQ(ierr); mover->meshVelocitySol = PETSC_NULL; ierr = SLESDestroy(mover->meshVelocitySles); CHKERRQ(ierr); mover->meshVelocitySles = PETSC_NULL; } if (mover->nodeAccelerations != PETSC_NULL) { ierr = VecDestroy(mover->nodeAccelerations); CHKERRQ(ierr); ierr = VecDestroy(mover->nodeAccelerationsGhost); CHKERRQ(ierr); } if (mover->meshAccelerationGrid != PETSC_NULL) { ierr = GridDestroy(mover->meshAccelerationGrid); CHKERRQ(ierr); mover->meshAccelerationGrid = PETSC_NULL; ierr = GMatDestroy(mover->meshAccelerationMat); CHKERRQ(ierr); mover->meshAccelerationMat = PETSC_NULL; ierr = GVecDestroy(mover->meshAccelerationRhs); CHKERRQ(ierr); mover->meshAccelerationRhs = PETSC_NULL; ierr = GVecDestroy(mover->meshAccelerationSol); CHKERRQ(ierr); mover->meshAccelerationSol = PETSC_NULL; ierr = SLESDestroy(mover->meshAccelerationSles); CHKERRQ(ierr); mover->meshAccelerationSles = PETSC_NULL; } if (mover->meshVelocityScatter != PETSC_NULL) { ierr = VecScatterDestroy(mover->meshVelocityScatter); CHKERRQ(ierr); } /* Determine solution method */ ierr = PetscOptionsHasName(mover->prefix, "-mesh_move_lap", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { mover->meshVelocityMethod = MESH_LAPLACIAN; mover->meshAccelerationMethod = MESH_LAPLACIAN; } ierr = PetscOptionsHasName(mover->prefix, "-mesh_move_weighted_lap", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { mover->meshVelocityMethod = MESH_WEIGHTED_LAPLACIAN; mover->meshAccelerationMethod = MESH_WEIGHTED_LAPLACIAN; } /* Setup mesh node velocities vector */ ierr = MeshMoverGetMesh(mover, &mesh); CHKERRQ(ierr); ierr = MeshGetPartition(mesh, &part); CHKERRQ(ierr); ierr = MeshGetDimension(mesh, &dim); CHKERRQ(ierr); ierr = MeshGetInfo(mesh, PETSC_NULL, &numNodes, PETSC_NULL, PETSC_NULL); CHKERRQ(ierr); ierr = VecCreateMPI(mover->comm, numNodes*dim, PETSC_DECIDE, &mover->nodeVelocities); CHKERRQ(ierr); ierr = PartitionGetNumOverlapNodes(part, &numOverlapNodes); CHKERRQ(ierr); ierr = VecCreateSeq(PETSC_COMM_SELF, numOverlapNodes*dim, &mover->nodeVelocitiesGhost); CHKERRQ(ierr); /* Setup mesh node accelerations vector */ ierr = VecCreateMPI(mover->comm, numNodes*dim, PETSC_DECIDE, &mover->nodeAccelerations); CHKERRQ(ierr); ierr = VecCreateSeq(PETSC_COMM_SELF, numOverlapNodes*dim, &mover->nodeAccelerationsGhost); CHKERRQ(ierr); /* Setup grids and other implementation specific stuff */ ierr = (*mover->ops->setup)(mover); CHKERRQ(ierr); /* Setup scatter from node values vector to ghost vector */ ierr = ISCreateStride(mover->comm, numOverlapNodes*dim, 0, 1, &localIS); CHKERRQ(ierr); ierr = PetscMalloc(numOverlapNodes*dim * sizeof(int), &indices); CHKERRQ(ierr); for(node = 0; node < numOverlapNodes; node++) { ierr = PartitionLocalToGlobalNodeIndex(part, node, &gNode); CHKERRQ(ierr); for(j = 0; j < dim; j++) indices[node*dim+j] = gNode*dim+j; } ierr = ISCreateGeneral(mover->comm, numOverlapNodes*dim, indices, &globalIS); CHKERRQ(ierr); ierr = VecScatterCreate(mover->nodeVelocities, globalIS, mover->nodeVelocitiesGhost, localIS, &mover->meshVelocityScatter); CHKERRQ(ierr); ierr = ISDestroy(localIS); CHKERRQ(ierr); ierr = ISDestroy(globalIS); CHKERRQ(ierr); ierr = PetscFree(indices); CHKERRQ(ierr); /* Setup viewer */ ierr = PetscOptionsHasName(mover->prefix, "-mesh_moving_view", &opt); CHKERRQ(ierr); if ((opt == PETSC_TRUE) && (mover->movingMeshViewer == PETSC_NULL)) { mover->movingMeshViewer = PETSC_VIEWER_STDOUT_WORLD; } ierr = PetscOptionsHasName(mover->prefix, "-mesh_moving_view_draw", &opt); CHKERRQ(ierr); if ((opt == PETSC_TRUE) && (mover->movingMeshViewer == PETSC_NULL)) { PetscReal startX, startY, endX, endY, sizeX, sizeY; int xsize, ysize; ierr = MeshGetBoundingBox(mesh, &startX, &startY, 0, &endX, &endY, 0); CHKERRQ(ierr); sizeX = endX - startX; sizeY = endY - startY; if (sizeX > sizeY) { ysize = 300; xsize = ysize * (int) (sizeX/sizeY); } else { xsize = 300; ysize = xsize * (int) (sizeY/sizeX); } ierr = PetscViewerDrawOpen(mover->comm, 0, 0, 0, 315, xsize, ysize, &mover->movingMeshViewer); CHKERRQ(ierr); } /* Check grids */ ierr = GridGetReduceSystem(mover->meshVelocityGrid, &reduceVel); CHKERRQ(ierr); ierr = GridGetReduceSystem(mover->meshAccelerationGrid, &reduceAcc); CHKERRQ(ierr); if (reduceVel != reduceAcc) { SETERRQ(PETSC_ERR_ARG_WRONG, "Velocity and acceleration grids must have same BC reduction"); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetTypeFromOptions" /* MeshMoverSetTypeFromOptions - Sets the type of mesh movement from user options. Collective on MeshMover Input Parameter: . mover - The mover Level: intermediate .keywords: MeshMover, set, options, database, type .seealso: MeshMoverSetFromOptions(), MeshMoverSetType() */ static int MeshMoverSetTypeFromOptions(MeshMover mover) { Mesh mesh; PetscTruth opt; char *defaultType; char typeName[256]; int dim; int ierr; PetscFunctionBegin; ierr = MeshMoverGetMesh(mover, &mesh); CHKERRQ(ierr); ierr = MeshGetDimension(mesh, &dim); CHKERRQ(ierr); if (mover->type_name != PETSC_NULL) { defaultType = mover->type_name; } else { switch(dim) { case 2: defaultType = MESH_TRIANGULAR_2D; break; default: SETERRQ1(PETSC_ERR_SUP, "MeshMover dimension %d is not supported", dim); } } if (!MeshMoverRegisterAllCalled) { ierr = MeshMoverRegisterAll(PETSC_NULL); CHKERRQ(ierr); } ierr = PetscOptionsList("-mesh_mover_type", "Mesh movement method", "MeshMoverSetType", MeshMoverList, defaultType, typeName, 256, &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = MeshMoverSetType(mover, typeName); CHKERRQ(ierr); } else { ierr = MeshMoverSetType(mover, defaultType); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetFromOptions" /*@ MeshMoverSetFromOptions - Sets various MeshMover parameters from user options. Collective on MeshMover Input Parameter: . mover - The mover Notes: To see all options, run your program with the -help option, or consult the users manual. Must be called after MeshMoverCreate() but before the vector is used. Level: intermediate .keywords: MeshMover, set, options, database .seealso: MeshMoverCreate(), MeshMoverPrintHelp(), MeshMoverSetOptionsPrefix() @*/ int MeshMoverSetFromOptions(MeshMover mover) { PetscTruth opt; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); ierr = PetscOptionsBegin(mover->comm, mover->prefix, "MeshMover options", "MeshMover"); CHKERRQ(ierr); /* Handle generic mesh options */ ierr = PetscOptionsHasName(PETSC_NULL, "-help", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = MeshMoverPrintHelp(mover); CHKERRQ(ierr); } /* Handle mesh type options */ ierr = MeshMoverSetTypeFromOptions(mover); CHKERRQ(ierr); /* Handle specific mesh options */ if (mover->ops->setfromoptions != PETSC_NULL) { ierr = (*mover->ops->setfromoptions)(mover); CHKERRQ(ierr); } ierr = PetscOptionsEnd(); CHKERRQ(ierr); /* Handle subobject options */ PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverPrintHelp" /*@ MeshMoverPrintHelp - Prints all options for the MeshMover. Input Parameter: . mover - The mover Options Database Keys: $ -help, -h Level: intermediate .keywords: MeshMover, help .seealso: MeshMoverSetFromOptions() @*/ int MeshMoverPrintHelp(MeshMover mover) { char p[64]; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); ierr = PetscStrcpy(p, "-"); CHKERRQ(ierr); if (mover->prefix != PETSC_NULL) { ierr = PetscStrcat(p, mover->prefix); CHKERRQ(ierr); } (*PetscHelpPrintf)(mover->comm, "MeshMover options ------------------------------------------------\n"); (*PetscHelpPrintf)(mover->comm," %smesh_mover_type : Sets the mesh movement type\n", p); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverDuplicate" /*@ MeshMoverDuplicate - This function returns an identical MeshMover. Not collective Input Parameter: . mover - The mover Output Parameter: . newMover - The new mover Level: intermediate .keywords: Mesh, movement, mover .seealso: MeshMoverCopy() @*/ int MeshMoverDuplicate(MeshMover mover, MeshMover *newMover) { Mesh mesh; int ierr; PetscFunctionBegin; ierr = MeshMoverGetMesh(mover, &mesh); CHKERRQ(ierr); ierr = MeshMoverConvert(mover, mesh, newMover); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverConvert" /*@ MeshMoverConvert - This function returns an identical MeshMover on the new Mesh. Not collective Input Parameters: + mover - The mover - newMesh - The new Mesh Output Parameter: . newMover - The new mover Level: intermediate .keywords: Mesh, movement, mover .seealso: MeshMoverCopy() @*/ int MeshMoverConvert(MeshMover mover, Mesh newMesh, MeshMover *newMover) { MeshMover m; Mesh mesh; int ierr; PetscFunctionBegin; ierr = MeshMoverGetMesh(mover, &mesh); CHKERRQ(ierr); if (mesh == newMesh) { /* Need to play reference games since MeshMoverCreate automatically sets it in the Mesh */ ierr = PetscObjectReference((PetscObject) mover); CHKERRQ(ierr); ierr = MeshMoverCreate(newMesh, &m); CHKERRQ(ierr); ierr = PetscObjectReference((PetscObject) m); CHKERRQ(ierr); ierr = MeshSetMover(mesh, mover); CHKERRQ(ierr); } else { ierr = MeshMoverCreate(newMesh, &m); CHKERRQ(ierr); } ierr = MeshMoverSetFromOptions(m); CHKERRQ(ierr); m->meshVelocityMethod = mover->meshVelocityMethod; m->meshAccelerationMethod = mover->meshAccelerationMethod; m->movingMeshViewer = mover->movingMeshViewer; if (mover->movingMeshViewer != PETSC_NULL) PetscObjectReference((PetscObject) mover->movingMeshViewer); m->movingMeshViewerCaption = mover->movingMeshViewerCaption; m->movingMeshCtx = mover->movingMeshCtx; ierr = MeshMoverSetup(m); CHKERRQ(ierr); ierr = GridDuplicateBC(mover->meshVelocityGrid, m->meshVelocityGrid); CHKERRQ(ierr); ierr = GridDuplicateBC(mover->meshAccelerationGrid, m->meshAccelerationGrid); CHKERRQ(ierr); *newMover = m; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverView" /*@ MeshMoverView - Views a mesh mover object. Collective on MeshMover Input Parameters: + mover - The mover - viewer - The viewer with which to view the mesh mover Level: beginner .keywords: mesh, view, mover .seealso: MeshMoverDestroy(), PetscViewerDrawOpen() @*/ int MeshMoverView(MeshMover mover, PetscViewer viewer) { PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); if (viewer == PETSC_NULL) { viewer = PETSC_VIEWER_STDOUT_SELF; } else { PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverDestroy" /*@ MeshMoverDestroy - Destroys a mesh mover object. Collective on MeshMover Input Parameter: . mover - The mover Level: beginner .keywords: mesh, destroy, mover .seealso MeshMoverCreate() @*/ int MeshMoverDestroy(MeshMover mover) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); if (--mover->refct > 0) PetscFunctionReturn(0); if (mover->nodeVelocities != PETSC_NULL) { ierr = GVecDestroy(mover->nodeVelocities); CHKERRQ(ierr); ierr = GVecDestroy(mover->nodeVelocitiesGhost); CHKERRQ(ierr); } if (mover->meshVelocitySles != PETSC_NULL) { ierr = GMatDestroy(mover->meshVelocityMat); CHKERRQ(ierr); ierr = GVecDestroy(mover->meshVelocityRhs); CHKERRQ(ierr); ierr = GVecDestroy(mover->meshVelocitySol); CHKERRQ(ierr); ierr = SLESDestroy(mover->meshVelocitySles); CHKERRQ(ierr); } if (mover->meshVelocityGrid != PETSC_NULL) { ierr = GridFinalizeBC(mover->meshVelocityGrid); CHKERRQ(ierr); ierr = GridDestroy(mover->meshVelocityGrid); CHKERRQ(ierr); } if (mover->nodeAccelerations != PETSC_NULL) { ierr = GVecDestroy(mover->nodeAccelerations); CHKERRQ(ierr); ierr = GVecDestroy(mover->nodeAccelerationsGhost); CHKERRQ(ierr); } if (mover->meshAccelerationSles != PETSC_NULL) { ierr = GMatDestroy(mover->meshAccelerationMat); CHKERRQ(ierr); ierr = GVecDestroy(mover->meshAccelerationRhs); CHKERRQ(ierr); ierr = GVecDestroy(mover->meshAccelerationSol); CHKERRQ(ierr); ierr = SLESDestroy(mover->meshAccelerationSles); CHKERRQ(ierr); } if (mover->meshAccelerationGrid != PETSC_NULL) { ierr = GridFinalizeBC(mover->meshAccelerationGrid); CHKERRQ(ierr); ierr = GridDestroy(mover->meshAccelerationGrid); CHKERRQ(ierr); } if (mover->meshVelocityScatter != PETSC_NULL) { ierr = VecScatterDestroy(mover->meshVelocityScatter); CHKERRQ(ierr); } if (mover->movingMeshViewer != PETSC_NULL) { ierr = PetscViewerDestroy(mover->movingMeshViewer); CHKERRQ(ierr); } ierr = (*mover->ops->destroy)(mover); CHKERRQ(ierr); PetscLogObjectDestroy(mover); PetscHeaderDestroy(mover); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetMesh" /*@ MeshMoverSetMesh - This function returns the associated Mesh. Not collective Input Parameter: . mover - The mover Output Parameter: . mesh - The mesh, or PETSC_NULL if it does not exist Level: intermediate .keywords: Mesh, movement, mover .seealso: MeshSetMover() @*/ int MeshMoverSetMesh(MeshMover mover, Mesh mesh) { PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); PetscValidHeaderSpecific(mesh, MESH_COOKIE); mover->mesh = mesh; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverGetMesh" /*@ MeshMoverGetMesh - This function returns the associated Mesh. Not collective Input Parameter: . mover - The mover Output Parameter: . mesh - The mesh, or PETSC_NULL if it does not exist Level: intermediate .keywords: Mesh, movement, mover .seealso: MeshGetMover() @*/ int MeshMoverGetMesh(MeshMover mover, Mesh *mesh) { PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); PetscValidPointer(mesh); *mesh = mover->mesh; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverGetMovement" /*@ MeshMoverGetMovement - Returns the type of mesh movement calculation. Not collective Input Parameter: . mover - The mover Output Parameters: + vtype - The type of mesh velocity calculation, like MESH_LAPLACIAN . atype - The type of mesh acceleration calculation, like MESH_LAPLACIAN - ctx - A user context Level: intermediate .keywords: mesh, movement .seealso: MeshMoverSetMovement() @*/ int MeshMoverGetMovement(MeshMover mover, MeshSolveMethod *vtype, MeshSolveMethod *atype, PetscObject *ctx) { PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); if (vtype != PETSC_NULL) { PetscValidPointer(vtype); *vtype = mover->meshVelocityMethod; } if (atype != PETSC_NULL) { PetscValidPointer(atype); *atype = mover->meshAccelerationMethod; } if (ctx != PETSC_NULL) { PetscValidPointer(ctx); *ctx = mover->movingMeshCtx; } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetMovement" /*@ MeshMoverSetMovement - Determines the type of mesh mvoement calculation Collective on MeshMover Input Parameters: + mover - The mover . vtype - The type of mesh velocity calculation, like MESH_LAPLACIAN . atype - The type of mesh acceleration calculation, like MESH_LAPLACIAN - ctx - A user context Level: intermediate .keywords: mesh, movement .seealso: MeshMoverGetMovement() @*/ int MeshMoverSetMovement(MeshMover mover, MeshSolveMethod vtype, MeshSolveMethod atype, PetscObject ctx) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); if (vtype != PETSC_DECIDE) mover->meshVelocityMethod = vtype; if (atype != PETSC_DECIDE) mover->meshAccelerationMethod = atype; if (ctx != PETSC_NULL) mover->movingMeshCtx = ctx; if (mover->meshVelocityGrid == PETSC_NULL) { ierr = MeshMoverSetup(mover); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverGetVelocityGrid" /*@ MeshMoverGetVelocityGrid - Returns the grid which defines the mesh velocity. Not collective Input Parameter: . mover - The mover Output Parameter: . velGrid - The mesh velocity grid Level: intermediate .keywords: grid, mesh, mesh velocity, movement .seealso: MeshMoverGetAccelerationGrid(), GridGetMeshMovement() @*/ int MeshMoverGetVelocityGrid(MeshMover mover, Grid *velGrid) { PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); *velGrid = mover->meshVelocityGrid; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverGetAccelerationGrid" /*@ MeshMoverGetAccelerationGrid - Returns the grid which defines the mesh acceleration. Not collective Input Parameter: . mover - The mover Output Parameter: . accGrid - The mesh acceleration grid Level: intermediate .keywords: grid, mesh, mesh acceleration, movement .seealso: MeshMoverGetVelocityGrid(), GridGetMeshMovement() @*/ int MeshMoverGetAccelerationGrid(MeshMover mover, Grid *accGrid) { PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); *accGrid = mover->meshAccelerationGrid; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverCalcVelocityBCValues" /*@ MeshMoverCalcVelocityBCValues - This function calculates the boundary values for the mesh velocity. It is normally called once a timestep when using time dependent boundary conditions. Collective on MeshMover Input Parameters: + mover - The MeshMover . save - A flag used to store old values, usually for timestepping - ctx - The context Level: advanced .keywords: MeshMover, reduction, boundary conditions .seealso: MeshMoverSetVelocityBCContext(), MeshMoverSetVelocityBC() @*/ int MeshMoverCalcVelocityBCValues(MeshMover mover) { Grid g; PetscTruth reduce; void *reduceCtx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); ierr = MeshMoverGetMovement(mover, PETSC_NULL, PETSC_NULL, PETSC_NULL); CHKERRQ(ierr); ierr = MeshMoverGetVelocityGrid(mover, &g); CHKERRQ(ierr); PetscValidHeaderSpecific(g, GRID_COOKIE); ierr = GridGetReduceSystem(g, &reduce); CHKERRQ(ierr); if (reduce == PETSC_FALSE) PetscFunctionReturn(0); ierr = GridSetUp(g); CHKERRQ(ierr); ierr = GridGetBCContext(g, &reduceCtx); CHKERRQ(ierr); ierr = GridCalcBCValues(g, PETSC_FALSE, reduceCtx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverCalcAccelerationBCValues" /*@ MeshMoverCalcAccelerationBCValues - This function calculates the boundary values for the mesh acceleration. It is normally called once a timestep when using time dependent boundary conditions. Collective on MeshMover Input Parameters: + mover - The MeshMover . save - A flag used to store old values, usually for timestepping - ctx - The context Level: advanced .keywords: MeshMover, reduction, boundary conditions .seealso: MeshMoverSetAccelerationBCContext(), MeshMoverSetAccelerationBC() @*/ int MeshMoverCalcAccelerationBCValues(MeshMover mover) { Grid g; PetscTruth reduce; void *reduceCtx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); ierr = MeshMoverGetMovement(mover, PETSC_NULL, PETSC_NULL, PETSC_NULL); CHKERRQ(ierr); ierr = MeshMoverGetAccelerationGrid(mover, &g); CHKERRQ(ierr); PetscValidHeaderSpecific(g, GRID_COOKIE); ierr = GridGetReduceSystem(g, &reduce); CHKERRQ(ierr); if (reduce == PETSC_FALSE) PetscFunctionReturn(0); ierr = GridSetUp(g); CHKERRQ(ierr); ierr = GridGetBCContext(g, &reduceCtx); CHKERRQ(ierr); ierr = GridCalcBCValues(g, PETSC_FALSE, reduceCtx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverCalcNodeVelocities" /*@ MeshMoverCalcNodeVelocities - Recalculates the node velocities. Collective on MeshMover Input Parameters: + mover - The mover - flag - Flag for reuse of old operator in calculation Level: advanced .keywords: mesh, movement .seealso: MeshMoverCalcNodeAccelerations() @*/ int MeshMoverCalcNodeVelocities(MeshMover mover, PetscTruth flag) { Grid grid = mover->meshVelocityGrid; MatStructure matStructFlag = DIFFERENT_NONZERO_PATTERN; Mesh mesh; SLES sles; GMat mat; GVec rhs, sol; PetscScalar zero = 0.0; int its; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); ierr = PetscLogEventBegin(MESH_CalcNodeVelocities, mover, 0, 0, 0); CHKERRQ(ierr); /* Allocation */ if (mover->meshVelocitySles == PETSC_NULL) { ierr = GVecCreate(grid, &mover->meshVelocityRhs); CHKERRQ(ierr); PetscLogObjectParent(mover, mover->meshVelocityRhs); ierr = GVecCreate(grid, &mover->meshVelocitySol); CHKERRQ(ierr); PetscLogObjectParent(mover, mover->meshVelocitySol); ierr = GMatCreate(grid, &mover->meshVelocityMat); CHKERRQ(ierr); PetscLogObjectParent(mover, mover->meshVelocityMat); ierr = SLESCreate(mover->comm, &mover->meshVelocitySles); CHKERRQ(ierr); ierr = SLESAppendOptionsPrefix(mover->meshVelocitySles, "mesh_vel_"); CHKERRQ(ierr); ierr = SLESSetFromOptions(mover->meshVelocitySles); CHKERRQ(ierr); PetscLogObjectParent(mover, mover->meshVelocitySles); } sles = mover->meshVelocitySles; mat = mover->meshVelocityMat; rhs = mover->meshVelocityRhs; sol = mover->meshVelocitySol; /* Form System Matrix */ if (flag == PETSC_TRUE) { ierr = MatZeroEntries(mat); CHKERRQ(ierr); ierr = GridEvaluateSystemMatrix(grid, PETSC_NULL, &mat, &mat, &matStructFlag, mover->movingMeshCtx); CHKERRQ(ierr); ierr = GMatSetBoundary(mat, 1.0, mover->movingMeshCtx); CHKERRQ(ierr); ierr = MatSetOption(mat, MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRQ(ierr); } /* Form Rhs */ ierr = VecSet(&zero, rhs); CHKERRQ(ierr); ierr = GridEvaluateRhs(grid, PETSC_NULL, rhs, mover->movingMeshCtx); CHKERRQ(ierr); ierr = GVecSetBoundary(rhs, mover->movingMeshCtx); CHKERRQ(ierr); /* Solve system */ ierr = SLESSetOperators(sles, mat, mat, SAME_NONZERO_PATTERN); CHKERRQ(ierr); ierr = SLESSolve(sles, rhs, sol, &its); CHKERRQ(ierr); /* Copy values into the mesh vector */ ierr = MeshMoverGetMesh(mover, &mesh); CHKERRQ(ierr); ierr = MeshUpdateNodeValues(mesh, sol, mover->nodeVelocities, mover->nodeVelocitiesGhost); CHKERRQ(ierr); ierr = PetscLogEventEnd(MESH_CalcNodeVelocities, mover, 0, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverCalcNodeAccelerations" /*@ MeshMoverCalcNodeAccelerations - Recalculates the node accelerations. Collective on MeshMover Input Parameters: + mover - The mover - flag - Flag for reuse of old operator in calculation Level: advanced .keywords: mesh, movement .seealso: MeshMoverCalcNodeVelocities() @*/ int MeshMoverCalcNodeAccelerations(MeshMover mover, PetscTruth flag) { Grid grid = mover->meshAccelerationGrid; MatStructure matStructFlag = DIFFERENT_NONZERO_PATTERN; Mesh mesh; SLES sles; GMat mat; GVec rhs, sol; PetscScalar zero = 0.0; int its; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); ierr = PetscLogEventBegin(MESH_CalcNodeAccelerations, mover, 0, 0, 0); CHKERRQ(ierr); /* Allocation */ if (mover->meshAccelerationSles == PETSC_NULL) { ierr = GVecCreate(grid, &mover->meshAccelerationRhs); CHKERRQ(ierr); PetscLogObjectParent(mover, mover->meshAccelerationRhs); ierr = GVecCreate(grid, &mover->meshAccelerationSol); CHKERRQ(ierr); PetscLogObjectParent(mover, mover->meshAccelerationSol); ierr = GMatCreate(grid, &mover->meshAccelerationMat); CHKERRQ(ierr); PetscLogObjectParent(mover, mover->meshAccelerationMat); ierr = SLESCreate(mover->comm, &mover->meshAccelerationSles); CHKERRQ(ierr); ierr = SLESAppendOptionsPrefix(mover->meshAccelerationSles, "mesh_accel_"); CHKERRQ(ierr); ierr = SLESSetFromOptions(mover->meshAccelerationSles); CHKERRQ(ierr); PetscLogObjectParent(mover, mover->meshAccelerationSles); } sles = mover->meshAccelerationSles; mat = mover->meshAccelerationMat; rhs = mover->meshAccelerationRhs; sol = mover->meshAccelerationSol; /* Form System Matrix */ if (flag == PETSC_TRUE) { ierr = MatZeroEntries(mat); CHKERRQ(ierr); ierr = GridEvaluateSystemMatrix(grid, PETSC_NULL, &mat, &mat, &matStructFlag, mover->movingMeshCtx); CHKERRQ(ierr); ierr = GMatSetBoundary(mat, 1.0, mover->movingMeshCtx); CHKERRQ(ierr); ierr = MatSetOption(mat, MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRQ(ierr); } /* Form Rhs */ ierr = VecSet(&zero, rhs); CHKERRQ(ierr); ierr = GridEvaluateRhs(grid, PETSC_NULL, rhs, mover->movingMeshCtx); CHKERRQ(ierr); ierr = GVecSetBoundary(rhs, mover->movingMeshCtx); CHKERRQ(ierr); /* Solve system */ ierr = SLESSetOperators(sles, mat, mat, SAME_NONZERO_PATTERN); CHKERRQ(ierr); ierr = SLESSolve(sles, rhs, sol, &its); CHKERRQ(ierr); /* Copy values into the mesh vector */ ierr = MeshMoverGetMesh(mover, &mesh); CHKERRQ(ierr); ierr = MeshUpdateNodeValues(mesh, sol, mover->nodeAccelerations, mover->nodeAccelerationsGhost); CHKERRQ(ierr); ierr = PetscLogEventEnd(MESH_CalcNodeAccelerations, mover, 0, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetNodeVelocities" /*@ MeshMoverSetNodeVelocities - Allows the user to specify node velocities. Collective on MeshMover Input Parameters: + mover - The mover . func - A PointFunction defining the node velocities . alpha - A scalar multiplier - ctx - A user context Level: advanced .keywords: mesh, movement .seealso: MeshMoverSetNodeAccelerations() @*/ int MeshMoverSetNodeVelocities(MeshMover mover, PointFunction func, PetscScalar alpha, void *ctx) { #if 0 int field = 0; int ierr; PetscFunctionBegin; ierr = GVecEvaluateFunction(mover->nodeVelocities, 1, &field, func, alpha, ctx); CHKERRQ(ierr); /* Now need a way to get mesh velocities back to the system */ PetscFunctionReturn(0); #else SETERRQ(PETSC_ERR_SUP, "Not yet a grid vector"); #endif } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetNodeAccelerations" /*@ MeshMoverSetNodeAccelerations - Allows the user to specify node accelerations. Collective on MeshMover Input Parameters: + mover - The mover . func - A function defining the node accelerations . alpha - A scalar multiplier - ctx - A user context Level: intermediate .keywords: mesh, movement .seealso: MeshMoverSetNodeVelocities() @*/ int MeshMoverSetNodeAccelerations(MeshMover mover, PointFunction func, PetscScalar alpha, void *ctx) { #if 0 int field = 0; int ierr; PetscFunctionBegin; ierr = GVecEvaluateFunction(mover->nodeAccelerations, 1, &field, func, alpha, ctx); CHKERRQ(ierr); /* Now need a way to get mesh accelerations back to the system */ PetscFunctionReturn(0); #else SETERRQ(PETSC_ERR_SUP, "Not yet a grid vector"); #endif } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetVelocityBC" /*@C MeshMoverSetVelocityBC - This function sets the boundary conditions to use for the mesh velocity calculation. Collective on MeshMover Input Parameters: + mover - The mover . bd - The marker for the boundary along which conditions are applied . f - The function which defines the boundary condition - reduce - The flag for explicit reduction of the system Note: This function clears any previous conditions, whereas MeshMoverAddVelocityBC() appends the condition. Level: intermediate .keywords mesh, velocity, boundary condition .seealso MeshMoverAddVelocityBC(), MeshMoverSetAccelerationBC(), GridSetBC(), GridAddBC() @*/ int MeshMoverSetVelocityBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); if (mover->meshVelocityGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Velocity grid not created"); ierr = GridSetBC(mover->meshVelocityGrid, bd, 0, f, reduce); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverAddVelocityBC" /*@C MeshMoverAddVelocityBC - This function adds a boundary condition to use for the mesh velocity calculation. Collective on MeshMover Input Parameters: + mover - The mover . bd - The marker for the boundary along which conditions are applied . f - The function which defines the boundary condition - reduce - The flag for explicit reduction of the system Note: This function appends the condition, whereas MeshMoverSetVelocityBC() clears any previous conditions. Level: intermediate .keywords mesh, velocity, boundary condition .seealso MeshMoverSetVelocityBC(), MeshMoverAddAccelerationBC(), GridSetBC(), GridAddBC() @*/ int MeshMoverAddVelocityBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); if (mover->meshVelocityGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Velocity grid not created"); ierr = GridAddBC(mover->meshVelocityGrid, bd, 0, f, reduce); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetVelocityBCContext" /*@C MeshMoverSetVelocityBCContext - This function sets the boundary condition context to use for the mesh velocity calculation. Collective on MeshMover Input Parameters: + mover - The mover - ctx - The user context Level: intermediate .keywords mesh, velocity, boundary condition .seealso MeshMoverAddVelocityBC(), MeshMoverSetVelocityBC() @*/ int MeshMoverSetVelocityBCContext(MeshMover mover, void *ctx) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); if (mover->meshVelocityGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Velocity grid not created"); ierr = GridSetBCContext(mover->meshVelocityGrid, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetAccelerationBC" /*@C MeshMoverSetAccelerationBC - This function sets the boundary condition to use for the mesh acceleration calculation. Collective on MeshMover Input Parameters: + mover - The mover . bd - The marker for the boundary along which conditions are applied . f - The function which defines the boundary condition - reduce - The flag for explicit reduction of the system Note: This function clears any previous conditions, whereas MeshMoverAddAccelerationBC() appends the condition. Level: intermediate .keywords mesh, velocity, boundary condition .seealso MeshMoverAddAccelerationBC(), MeshMoverSetVelocityBC(), GridSetBC(), GridAddBC() @*/ int MeshMoverSetAccelerationBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); if (mover->meshAccelerationGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Acceleration grid not created"); ierr = GridSetBC(mover->meshAccelerationGrid, bd, 0, f, reduce); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverAddAccelerationBC" /*@C MeshMoverAddAccelerationBC - This function adds a boundary condition to use for the mesh acceleration calculation. Collective on MeshMover Input Parameters: + mover - The mover . bd - The marker for the boundary along which conditions are applied . f - The function which defines the boundary condition - reduce - The flag for explicit reduction of the system Note: This function appends the condition, whereas MeshMoverSetVelocityBC() clears any previous conditions. Level: intermediate .keywords mesh, velocity, boundary condition .seealso MeshMoverSetAccelerationBC(), MeshMoverAddVelocityBC(), GridSetBC(), GridAddBC() @*/ int MeshMoverAddAccelerationBC(MeshMover mover, int bd, PointFunction f, PetscTruth reduce) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); if (mover->meshAccelerationGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Acceleration grid not created"); ierr = GridAddBC(mover->meshAccelerationGrid, bd, 0, f, reduce); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetAccelerationBCContext" /*@C MeshMoverSetAccelerationBCContext - This function sets the boundary condition context to use for the mesh acceleration calculation. Collective on MeshMover Input Parameter: + mover - The mover - ctx - The user context Level: intermediate .keywords mesh, velocity, boundary condition .seealso MeshMoverAddAccelerationBC(), MeshMoverSetAccelerationBC() @*/ int MeshMoverSetAccelerationBCContext(MeshMover mover, void *ctx) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); if (mover->meshAccelerationGrid == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Acceleration grid not created"); ierr = GridSetBCContext(mover->meshAccelerationGrid, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetMovementCaption" /*@ MeshMoverSetMovementCaption - Sets the caption of the moving mesh viewer. Can be used to provide auxilliary information, e.g. time. Collective on MeshMover Input Parameters: + mover - The mover - caption - The caption Level: intermediate .keywords: mesh, movement .seealso: MeshMoverMoveMesh() @*/ int MeshMoverSetMovementCaption(MeshMover mover, char *caption) { PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); PetscValidPointer(caption); mover->movingMeshViewerCaption = caption; PetscFunctionReturn(0); } /*----------------------------------------- MeshMover Registration Functions -----------------------------------------*/ PetscFList MeshMoverList = 0; int MeshMoverRegisterAllCalled = 0; #undef __FUNCT__ #define __FUNCT__ "MeshMoverSetType" /*@C MeshMoverSetType - Builds a MeshMover, for a particular Mesh implementation. Collective on MeshMover Input Parameters: + mover - The MeshMover - type - The type name Options Database Key: . -mesh_mover_type - Sets the mesh mvoer type Level: intermediate .keywords: vector, set, type .seealso: MeshMoverCreate() @*/ int MeshMoverSetType(MeshMover mover, MeshMoverType method) { int (*r)(MeshMover); PetscTruth match; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); ierr = PetscTypeCompare((PetscObject) mover, method, &match); CHKERRQ(ierr); if (match == PETSC_TRUE) PetscFunctionReturn(0); /* Get the function pointers for the vector requested */ if (!MeshMoverRegisterAllCalled) { ierr = MeshMoverRegisterAll(PETSC_NULL); CHKERRQ(ierr); } ierr = PetscFListFind(mover->comm, MeshMoverList, method,(void (**)(void)) &r); CHKERRQ(ierr); if (!r) SETERRQ1(PETSC_ERR_ARG_WRONG, "Unknown mesh mover type: %s", method); if (mover->ops->destroy != PETSC_NULL) { ierr = (*mover->ops->destroy)(mover); CHKERRQ(ierr); } ierr = (*r)(mover); CHKERRQ(ierr); ierr = PetscObjectChangeTypeName((PetscObject) mover, method); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverGetType" /*@C MeshMoverGetType - Gets the MeshMover method type name (as a string). Not collective Input Parameter: . mover - The mover Output Parameter: . type - The name of MeshMover method Level: intermediate .keywords: mover, get, type, name .seealso MeshMoverSetType() @*/ int MeshMoverGetType(MeshMover mover, MeshMoverType *type) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mover, MESH_COOKIE); PetscValidPointer(type); if (!MeshMoverRegisterAllCalled) { ierr = MeshMoverRegisterAll(PETSC_NULL); CHKERRQ(ierr); } *type = mover->type_name; PetscFunctionReturn(0); } /*@C MeshMoverRegister - Adds a creation method to the mesh movement package. Synopsis: MeshMoverRegister(char *name, char *path, char *func_name, int (*create_func)(MeshMover)) Not Collective Input Parameters: + name - The name of a new user-defined creation routine . path - The path (either absolute or relative) of the library containing this routine . func_name - The name of the creation routine - create_func - The creation routine itself Notes: MeshMoverRegister() may be called multiple times to add several user-defined meshes. If dynamic libraries are used, then the fourth input argument (create_func) is ignored. Sample usage: .vb MeshMoverRegisterDynamic("my_mesh", "/home/username/my_lib/lib/libO/solaris/libmy.a", "MyMeshMoverCreate", MyMeshMoverCreate); .ve Then, your mesh type can be chosen with the procedural interface via .vb MeshMoverCreate(MPI_Comm, MeshMover *); MeshMoverSetType(vec, "my_mesh_mover") .ve or at runtime via the option .vb -mesh_mover_type my_mesh_mover .ve Note: $PETSC_ARCH and $BOPT occuring in pathname will be replaced with appropriate values. Level: advanced .keywords: MeshMover, register .seealso: MeshMoverRegisterAll(), MeshMoverRegisterDestroy() @*/ #undef __FUNCT__ #define __FUNCT__ "MeshMoverRegister" int MeshMoverRegister(const char sname[], const char path[], const char name[], int (*function)(MeshMover)) { char fullname[256]; int ierr; PetscFunctionBegin; ierr = PetscStrcpy(fullname, path); CHKERRQ(ierr); ierr = PetscStrcat(fullname, ":"); CHKERRQ(ierr); ierr = PetscStrcat(fullname, name); CHKERRQ(ierr); ierr = PetscFListAdd(&MeshMoverList, sname, fullname, (void (*)(void)) function); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshMoverRegisterDestroy" /*@C MeshMoverRegisterDestroy - Frees the list of creation routines for meshes that were registered by MeshMoverRegister(). Not collective Level: advanced .keywords: mesh, register, destroy .seealso: MeshMoverRegister(), MeshMoverRegisterAll() @*/ int MeshMoverRegisterDestroy() { int ierr; PetscFunctionBegin; if (MeshMoverList != PETSC_NULL) { ierr = PetscFListDestroy(&MeshMoverList); CHKERRQ(ierr); MeshMoverList = PETSC_NULL; } MeshMoverRegisterAllCalled = 0; PetscFunctionReturn(0); } EXTERN_C_BEGIN extern int MeshMoverCreate_Triangular_2D(MeshMover); EXTERN_C_END #undef __FUNCT__ #define __FUNCT__ "MeshMoverRegisterAll" /*@C MeshMoverRegisterAll - Registers all of the generation routines in the MeshMover package. Not Collective Input parameter: . path - The dynamic library path Level: advanced .keywords: MeshMover, register, all .seealso: MeshMoverRegister(), MeshMoverRegisterDestroy() @*/ int MeshMoverRegisterAll(const char path[]) { int ierr; PetscFunctionBegin; MeshMoverRegisterAllCalled = 1; ierr = MeshMoverRegister(MESH_MOVER_TRIANGULAR_2D, path, "MeshMoverCreate_Triangular_2D", MeshMoverCreate_Triangular_2D); CHKERRQ(ierr); PetscFunctionReturn(0); } /*-------------------------------------------------- Grid Functions --------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "GridReformMesh" /*@ GridReformMesh - Reconstruct the grid from the current description and underlying mesh boundary. The old mesh is returned if requested. This allows postprocessing of the mesh before GridReform() is called. Collective on Grid Input Parameters: + grid - The grid . refine - Flag indicating whether to refine or recreate the mesh . area - [Optional] A function which gives area constraints for refinement after reform, PETSC_NULL for no refinement - newBd - Flag to determine whether the boundary should be generated (PETSC_TRUE) or taken from storage Output Parameter: . oldMesh - The old mesh Level: advanced .keywords: Grid, mesh, reform, refine .seealso: MeshReform(), GridReform(), GridRefineMesh() @*/ int GridReformMesh(Grid grid, PetscTruth refine, PointFunction area, PetscTruth newBd, Mesh *oldMesh) { Mesh mesh; Mesh newMesh; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); /* Reform mesh */ ierr = GridGetMesh(grid, &mesh); CHKERRQ(ierr); ierr = MeshReform(mesh, refine, area, newBd, &newMesh); CHKERRQ(ierr); /* Save old information */ if (oldMesh != PETSC_NULL){ PetscValidPointer(oldMesh); *oldMesh = mesh; } grid->mesh = newMesh; /* Recalculate grid structures for mesh */ ierr = GridCalcPointBCNodes(grid); CHKERRQ(ierr); /* Recreate class structure */ ierr = FieldClassMapDestroy(grid->cm); CHKERRQ(ierr); ierr = FieldClassMapCreateTriangular2D(grid, grid->numActiveFields, grid->defaultFields, &grid->cm); CHKERRQ(ierr); /* Specific setup */ if (grid->ops->gridreformmesh != PETSC_NULL) { ierr = (*grid->ops->gridreformmesh)(grid); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GridReform" /*@ GridReform - Reconstruct the variable structures, usually after a change in mesh structure. The optional arguments allow projection of the previous solution onto the new mesh. Collective on Grid Input Parameters: + grid - The grid . oldMesh - [Optional] The old mesh - x - [Optional] A grid vector to project onto the new mesh Output Parameter: . y - [Optional] The projected grid vector Level: advanced .keywords: grid, mesh, reform .seealso: GridReformMesh(), MeshReform() @*/ int GridReform(Grid grid, Mesh oldMesh, GVec x, GVec *y) { MPI_Comm comm; VarOrdering order; FieldClassMap map; int numFields; int *fields; PetscTruth isConstrained, explicitConstraints; InterpolationContext iCtx; /* Used to form the rhs */ /* L_2 projection variables */ #if 0 GMat M; /* The (local) mass matrix M */ GMat m; /* The restriction of M to the boundary m */ GMat S; /* The Schur complement S = m^T M^{-1} m */ GVec sol, rhs; /* The solution x and rhs f */ GVec projRhs; /* The projection rhs g */ GVec uzawaRhs; /* The Uzawa rhs m^T M^{-1} f - g */ GVec uzawaSol; /* The Uzawa solution \lambda = (m^T M^{-1} m)^{-1} (m^T M^{-1} f - g) */ KSP ksp; /* The Krylov solver for M */ PC pc; /* The preconditioner for M */ SLES sles; /* The solver for M */ SLES uzawaSles; /* The solver for S */ PetscScalar minusOne = -1.0; PetscScalar zero = 0.0; int its; #endif /* Local Interpolation variables */ VarOrdering reduceOrder, reduceOrderOld; /* TEMPORARY, until I fix the interpolation */ GVec bdReduceVec, prevBdReduceVec; /* TEMPORARY, until I fix the interpolation */ GVec bdReduceVecOld, prevBdReduceVecOld; /* TEMPORARY, until I fix the interpolation */ GVec bdReduceVecDiff, prevBdReduceVecDiff; /* TEMPORARY, until I fix the interpolation */ PetscReal norm; int maxComp; int field, f; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = PetscLogEventBegin(GRID_Reform, grid, 0, 0, 0); CHKERRQ(ierr); /* Make sure the vector is saved in local storage */ ierr = GridGlobalToLocal(grid, INSERT_VALUES, x); CHKERRQ(ierr); /* Need to setup grid again */ grid->setupcalled = PETSC_FALSE; grid->bdSetupCalled = PETSC_FALSE; /* Check arguments */ if ((x == PETSC_NULL) || (y == PETSC_NULL)) { ierr = GridSetupGhostScatter(grid); CHKERRQ(ierr); PetscFunctionReturn(0); } PetscValidHeaderSpecific(x, VEC_COOKIE); PetscValidPointer(y); ierr = PetscObjectGetComm((PetscObject) oldMesh, &comm); CHKERRQ(ierr); ierr = GridIsConstrained(grid, &isConstrained); CHKERRQ(ierr); ierr = GridGetExplicitConstraints(grid, &explicitConstraints); CHKERRQ(ierr); if (explicitConstraints == PETSC_FALSE) { order = grid->order; } else { order = grid->constraintOrder; } ierr = VarOrderingGetClassMap(order, &map); CHKERRQ(ierr); #ifdef INTERFACE_UPDATE ierr = FieldClassMapGetNumFields(map, &numFields); CHKERRQ(ierr); ierr = FieldClassMapGetFields(map, &fields); CHKERRQ(ierr); #endif numFields = map->numFields; fields = map->fields; /* Project fields onto new mesh */ switch(grid->interpolationType) { #if 0 case INTERPOLATION_L2: /* Interpolation Procedure: We must locally interpolate the value of the old field at each Gaussian intergration point in an element of the new field. The function PointFunctionInterpolateField() calls GVecInterpolateField(), which calls DiscretizationInterpolateField(). The vector, field, and old mesh are stored in a context so that PointFunctionInterpolateField() can pass them to GVecInterpolateField(). The indices and values for DiscretizationInterpolateField() are calculated by GridCalcLocalElementVecIndices() and GridLocalToElement(). These functions need the old mesh, classes, variable ordering, and local vector ghostVec in order to correctly calculate the old vector indices and values. However, the GVecEvaluateFunctionGalerkin() uses the new values to get the weak form. The old mesh we have, and the classes are stored in classesOld[]. We will keep the variable ordering and local vector in temporary veriables. The old copies must be set to NULL or they will be deleted. */ /* Save old structures */ iCtx.vec = x; iCtx.mesh = oldMesh; iCtx.order = order; iCtx.ghostVec = grid->ghostVec; /* Setup context */ for(f = 0, maxComp = 1; f < numFields; f++) { maxComp = PetscMax(maxComp, grid->comp[fields[f]]); } ierr = MPI_Comm_size(comm, &iCtx.numProcs); CHKERRQ(ierr); ierr = MPI_Comm_rank(comm, &iCtx.rank); CHKERRQ(ierr); iCtx.batchMode = PETSC_FALSE; iCtx.curPoint = 0; iCtx.maxPoints = iCtx.numProcs; ierr = PetscMalloc((iCtx.numProcs+1) * sizeof(int), &iCtx.numPoints); CHKERRQ(ierr); ierr = PetscMalloc(iCtx.numProcs * sizeof(int), &iCtx.activeProcs); CHKERRQ(ierr); ierr = PetscMalloc(iCtx.numProcs * sizeof(int), &iCtx.calcProcs); CHKERRQ(ierr); ierr = PetscMalloc(iCtx.maxPoints*3 * sizeof(PetscScalar), &iCtx.coords); CHKERRQ(ierr); ierr = PetscMalloc(iCtx.numProcs*maxComp * sizeof(PetscScalar), &iCtx.values); CHKERRQ(ierr); /* Create new variable numbering and ghost structures */ if (explicitConstraints == PETSC_TRUE) { grid->constraintOrder = PETSC_NULL; } else { grid->order = PETSC_NULL; } grid->ghostVec = PETSC_NULL; reduceOrderOld = grid->reduceOrder; prevBdReduceVec = grid->bdReduceVec; prevBdReduceVecOld = grid->bdReduceVecOld; prevBdReduceVecDiff = grid->bdReduceVecDiff; grid->reduceOrder = PETSC_NULL; grid->bdReduceVec = PETSC_NULL; grid->bdReduceVecOld = PETSC_NULL; grid->bdReduceVecDiff = PETSC_NULL; ierr = GridSetupGhostScatter(grid); CHKERRQ(ierr); reduceOrder = grid->reduceOrder; bdReduceVec = grid->bdReduceVec; bdReduceVecOld = grid->bdReduceVecOld; bdReduceVecDiff = grid->bdReduceVecDiff; grid->reduceOrder = reduceOrderOld; grid->bdReduceVec = prevBdReduceVec; grid->bdReduceVecOld = prevBdReduceVecOld; grid->bdReduceVecDiff = prevBdReduceVecDiff; if (grid->bdReduceVecCur == bdReduceVec) { grid->bdReduceVecCur = grid->bdReduceVec; } else if (grid->bdReduceVecCur == bdReduceVecOld) { grid->bdReduceVecCur = grid->bdReduceVecOld; } else if (grid->bdReduceVecCur == bdReduceVecDiff) { grid->bdReduceVecCur = grid->bdReduceVecDiff; } else { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector"); } /* Allocate structures */ if (explicitConstraints) { ierr = GVecCreateConstrained(grid, y); CHKERRQ(ierr); } else { ierr = GVecCreate(grid, y); CHKERRQ(ierr); } ierr = GVecDuplicate(*y, &rhs); CHKERRQ(ierr); ierr = GVecDuplicate(*y, &sol); CHKERRQ(ierr); ierr = GVecCreateBoundaryRestriction(grid, &uzawaSol); CHKERRQ(ierr); ierr = GVecDuplicate(uzawaSol, &projRhs); CHKERRQ(ierr); ierr = GVecDuplicate(uzawaSol, &uzawaRhs); CHKERRQ(ierr); /* Create the rhs */ for(f = 0; f < numFields; f++) { field = fields[f]; iCtx.field = field; if (iCtx.batchMode == PETSC_FALSE) { ierr = GVecEvaluateFunctionGalerkinCollective(rhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx); CHKERRQ(ierr); #if 1 ierr = GVecEvaluateBoundaryFunctionGalerkinCollective(projRhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx); CHKERRQ(ierr); #else ierr = VecSet(&zero, projRhs); #endif } else { ierr = GVecEvaluateFunctionGalerkin(rhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx); CHKERRQ(ierr); ierr = GVecInterpolateFieldBatch(rhs, field, &iCtx); CHKERRQ(ierr); iCtx.curPoint = iCtx.numPoints[iCtx.rank]; ierr = GVecEvaluateFunctionGalerkin(rhs, 1, &field, PointFunctionInterpolateFieldBatch, 1.0, &iCtx); CHKERRQ(ierr); if (iCtx.curPoint != iCtx.numPoints[iCtx.rank+1]) SETERRQ(PETSC_ERR_PLIB, "Invalid batch point processing"); ierr = GVecEvaluateBoundaryFunctionGalerkin(projRhs, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx); CHKERRQ(ierr); ierr = GVecInterpolateFieldBatch(projRhs, field, &iCtx); CHKERRQ(ierr); iCtx.curPoint = iCtx.numPoints[iCtx.rank]; ierr = GVecEvaluateFunctionGalerkin(projRhs, 1, &field, PointFunctionInterpolateFieldBatch, 1.0, &iCtx); CHKERRQ(ierr); if (iCtx.curPoint != iCtx.numPoints[iCtx.rank+1]) SETERRQ(PETSC_ERR_PLIB, "Invalid batch point processing"); } } if (grid->bdReduceVecCur == grid->bdReduceVec) { grid->bdReduceVecCur = bdReduceVec; } else if (grid->bdReduceVecCur == grid->bdReduceVecOld) { grid->bdReduceVecCur = bdReduceVecOld; } else if (grid->bdReduceVecCur == grid->bdReduceVecDiff) { grid->bdReduceVecCur = bdReduceVecDiff; } else { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector"); } grid->reduceOrder = reduceOrder; grid->bdReduceVec = bdReduceVec; grid->bdReduceVecOld = bdReduceVecOld; grid->bdReduceVecDiff = bdReduceVecDiff; if (grid->reduceSystem == PETSC_TRUE) { ierr = VarOrderingDestroy(reduceOrderOld); CHKERRQ(ierr); ierr = VecDestroy(prevBdReduceVec); CHKERRQ(ierr); ierr = VecDestroy(prevBdReduceVecOld); CHKERRQ(ierr); ierr = VecDestroy(prevBdReduceVecDiff); CHKERRQ(ierr); } /* Destroy old structures */ ierr = VarOrderingDestroy(iCtx.order); CHKERRQ(ierr); ierr = VecDestroy(iCtx.ghostVec); CHKERRQ(ierr); /* Cleanup */ ierr = PetscFree(iCtx.numPoints); CHKERRQ(ierr); ierr = PetscFree(iCtx.activeProcs); CHKERRQ(ierr); ierr = PetscFree(iCtx.calcProcs); CHKERRQ(ierr); ierr = PetscFree(iCtx.values); CHKERRQ(ierr); ierr = PetscFree(iCtx.coords); CHKERRQ(ierr); /* Create the inner solver */ ierr = SLESCreate(PETSC_COMM_SELF, &sles); CHKERRQ(ierr); ierr = SLESAppendOptionsPrefix(sles, "grid_proj_"); CHKERRQ(ierr); ierr = SLESGetKSP(sles, &ksp); CHKERRQ(ierr); ierr = SLESGetPC(sles, &pc); CHKERRQ(ierr); ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr); ierr = PCSetType(pc, PCJACOBI); CHKERRQ(ierr); ierr = SLESSetFromOptions(sles); CHKERRQ(ierr); /* Create the outer solver */ ierr = SLESCreate(PETSC_COMM_SELF, &uzawaSles); CHKERRQ(ierr); ierr = SLESAppendOptionsPrefix(uzawaSles, "grid_proj_uzawa_"); CHKERRQ(ierr); ierr = SLESGetKSP(uzawaSles, &ksp); CHKERRQ(ierr); ierr = SLESGetPC(uzawaSles, &pc); CHKERRQ(ierr); ierr = KSPSetType(ksp, KSPGMRES); CHKERRQ(ierr); ierr = PCSetType(pc, PCNONE); CHKERRQ(ierr); ierr = SLESSetFromOptions(uzawaSles); CHKERRQ(ierr); /* Create M, the mass matrix */ ierr = GMatCreate(grid, &M); CHKERRQ(ierr); ierr = GMatEvaluateOperatorGalerkin(M, PETSC_NULL, numFields, fields, fields, IDENTITY, 1.0, MAT_FINAL_ASSEMBLY, PETSC_NULL); CHKERRQ(ierr); ierr = MatNorm(M, NORM_FROBENIUS, &norm); PetscPrintf(PETSC_COMM_SELF, "Mass matrix norm %lf\n", norm); /* Create m, the restriction of M to the boundary */ ierr = GMatCreateBoundaryRestriction(grid, &m); CHKERRQ(ierr); #if 1 ierr = GMatEvaluateBoundaryOperatorGalerkin(m, numFields, grid->cm->fields, grid->bdOrder, grid->bdLocOrder, fields, grid->order, grid->locOrder, IDENTITY, 1.0, MAT_FINAL_ASSEMBLY, PETSC_NULL); CHKERRQ(ierr); #else ierr = MatAssemblyBegin(m, MAT_FINAL_ASSEMBLY); ierr = MatAssemblyEnd(m, MAT_FINAL_ASSEMBLY); ierr = MatZeroEntries(m); #endif ierr = MatNorm(m, NORM_FROBENIUS, &norm); PetscPrintf(PETSC_COMM_SELF, "Restricted mass matrix norm %lf\n", norm); ierr = VecNorm(rhs, NORM_2, &norm); PetscPrintf(PETSC_COMM_SELF, "Rhs norm %lf\n", norm); ierr = VecNorm(projRhs, NORM_2, &norm); PetscPrintf(PETSC_COMM_SELF, "Restricted rhs norm %lf\n", norm); /* Create the Schur complement */ ierr = SLESSetOperators(sles, M, M, SAME_NONZERO_PATTERN); CHKERRQ(ierr); ierr = GMatCreateUzawa(M, m, sles, &S); CHKERRQ(ierr); /* Create the restriction of the rhs */ ierr = SLESSolve(sles, rhs, sol, &its); CHKERRQ(ierr); ierr = MatMultTranspose(m, sol, uzawaRhs); CHKERRQ(ierr); ierr = VecAXPY(&minusOne, projRhs, uzawaRhs); CHKERRQ(ierr); ierr = VecNorm(uzawaRhs, NORM_2, &norm); PetscPrintf(PETSC_COMM_SELF, "Uzawa rhs norm %lf\n", norm); /* Solve system */ ierr = SLESSetOperators(uzawaSles, S, S, SAME_NONZERO_PATTERN); CHKERRQ(ierr); ierr = SLESSolve(uzawaSles, uzawaRhs, uzawaSol, &its); CHKERRQ(ierr); PetscPrintf(PETSC_COMM_SELF, "Uzawa solution took %d iterations\n", its); ierr = VecNorm(uzawaSol, NORM_2, &norm); PetscPrintf(PETSC_COMM_SELF, "Uzawa solution norm %lf\n", norm); /* Get solution */ ierr = MatMult(m, uzawaSol, sol); CHKERRQ(ierr); ierr = VecNorm(sol, NORM_2, &norm); PetscPrintf(PETSC_COMM_SELF, "Boundary solution update norm %lf\n", norm); ierr = VecAYPX(&minusOne, rhs, sol); CHKERRQ(ierr); ierr = VecNorm(sol, NORM_2, &norm); PetscPrintf(PETSC_COMM_SELF, "Solution update norm %lf\n", norm); ierr = SLESSolve(sles, sol, *y, &its); CHKERRQ(ierr); PetscPrintf(PETSC_COMM_SELF, "Projection solution took %d iterations\n", its); ierr = VecNorm(*y, NORM_2, &norm); PetscPrintf(PETSC_COMM_SELF, "Projected solution norm %lf\n", norm); /* Cleanup */ ierr = SLESDestroy(sles); CHKERRQ(ierr); ierr = SLESDestroy(uzawaSles); CHKERRQ(ierr); ierr = GMatDestroyUzawa(S); CHKERRQ(ierr); ierr = GMatDestroy(M); CHKERRQ(ierr); ierr = GMatDestroy(m); CHKERRQ(ierr); ierr = GVecDestroy(rhs); CHKERRQ(ierr); ierr = GVecDestroy(sol); CHKERRQ(ierr); ierr = GVecDestroy(projRhs); CHKERRQ(ierr); ierr = GVecDestroy(uzawaRhs); CHKERRQ(ierr); ierr = GVecDestroy(uzawaSol); CHKERRQ(ierr); break; #endif case INTERPOLATION_LOCAL: /* Interpolation Procedure: We must locally interpolate the value of the old field at each node in an element of the new field. The function PointFunctionInterpolateField()calls GVecInterpolateField(), which calls DiscretizationInterpolateField(). The vector, field, and old mesh are stored in a context so that PointFunctionInterpolateField() can pass them to GVecInterpolateField(). The indices and values for DiscretizationInterpolateField() are calculated by GridCalcLocalElementVecIndices() and GridLocalToElement(). These functions need the old mesh, classes, variable ordering, and local vector ghostVec in order to correctly calculate the old vector indices and values. The old mesh we have, and the classes are stored in classesOld[]. We will keep the variable ordering and local vector in temporary veriables. The old copies must be set to NULL or they will be deleted. */ /* Save old structures */ iCtx.vec = x; iCtx.mesh = oldMesh; iCtx.order = order; iCtx.ghostVec = grid->ghostVec; /* Setup context */ for(f = 0, maxComp = 1; f < numFields; f++) { maxComp = PetscMax(maxComp, grid->fields[fields[f]].numComp); } ierr = MPI_Comm_size(comm, &iCtx.numProcs); CHKERRQ(ierr); ierr = MPI_Comm_rank(comm, &iCtx.rank); CHKERRQ(ierr); iCtx.batchMode = PETSC_FALSE; iCtx.curPoint = 0; iCtx.maxPoints = iCtx.numProcs; ierr = PetscMalloc((iCtx.numProcs+1) * sizeof(int), &iCtx.numPoints); CHKERRQ(ierr); ierr = PetscMalloc(iCtx.numProcs * sizeof(int), &iCtx.activeProcs); CHKERRQ(ierr); ierr = PetscMalloc(iCtx.numProcs * sizeof(int), &iCtx.calcProcs); CHKERRQ(ierr); ierr = PetscMalloc(iCtx.maxPoints*3 * sizeof(PetscScalar), &iCtx.coords); CHKERRQ(ierr); ierr = PetscMalloc(iCtx.numProcs*maxComp * sizeof(PetscScalar), &iCtx.values); CHKERRQ(ierr); ierr = PetscMemzero(iCtx.numPoints, (iCtx.numProcs+1) * sizeof(int)); CHKERRQ(ierr); ierr = PetscOptionsHasName(PETSC_NULL, "-grid_int_mode", &iCtx.batchMode); CHKERRQ(ierr); /* Create new variable numbering and ghost structures */ if (explicitConstraints == PETSC_TRUE) { grid->constraintOrder = PETSC_NULL; } else { grid->order = PETSC_NULL; } grid->ghostVec = PETSC_NULL; reduceOrderOld = grid->reduceOrder; prevBdReduceVec = grid->bdReduceVec; prevBdReduceVecOld = grid->bdReduceVecOld; prevBdReduceVecDiff = grid->bdReduceVecDiff; grid->reduceOrder = PETSC_NULL; grid->bdReduceVec = PETSC_NULL; grid->bdReduceVecOld = PETSC_NULL; grid->bdReduceVecDiff = PETSC_NULL; ierr = GridSetupGhostScatter(grid); CHKERRQ(ierr); reduceOrder = grid->reduceOrder; bdReduceVec = grid->bdReduceVec; bdReduceVecOld = grid->bdReduceVecOld; bdReduceVecDiff = grid->bdReduceVecDiff; grid->reduceOrder = reduceOrderOld; grid->bdReduceVec = prevBdReduceVec; grid->bdReduceVecOld = prevBdReduceVecOld; grid->bdReduceVecDiff = prevBdReduceVecDiff; if (grid->bdReduceVecCur == bdReduceVec) { grid->bdReduceVecCur = grid->bdReduceVec; } else if (grid->bdReduceVecCur == bdReduceVecOld) { grid->bdReduceVecCur = grid->bdReduceVecOld; } else if (grid->bdReduceVecCur == bdReduceVecDiff) { grid->bdReduceVecCur = grid->bdReduceVecDiff; } else { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector"); } /* Allocate structures */ if (explicitConstraints) { ierr = GVecCreateConstrained(grid, y); CHKERRQ(ierr); } else { ierr = GVecCreate(grid, y); CHKERRQ(ierr); } /* Locally interpolate the solution values */ for(f = 0; f < numFields; f++) { field = fields[f]; iCtx.field = field; if (iCtx.batchMode == PETSC_FALSE) { ierr = GVecEvaluateFunctionCollective(*y, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx); CHKERRQ(ierr); } else { ierr = GVecEvaluateFunction(*y, 1, &field, PointFunctionInterpolateField, 1.0, &iCtx); CHKERRQ(ierr); ierr = GVecInterpolateFieldBatch(*y, field, &iCtx); CHKERRQ(ierr); iCtx.curPoint = iCtx.numPoints[iCtx.rank]; ierr = GVecEvaluateFunction(*y, 1, &field, PointFunctionInterpolateFieldBatch, 1.0, &iCtx); CHKERRQ(ierr); if (iCtx.curPoint != iCtx.numPoints[iCtx.rank+1]) SETERRQ(PETSC_ERR_PLIB, "Invalid batch point processing"); } } if (grid->bdReduceVecCur == grid->bdReduceVec) { grid->bdReduceVecCur = bdReduceVec; } else if (grid->bdReduceVecCur == grid->bdReduceVecOld) { grid->bdReduceVecCur = bdReduceVecOld; } else if (grid->bdReduceVecCur == grid->bdReduceVecDiff) { grid->bdReduceVecCur = bdReduceVecDiff; } else { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid reduction vector"); } grid->reduceOrder = reduceOrder; grid->bdReduceVec = bdReduceVec; grid->bdReduceVecOld = bdReduceVecOld; grid->bdReduceVecDiff = bdReduceVecDiff; if (grid->reduceSystem == PETSC_TRUE) { ierr = VarOrderingDestroy(reduceOrderOld); CHKERRQ(ierr); ierr = VecDestroy(prevBdReduceVec); CHKERRQ(ierr); ierr = VecDestroy(prevBdReduceVecOld); CHKERRQ(ierr); ierr = VecDestroy(prevBdReduceVecDiff); CHKERRQ(ierr); } ierr = VecNorm(*y, NORM_2, &norm); PetscPrintf(PETSC_COMM_SELF, "Projected solution norm %g\n", norm); /* Destroy old structures */ ierr = VarOrderingDestroy(iCtx.order); CHKERRQ(ierr); ierr = VecDestroy(iCtx.ghostVec); CHKERRQ(ierr); /* Cleanup */ ierr = PetscFree(iCtx.numPoints); CHKERRQ(ierr); ierr = PetscFree(iCtx.activeProcs); CHKERRQ(ierr); ierr = PetscFree(iCtx.calcProcs); CHKERRQ(ierr); ierr = PetscFree(iCtx.values); CHKERRQ(ierr); ierr = PetscFree(iCtx.coords); CHKERRQ(ierr); break; case INTERPOLATION_HALO: case INTERPOLATION_L2_LOCAL: SETERRQ1(PETSC_ERR_SUP, "Interpolation type %d not yet supported", grid->interpolationType); default: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid interpolation type %d", grid->interpolationType); } /* Anything specific to the implementation */ if (grid->ops->gridreform != PETSC_NULL) { ierr = (*grid->ops->gridreform)(grid, oldMesh, x, y); CHKERRQ(ierr); } ierr = PetscLogEventEnd(GRID_Reform, grid, 0, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); }