SizeField.cpp

Go to the documentation of this file.

#include "MeshAdaptPUMI.h"
#include <apf.h>
#include <apfVector.h>
#include <apfMesh.h>
#include <apfShape.h>
#include <apfDynamicVector.h>
#include <apfCavityOp.h>
#include <string>
#include <iostream>
#include <sstream>
#include <PCU.h>
#include <samElementCount.h>
#include <queue>
#include <algorithm> 
 
static void SmoothField(apf::Field *f);
void gradeAnisoMesh(apf::Mesh* m,double gradingFactor,PCU_t PCUObj);
void gradeAspectRatio(apf::Mesh* m, int idx, double gradingFactor,PCU_t PCUObj);
 
/* Based on the distance from the interface epsilon can be controlled to determine
   thickness of refinement near the interface */
static double isotropicFormula(double phi, double dphi, double verr, double hmin, double hmax, double phi_s = 0, double epsFact = 0)
{
  double size;
  double dphi_size_factor;
  double v_size_factor;
  if (fabs(phi_s) < (epsFact*7.5) * hmin)
    return hmin;
  else
    return hmax;
}
 
static void setSizeField(apf::Mesh2 *m,apf::MeshEntity *vertex,double h,apf::MeshTag *marker,apf::Field* sizeField,PCU_t PCUObj)
//helper function for banded interface to facilitate with setting the proper mesh size and parallel communication
{
  int isMarked=0;
  if(m->hasTag(vertex,marker))
    isMarked=1;
  double h_new;
  if(isMarked)
    h_new = std::min(h,apf::getScalar(sizeField,vertex,0));
  else
  {
    h_new = h;
    int newMark = 1;
    m->setIntTag(vertex,marker,&newMark);
  }
  apf::setScalar(sizeField,vertex,0,h_new);
 
  //Parallel Communication with owning copy
  if(!m->isOwned(vertex))
  {
    apf::Copies remotes;
    m->getRemotes(vertex,remotes);
    int owningPart=m->getOwner(vertex);
    PCU_COMM_PACK(PCUObj, owningPart, remotes[owningPart]);
    PCU_COMM_PACK(PCUObj, owningPart, h_new);
  }
}
 
int MeshAdaptPUMIDrvr::setSphereSizeField()
{
  freeField(size_iso);
  size_iso = apf::createLagrangeField(m, "proteus_size", apf::SCALAR, 1);
 
  apf::MeshIterator *it = m->begin(0);
  apf::MeshEntity* ent;
  while ((ent = m->iterate(it)))
  {
    int modelTag = m->getModelTag(m->toModel(ent));
    //std::cout<<"This is the model tag "<<modelTag<<std::endl;
    double sizeDesired;
    if(modelTag==123)
        sizeDesired=hmin;
    else
        sizeDesired=hmax;
    apf::setScalar(size_iso,ent,0,sizeDesired);
  }
  m->end(it);
  gradeMesh(1.5);
}
 
 
int MeshAdaptPUMIDrvr::calculateSizeField(double L_band)
//Implementation of banded interface, edge intersection algorithm
//If mesh edge intersects the 0 level-set, then the adjacent edges need to be refined 
{
  if(m->findField("interfaceBand"))
    apf::destroyField(m->findField("interfaceBand"));
 
  apf::Field* interfaceBand = apf::createLagrangeField(m, "interfaceBand", apf::SCALAR, 1);
  apf::Field *phif = m->findField("phi");
  assert(phif);
 
  apf::MeshTag* vertexMarker = m->createIntTag("vertexMarker",1);
  apf::MeshIterator *it = m->begin(1);
  apf::MeshEntity *edge;
 
  //double L_band = (numAdaptSteps+N_interface_band)*hPhi;
 
  PCU_Comm_Begin(PCUObj);
  while ((edge = m->iterate(it)))
  {
    apf::Adjacent edge_adjVerts;
    m->getAdjacent(edge,0,edge_adjVerts);
    apf::MeshEntity *vertex1 = edge_adjVerts[0];
    apf::MeshEntity *vertex2 = edge_adjVerts[1];
    double phi1 = apf::getScalar(phif,vertex1,0);
    double phi2 = apf::getScalar(phif,vertex2,0);
    int caseNumber = 1;
    if(std::fabs(phi1)>L_band)
      caseNumber++;
    if(std::fabs(phi2)>L_band)
      caseNumber++;
 
    if(caseNumber==1 || caseNumber == 2)
    {
      setSizeField(m,vertex1,hPhi,vertexMarker,interfaceBand,PCUObj);
      setSizeField(m,vertex2,hPhi,vertexMarker,interfaceBand,PCUObj);
    }
    else
    {
      if (phi1*phi2 <0)
      {
        setSizeField(m,vertex1,hPhi,vertexMarker,interfaceBand,PCUObj);
        setSizeField(m,vertex2,hPhi,vertexMarker,interfaceBand,PCUObj);
      }
      else
      {
        setSizeField(m,vertex1,hmax,vertexMarker,interfaceBand,PCUObj);
        setSizeField(m,vertex2,hmax,vertexMarker,interfaceBand,PCUObj);
      }
    }
 
  }//end while
 
  PCU_Comm_Send(PCUObj);
 
  //Take minimum between received value and current value
  apf::MeshEntity *ent;
  double h_received;
  while(PCU_Comm_Receive(PCUObj))
  {
    //Note: the only receiving entities should be owning copies
    PCU_COMM_UNPACK(PCUObj, ent);
    PCU_COMM_UNPACK(PCUObj, h_received);
    //take minimum of received values
    double h_current = apf::getScalar(interfaceBand,ent,0);
    double h_final = std::min(h_current,h_received);
    apf::setScalar(interfaceBand,ent,0,h_final);
  }
 
  //Synchronization has all remote copies track the owning copy value
  apf::synchronize(interfaceBand);
  m->end(it);
 
  m->destroyTag(vertexMarker);
 
  //add to queue
  sizeFieldList.push(interfaceBand);
  return 0;
}
 
int intersectsInterface(apf::MeshEntity* edge, apf::Field* levelSet)
{
    apf::Mesh* m = apf::getMesh(levelSet);
    apf::Adjacent edge_adjVerts;
    m->getAdjacent(edge,0,edge_adjVerts);
    apf::MeshEntity *vertex1 = edge_adjVerts[0];
    apf::MeshEntity *vertex2 = edge_adjVerts[1];
    double phi1 = apf::getScalar(levelSet,vertex1,0);
    double phi2 = apf::getScalar(levelSet,vertex2,0);
    int doesIntersect = 0;
    if(phi1*phi2 < 0) //implies different signs and therefore intersects interface
        doesIntersect = 1;
    return doesIntersect;
}
 
//Struct definition
struct edgeWalkerInfo{
    apf::MeshEntity* vertex;  
    apf::Vector3 actualPosition;
    double direction;
    double L_local;
    apf::MeshTag* trackerTag;
    //const char* tagName;
    int edgeID;
    int initialRank;
};
 
int checkForPropagation(apf::Mesh* m, edgeWalkerInfo inputObject)
{
    apf::MeshEntity* vert = inputObject.vertex;
    apf::Vector3 actualPosition = inputObject.actualPosition;
    double L_local = inputObject.L_local;
    double direction = inputObject.direction;
 
    apf::MeshTag* vertexMaxTraverse = m->findTag("maximumTraversal");
    apf::Field* predictInterfaceBand = m->findField("predictInterfaceBand");
    apf::Field* levelSet = m->findField("phi");
 
    apf::Vector3 pt_vert;
    m->getPoint(vert,0,pt_vert);
    apf::Vector3 difference_vect = pt_vert-actualPosition;
 
    //check if vertex needs to be added to queue based on traversal distance
    int dontContinue = 0;
    if(m->hasTag(vert,vertexMaxTraverse))
    {
        double traversalDistance;
        m->getDoubleTag(vert,vertexMaxTraverse,&traversalDistance);    
        if((L_local-difference_vect.getLength()) < traversalDistance*1.01) //has to be 1% higher
            dontContinue=1;
    } 
 
    //directionality
    double phiCurrent = apf::getScalar(levelSet,vert,0);
    if((difference_vect.getLength() > L_local) || dontContinue || (phiCurrent*direction<=0))
        return 0;
    else
        return 1;
}
 
int BFS_propagation(apf::Mesh* m, std::queue<edgeWalkerInfo> &markedVertices, PCU_t PCUObj)
{   
//objects in the queue are assumed to be checked and needing to be modified
//the adjacencies are checked before adding into the queue
//for parallelism, what's important is to stop a constant back-and-forth communication on shared vertices; we do this by considering communication when we add an entry into the queue
 
    //get the latest object
    edgeWalkerInfo inputObject = markedVertices.front();
    markedVertices.pop();
 
    //set the variables from the inputObject
    apf::MeshEntity* vert = inputObject.vertex;
    apf::Vector3 actualPosition = inputObject.actualPosition;
    double L_local = inputObject.L_local;
    double direction = inputObject.direction;
 
    //get necessary fields
    apf::MeshTag* vertexMaxTraverse = m->findTag("maximumTraversal");
    apf::Field* predictInterfaceBand = m->findField("predictInterfaceBand");
    apf::Field* levelSet = m->findField("phi");
 
    int needsParallel=0;
    
    apf::Adjacent vertex_adjVerts; 
    apf::getBridgeAdjacent(m,vert,1,0,vertex_adjVerts);
    for(int i=0;i<vertex_adjVerts.getSize();i++)
    {
 
        apf::MeshEntity* newVert = vertex_adjVerts[i];
        inputObject.vertex=newVert;
 
        if(checkForPropagation(m,inputObject))
        {
 
            //set new traversal distance
            apf::Vector3 pt_vert;
            m->getPoint(newVert,0,pt_vert);
            apf::Vector3 difference_vect = pt_vert-actualPosition;
 
            double traversalDistance = L_local-difference_vect.getLength();
            m->setDoubleTag(newVert,vertexMaxTraverse,&traversalDistance);
 
            //set size
            apf::setScalar(predictInterfaceBand,newVert,0,apf::getScalar(predictInterfaceBand,vert,0));
            //add vertex to queue 
            inputObject.vertex = newVert;
            markedVertices.push(inputObject);
 
            if(m->isShared(newVert))
            {
                int initialRank = PCU_Comm_Self(PCUObj);
                double desiredSize = apf::getScalar(predictInterfaceBand,vert,0);
                apf::Copies remotes;
                m->getRemotes(newVert,remotes);
                for(apf::Copies::iterator iter=remotes.begin(); iter!=remotes.end();++iter)
                {
                    PCU_COMM_PACK(PCUObj, iter->first, iter->second);
                    PCU_COMM_PACK(PCUObj, iter->first, L_local);
                    PCU_COMM_PACK(PCUObj, iter->first, actualPosition);
                    PCU_COMM_PACK(PCUObj, iter->first, direction);
                    PCU_COMM_PACK(PCUObj, iter->first, initialRank);
                    PCU_COMM_PACK(PCUObj, iter->first, inputObject.edgeID);
                    PCU_COMM_PACK(PCUObj, iter->first, desiredSize);
                }
                needsParallel++;
            }
        }
    } //end of adjacent
 
    return needsParallel;
}
 
 
void MeshAdaptPUMIDrvr::predictiveInterfacePropagation()
//compute Lband
//edge walk
{
    apf::Field* interfaceBand = m->findField("interfaceBand");
    apf::Field* velocity = m->findField("velocity");
    apf::Field* levelSet = m->findField("phi");
 
    //get gradient field
    apf::Field *gradphi = apf::recoverGradientByVolume(levelSet);
 
    apf::Field* predictInterfaceBand = apf::createLagrangeField(m,"predictInterfaceBand",apf::SCALAR,1);
    apf::copyData(predictInterfaceBand,interfaceBand);
 
    //edge-walk to predict
    apf::MeshTag* vertexMaxTraverse = m->createDoubleTag("maximumTraversal",1); //define tag field for each vertex to store maximum distance that will be travelled
 
    apf::MeshEntity* edge;
    apf::MeshIterator* it = m->begin(1);
    
    std::queue <edgeWalkerInfo> markedVertices;
 
    PCU_Comm_Begin(PCUObj);
    while( (edge = m->iterate(it)) )
    {
        if( intersectsInterface(edge,levelSet))
        {
            //get the parameterized position of interface along edge
            apf::Adjacent edge_adjVerts;
            m->getAdjacent(edge,0,edge_adjVerts);
            apf::MeshEntity *vertex1 = edge_adjVerts[0];
            apf::MeshEntity *vertex2 = edge_adjVerts[1];
            double phi1 = apf::getScalar(levelSet,vertex1,0);
            double phi2 = apf::getScalar(levelSet,vertex2,0);
            apf::Vector3 pt_1, pt_2;
            m->getPoint(vertex1,0,pt_1);
            m->getPoint(vertex2,0,pt_2);
            apf::Vector3 edgeVector = pt_2-pt_1;
            double edgeLength = apf::measure(m,edge);
            double zeroPosition =  2*(-phi1/(phi2-phi1))-1.0; //parametric position of interface on the edge
            double relativePosition = -phi1/(phi2-phi1); //same as zeroPosition but in interval of [0,1]
            apf::Vector3 actualPosition = (pt_2-pt_1)*relativePosition + pt_1;           
 
            apf::Vector3 edgePoint(zeroPosition,0.0,0.0);
            apf::Element* phiElem = apf::createElement(levelSet,edge);
            apf::Element* gradPhiElem = apf::createElement(gradphi,edge);
            apf::Element* velocityElem = apf::createElement(velocity,edge);
            apf::Vector3 localVelocity;
            apf::getVector(velocityElem,edgePoint,localVelocity);
            apf::Vector3 localInterfaceNormal;
            apf::getVector(gradPhiElem,edgePoint,localInterfaceNormal);
            apf::destroyElement(phiElem);
            apf::destroyElement(velocityElem);
            apf::destroyElement(gradPhiElem);
 
            //get L_local 
            double L_local = localVelocity.getLength()*numAdaptSteps*delta_T;
        
            L_local += (N_interface_band)*hPhi; //add blending region   
 
            //get direction, multiply this with levelSet value to determine if in same direction
            double signValue = localVelocity*localInterfaceNormal;
 
            //find adjacent vertices and their adjacent edges
            for(int i=0; i<edge_adjVerts.getSize();i++)
            {
                edgeWalkerInfo inputObject;
                inputObject.vertex = edge_adjVerts[i];
                inputObject.actualPosition = actualPosition;
                inputObject.direction = signValue;
                inputObject.L_local = L_local;
                inputObject.edgeID = localNumber(edge);
                inputObject.initialRank = PCU_Comm_Self(PCUObj);
                if(checkForPropagation(m,inputObject))
                {
                    markedVertices.push(inputObject);
 
                    //check for parallel
                    if(m->isShared(inputObject.vertex))
                    {
                        int initialRank = PCU_Comm_Self(PCUObj);
                        double desiredSize = apf::getScalar(predictInterfaceBand,inputObject.vertex,0);
                        apf::Copies remotes;
                        m->getRemotes(inputObject.vertex,remotes);
                        for(apf::Copies::iterator iter=remotes.begin(); iter!=remotes.end();++iter)
                        {
                            PCU_COMM_PACK(PCUObj, iter->first, iter->second);
                            PCU_COMM_PACK(PCUObj, iter->first, L_local);
                            PCU_COMM_PACK(PCUObj, iter->first, actualPosition);
                            PCU_COMM_PACK(PCUObj, iter->first, inputObject.direction);
                            PCU_COMM_PACK(PCUObj, iter->first, initialRank);
                            PCU_COMM_PACK(PCUObj, iter->first, inputObject.edgeID);
                            PCU_COMM_PACK(PCUObj, iter->first, desiredSize);
                        }
                    }
 
                }
            }
        } //end if interface edge
    }
    m->end(it);
 
    //The following parallel code is just for the initialization step
    //There might be a way to put all of this under a function as this is repeated code later on
    PCU_Comm_Send(PCUObj);
    while(PCU_Comm_Receive(PCUObj))
    {
        apf::MeshEntity* vertex;
        double L_local;
        apf::Vector3 actualPosition;
        double direction;
        int initialRank;
        int edgeID;
        double desiredSize;
        PCU_COMM_UNPACK(PCUObj, vertex);
        PCU_COMM_UNPACK(PCUObj, L_local);
        PCU_COMM_UNPACK(PCUObj, actualPosition);
        PCU_COMM_UNPACK(PCUObj, direction);
        PCU_COMM_UNPACK(PCUObj, initialRank);
        PCU_COMM_UNPACK(PCUObj, edgeID);
        PCU_COMM_UNPACK(PCUObj, desiredSize);
            
        edgeWalkerInfo inputObject;
        inputObject.vertex = vertex;
        inputObject.L_local = L_local;
        inputObject.actualPosition = actualPosition;
        inputObject.direction = direction;
        inputObject.edgeID = edgeID;
        inputObject.initialRank=initialRank;
 
        //ensures that the size is the same across parts
        if(desiredSize < apf::getScalar(predictInterfaceBand,vertex,0))
            apf::setScalar(predictInterfaceBand,vertex,0,desiredSize);
 
        if(checkForPropagation(m,inputObject))
        {
            apf::Vector3 pt_vert;
            m->getPoint(vertex,0,pt_vert);
            apf::Vector3 difference_vect = pt_vert-actualPosition;
 
            double traversalDistance = L_local-difference_vect.getLength();
            m->setDoubleTag(vertex,vertexMaxTraverse,&traversalDistance);
 
            markedVertices.push(inputObject);
        }
    }
 
    //Parallel preparations
    int needsParallel=1;
 
    while(needsParallel>0)
    {
        needsParallel=0;
 
        PCU_Comm_Begin(PCUObj);
 
        //Handle the queue
        while(!markedVertices.empty())
        {
            needsParallel+=BFS_propagation(m,markedVertices,PCUObj);
        }
        PCU_Add_Ints(PCUObj, &needsParallel,1);
 
        PCU_Comm_Send(PCUObj);
        while( PCU_Comm_Receive(PCUObj) )
        {
            apf::MeshEntity* vertex;
            double L_local;
            apf::Vector3 actualPosition;
            double direction;
            int initialRank;
            int edgeID;
            double desiredSize;
            PCU_COMM_UNPACK(PCUObj, vertex);
            PCU_COMM_UNPACK(PCUObj, L_local);
            PCU_COMM_UNPACK(PCUObj, actualPosition);
            PCU_COMM_UNPACK(PCUObj, direction);
            PCU_COMM_UNPACK(PCUObj, initialRank);
            PCU_COMM_UNPACK(PCUObj, edgeID);
            PCU_COMM_UNPACK(PCUObj, desiredSize);
            
 
            edgeWalkerInfo inputObject;
            inputObject.vertex = vertex;
            inputObject.L_local = L_local;
            inputObject.actualPosition = actualPosition;
            inputObject.direction = direction;
            inputObject.edgeID = edgeID;
            inputObject.initialRank=initialRank;
 
            //ensures that the size is the same across parts
            if(desiredSize < apf::getScalar(predictInterfaceBand,vertex,0))
                apf::setScalar(predictInterfaceBand,vertex,0,desiredSize);
 
            if(checkForPropagation(m,inputObject))
            {
                apf::Vector3 pt_vert;
                m->getPoint(vertex,0,pt_vert);
                apf::Vector3 difference_vect = pt_vert-actualPosition;
 
                double traversalDistance = L_local-difference_vect.getLength();
                m->setDoubleTag(vertex,vertexMaxTraverse,&traversalDistance);
 
                markedVertices.push(inputObject);
            }
        }
    }
 
    apf::copyData(interfaceBand,predictInterfaceBand);
    apf::destroyField(predictInterfaceBand);
    apf::destroyField(gradphi);
 
    //get all tags for destruction
    apf::MeshIterator* tagIt = m->begin(0);
    apf::MeshEntity* taggedVertex;
    while( taggedVertex = m->iterate(tagIt) )
    {
        if(m->hasTag(taggedVertex,vertexMaxTraverse))
            m->removeTag(taggedVertex,vertexMaxTraverse);
    }
    m->end(tagIt);
    m->destroyTag(vertexMaxTraverse);
 
}
 
void MeshAdaptPUMIDrvr::isotropicIntersect()
{
  freeField(size_iso);
  if(m->findField("proteus_size"))
    apf::destroyField(m->findField("proteus_size"));
 
  size_iso = apf::createFieldOn(m, "proteus_size", apf::SCALAR);
 
  apf::MeshEntity *vert;
  apf::MeshIterator *it = m->begin(0);
  
  apf::Field *field = sizeFieldList.front();
  apf::copyData(size_iso,field);
  sizeFieldList.pop();
  //apf::destroyField(field);
  while(!sizeFieldList.empty())
  {
    field = sizeFieldList.front();
    while(vert = m->iterate(it))
    {
      double value1 = apf::getScalar(size_iso,vert,0);
      double value2 = apf::getScalar(field,vert,0);
      double minValue = std::min(value1,value2);
      apf::setScalar(size_iso,vert,0,minValue);
    } 
    sizeFieldList.pop();
    //apf::destroyField(field);
  }
  if(nAdapt < 2) //if first few adapts, need to make sure that gradation is low to allow interface to develop properly
    gradeMesh(1.1);
  else
    gradeMesh(gradingFactor);
}
 
//taken from Dan's superconvergent patch recovery code
void MeshAdaptPUMIDrvr::averageToEntity(apf::Field *ef, apf::Field *vf,
                                        apf::MeshEntity *ent)
//Function used to convert a region/element field into a vertex field via averaging.
//For each vertex, the average value of the adjacent elements is taken
//Input:
//  ef is the element field
//  ent is the target vertex
//Output:
//  vf is the vertex field
{
  apf::Mesh *m = apf::getMesh(ef);
  apf::Adjacent elements;
  m->getAdjacent(ent, m->getDimension(), elements);
  double s = 0;
  for (std::size_t i = 0; i < elements.getSize(); ++i)
    s += apf::getScalar(ef, elements[i], 0);
  s /= elements.getSize();
  apf::setScalar(vf, ent, 0, s);
  return;
}
 
void MeshAdaptPUMIDrvr::minToEntity(apf::Field* ef, apf::Field* vf,
    apf::MeshEntity* ent)
{
  apf::Mesh* m = apf::getMesh(ef);
  apf::Adjacent elements;
  m->getAdjacent(ent, m->getDimension(), elements);
  double s=0;
  for (std::size_t i=0; i < elements.getSize(); ++i){
    if(i==0)
      s = apf::getScalar(ef, elements[i], 0);
    else if(apf::getScalar(ef,elements[i],0) < s)
      s= apf::getScalar(ef,elements[i],0);
  }
  apf::setScalar(vf, ent, 0, s);
  return;
}
 
 
void MeshAdaptPUMIDrvr::volumeAverageToEntity(apf::Field *ef, apf::Field *vf,
                                              apf::MeshEntity *ent)
//Serves the same purpose as averageToEntity but considers a volume-weighted average
{
  apf::Mesh *m = apf::getMesh(ef);
  apf::Adjacent elements;
  apf::MeshElement *testElement;
  m->getAdjacent(ent, m->getDimension(), elements);
  double s = 0;
  double VolumeTotal = 0;
  for (std::size_t i = 0; i < elements.getSize(); ++i)
  {
    testElement = apf::createMeshElement(m, elements[i]);
    s += apf::getScalar(ef, elements[i], 0)*apf::measure(testElement);
    VolumeTotal += apf::measure(testElement);
    apf::destroyMeshElement(testElement);
  }
  s /= VolumeTotal;
  apf::setScalar(vf, ent, 0, s);
  return;
}
 
void errorAverageToEntity(apf::Field *ef, apf::Field *vf, apf::Field* err, apf::MeshEntity *ent)
//Serves the same purpose as averageToEntity but considers a error-weighted average
{
  apf::Mesh *m = apf::getMesh(ef);
  apf::Adjacent elements;
  m->getAdjacent(ent, m->getDimension(), elements);
  double s = 0;
  double errorTotal = 0;
  for (std::size_t i = 0; i < elements.getSize(); ++i)
  {
    s += apf::getScalar(ef, elements[i], 0)*apf::getScalar(err,elements[i],0);
    errorTotal += apf::getScalar(err,elements[i],0);
  }
  s /= errorTotal;
/*
  if (comm_rank == 0)
  {
    std::cout << "What is s final? " << s << std::endl;
  }
*/
  apf::setScalar(vf, ent, 0, s);
  return;
}
 
 
static apf::Field *extractSpeed(apf::Field *velocity)
//Function used to convert the velocity field into a speed field
{
  apf::Mesh *m = apf::getMesh(velocity);
  apf::Field *speedF = apf::createLagrangeField(m, "proteus_speed", apf::SCALAR, 1);
  apf::MeshIterator *it = m->begin(0);
  apf::MeshEntity *v;
  apf::Vector3 vel_vect;
  while ((v = m->iterate(it)))
  {
    apf::getVector(velocity, v, 0, vel_vect);
    double speed = vel_vect.getLength();
    apf::setScalar(speedF, v, 0, speed);
  }
  m->end(it);
  return speedF;
}
 
static apf::Matrix3x3 hessianFormula(apf::Matrix3x3 const &g2phi)
//Function used to output a symmetric hessian
{
  apf::Matrix3x3 g2phit = apf::transpose(g2phi);
  return (g2phi + g2phit) / 2;
}
 
static apf::Field *computeHessianField(apf::Field *grad2phi)
//Function that writes a hessian field
{
  apf::Mesh *m = apf::getMesh(grad2phi);
  apf::Field *hessf = createLagrangeField(m, "proteus_hess", apf::MATRIX, 1);
  apf::MeshIterator *it = m->begin(0);
  apf::MeshEntity *v;
  while ((v = m->iterate(it)))
  {
    apf::Matrix3x3 g2phi;
    apf::getMatrix(grad2phi, v, 0, g2phi);
    apf::Matrix3x3 hess = hessianFormula(g2phi);
    apf::setMatrix(hessf, v, 0, hess);
  }
  m->end(it);
  return hessf;
}
 
static apf::Field *computeMetricField(apf::Field *gradphi, apf::Field *grad2phi, apf::Field *size_iso, double eps_u)
//Function used to generate a field of metric tensors at vertices. It is meant to be an umbrella function that can compute Hessians too.
//Currently needs development.
{
  apf::Mesh *m = apf::getMesh(grad2phi);
  apf::Field *metricf = createLagrangeField(m, "proteus_metric", apf::MATRIX, 1);
  apf::MeshIterator *it = m->begin(0);
  apf::MeshEntity *v;
  while ((v = m->iterate(it)))
  {
    apf::Matrix3x3 g2phi;
    apf::getMatrix(grad2phi, v, 0, g2phi);
/*
    apf::Vector3 gphi;
    apf::getVector(gradphi, v, 0, gphi);
    apf::Matrix3x3 gphigphit(gphi[0] * gphi[0], gphi[0] * gphi[1], gphi[0] * gphi[2],
                             gphi[0] * gphi[1], gphi[1] * gphi[1], gphi[1] * gphi[2],
                             gphi[0] * gphi[2], gphi[1] * gphi[2], gphi[2] * gphi[2]);
*/
    apf::Matrix3x3 hess = hessianFormula(g2phi);
    apf::Matrix3x3 metric = hess;
    //apf::Matrix3x3 metric = gphigphit/(apf::getScalar(size_iso,v,0)*apf::getScalar(size_iso,v,0))+ hess/eps_u;
    //apf::Matrix3x3 metric(1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0);
    apf::setMatrix(metricf, v, 0, metric);
  }
  m->end(it);
  return metricf;
}
 
// Gaussian, Mean and principal curvatures based on Hessian and gradient of phi.
static void curveFormula(apf::Matrix3x3 const &h, apf::Vector3 const &g,
                         apf::Vector3 &curve)
{
  double a = (h[1][1] + h[2][2]) * g[0] * g[0] + (h[0][0] + h[2][2]) * g[1] * g[1] + (h[0][0] + h[1][1]) * g[2] * g[2];
 
  double b = g[0] * g[1] * h[0][1] + g[0] * g[2] * h[0][2] + g[1] * g[2] * h[1][2];
 
  double Km = (a - 2 * b) / pow(g * g, 1.5);
 
  double c = g[0] * g[0] * (h[1][1] * h[2][2] - h[1][2] * h[1][2]) + g[1] * g[1] * (h[0][0] * h[2][2] - h[0][2] * h[0][2]) + g[2] * g[2] * (h[0][0] * h[1][1] - h[0][1] * h[0][1]);
 
  double d = g[0] * g[1] * (h[0][2] * h[1][2] - h[0][1] * h[2][2]) + g[1] * g[2] * (h[0][1] * h[0][2] - h[1][2] * h[0][0]) + g[0] * g[2] * (h[0][1] * h[1][2] - h[0][2] * h[1][1]);
 
  double Kg = (c + 2 * d) / pow(g * g, 2);
 
  curve[0] = Km;                      //Mean curvature= (k_1+k_2)/2, Not used in normal direction
  curve[1] = Km + sqrt(Km * Km - Kg); //k_1, First principal curvature  (maximum curvature)
  curve[2] = Km - sqrt(Km * Km - Kg); //k_2, Second principal curvature (minimum curvature)
}
 
static apf::Field *getCurves(apf::Field *hessians, apf::Field *gradphi)
{
  apf::Mesh *m = apf::getMesh(hessians);
  apf::Field *curves;
  curves = apf::createLagrangeField(m, "proteus_curves", apf::VECTOR, 1);
  apf::MeshIterator *it = m->begin(0);
  apf::MeshEntity *v;
  while ((v = m->iterate(it)))
  {
    apf::Matrix3x3 hessian;
    apf::getMatrix(hessians, v, 0, hessian);
    apf::Vector3 gphi;
    apf::getVector(gradphi, v, 0, gphi);
    apf::Vector3 curve;
    curveFormula(hessian, gphi, curve);
    apf::setVector(curves, v, 0, curve);
  }
  m->end(it);
  return curves;
}
 
static void clamp(double &v, double min, double max)
//Function used to restrict the mesh size to user specified minimums and maximums
{
  v = std::min(max, v);
  v = std::max(min, v);
}
 
static void clampField(apf::Field *field, double min, double max)
//Function that loops through a field and clamps the values
{
  apf::Mesh *m = apf::getMesh(field);
  int numcomps = apf::countComponents(field);
  double components[numcomps];
  apf::MeshEntity *v;
  apf::MeshIterator *it = m->begin(0);
  while ((v = m->iterate(it)))
  {
    //double tempValue = apf::getScalar(field,v,0);
    apf::getComponents(field, v, 0, &components[0]);
    for (int i = 0; i < numcomps; i++)
      clamp(components[i], min, max);
    //apf::setScalar(field,v,0,tempValue);
    apf::setComponents(field, v, 0, &components[0]);
  }
  m->end(it);
}
 
static void scaleFormula(double phi, double hmin, double hmax,
                         int adapt_step,
                         apf::Vector3 const &curves,
                         apf::Vector3 &scale)
//Function used to set the size scale vector for the interface size field configuration
{
  double epsilon = 7.0 * hmin;
  if (fabs(phi) < epsilon)
  {
    scale[0] = hmin;
    scale[1] = sqrt(0.002 / fabs(curves[1]));
    scale[2] = sqrt(0.002 / fabs(curves[2]));
  }
  else if (fabs(phi) < 4 * epsilon)
  {
    scale[0] = 2 * hmin;
    scale[1] = 2 * sqrt(0.002 / fabs(curves[1]));
    scale[2] = 2 * sqrt(0.002 / fabs(curves[2]));
  }
  else
  {
    scale = apf::Vector3(1, 1, 1) * hmax;
  }
 
  //for (int i = 0; i < 3; ++i)
  //  clamp(scale[i], hmin, hmax);
}
 
static void scaleFormulaERM(double phi, double hmin, double hmax, double h_dest,
                            apf::Vector3 const &curves,
                            double lambda[3], double eps_u, apf::Vector3 &scale,int nsd,double maxAspect)
//Function used to set the size scale vector for the anisotropic ERM size field configuration
//Inputs:
// phi is is the distance to the interface
// hmin is the minimum mesh size
// hmax is the maximum mesh size
// h_dest is the computed new mesh size
// curves is a vector that denotes the curvature of the interface
// lambda is the ordered set of eigenvalues from an eigendecomposition of the metric tensor
// eps_u is a tolerance for the distance away from the interface
//Output:
// scale is the mesh size in each direction for a vertex
{
/*
  double epsilon = 7.0 * hmin;
  double lambdamin = 1.0 / (hmin * hmin);
  if (lambda[1] < 1e-10)
  {
    lambda[1] = lambdamin;
    lambda[2] = lambdamin;
  }
  if (lambda[2] < 1e-10)
  {
    lambda[2] = lambdamin;
  }
*/
  
/*
  scale[0] = h_dest*pow(lambda[1]/lambda[0],0.25);  
  scale[1] = sqrt(lambda[0]/lambda[1])*scale[0];
  scale[2] = 1.0;
*/
/*
  scale[0] = h_dest*pow(lambda[1]/lambda[0],0.25)*pow(3,0.25)*0.5;  
  scale[1] = sqrt(lambda[0]/lambda[1])*scale[0];
  scale[2] = 1.0;
*/
/*
  if(nsd == 2){
    scale[0] = h_dest;  
    scale[1] = sqrt(lambda[0]/lambda[1])*scale[0];
    scale[2] = 1.0;
  }
  else{
    scale[0] = h_dest;  
    scale[1] = sqrt(lambda[0]/lambda[1])*scale[0];
    scale[2] = sqrt(lambda[0]/lambda[2])*scale[0];
  }
*/
 
//3D
 
  if(nsd == 2){
    scale[0] = h_dest * pow((lambda[1] ) / (lambda[0]), 1.0 / 4.0);
    scale[1] = sqrt(lambda[0] / lambda[1]) * scale[0];
    scale[2] = 1.0;
  }
  else{
/*
    scale[0] = h_dest * pow((lambda[1] * lambda[2]) / (lambda[0] * lambda[0]), 1.0 / 6.0);
    scale[1] = sqrt(lambda[0] / lambda[1]) * scale[0];
    scale[2] = sqrt(lambda[0] / lambda[2]) * scale[0];
*/
    scale[0] = h_dest;
    scale[1] = sqrt(lambda[0] / lambda[1]) * scale[0];
    scale[2] = sqrt(lambda[0] / lambda[2]) * scale[0];
    if(scale[1]/scale[0] > maxAspect)
      scale[1] = maxAspect*scale[0];
    if(scale[2]/scale[0] > maxAspect)
      scale[2] = maxAspect*scale[0];
    if(scale[1]/scale[0] > maxAspect || scale[2]/scale[0] > maxAspect){
      logEvent("Scales reached maximum aspect ratio",4);
    }
  }
  //*/
  /*
    if(fabs(phi)<epsilon){
      scale[0] = h_dest;
      scale[1] = sqrt(eps_u/fabs(curves[2]));
      scale[2] = sqrt(eps_u/fabs(curves[1]));
    }
    else
      scale = apf::Vector3(1,1,1) * h_dest; 
*/
}
 
static apf::Field *getSizeScales(apf::Field *phif, apf::Field *curves,
                                 double hmin, double hmax, int adapt_step)
//Function that gets size scales in the interface size field configuration
{
  apf::Mesh *m = apf::getMesh(phif);
  apf::Field *scales;
  scales = apf::createLagrangeField(m, "proteus_size_scale", apf::VECTOR, 1);
  apf::MeshIterator *it = m->begin(0);
  apf::MeshEntity *v;
  while ((v = m->iterate(it)))
  {
    double phi = apf::getScalar(phif, v, 0);
    apf::Vector3 curve;
    apf::getVector(curves, v, 0, curve);
    apf::Vector3 scale;
    scaleFormula(phi, hmin, hmax, adapt_step, curve, scale);
    apf::setVector(scales, v, 0, scale);
  }
  m->end(it);
  return scales;
}
 
struct SortingStruct
{
  apf::Vector3 v;
  double wm;
  bool operator<(const SortingStruct &other) const
  {
    return wm < other.wm;
  }
};
 
static apf::Field *getSizeFrames(apf::Field *hessians, apf::Field *gradphi)
//Function that gets size frames, or the metric tensor, in the interface size field configuration
{
  apf::Mesh *m = apf::getMesh(gradphi);
  apf::Field *frames;
  frames = apf::createLagrangeField(m, "proteus_size_frame", apf::MATRIX, 1);
  apf::MeshIterator *it = m->begin(0);
  apf::MeshEntity *v;
  while ((v = m->iterate(it)))
  {
    apf::Vector3 gphi;
    apf::getVector(gradphi, v, 0, gphi);
    apf::Vector3 dir;
    if (gphi.getLength() > 1e-16)
      dir = gphi.normalize();
    else
      dir = apf::Vector3(1, 0, 0);
    apf::Matrix3x3 hessian;
    apf::getMatrix(hessians, v, 0, hessian);
    apf::Vector3 eigenVectors[3];
    double eigenValues[3];
    apf::eigen(hessian, eigenVectors, eigenValues);
    SortingStruct ssa[3];
    for (int i = 0; i < 3; ++i)
    {
      ssa[i].v = eigenVectors[i];
      ssa[i].wm = std::fabs(eigenValues[i]);
    }
    std::sort(ssa, ssa + 3);
    assert(ssa[2].wm >= ssa[1].wm);
    assert(ssa[1].wm >= ssa[0].wm);
    double firstEigenvalue = ssa[2].wm;
    apf::Matrix3x3 frame;
    frame[0] = dir;
    if (firstEigenvalue > 1e-16)
    {
      apf::Vector3 firstEigenvector = ssa[2].v;
      frame[1] = apf::reject(firstEigenvector, dir);
      frame[2] = apf::cross(frame[0], frame[1]);
      if (frame[2].getLength() < 1e-16)
        frame = apf::getFrame(dir);
    }
    else
      frame = apf::getFrame(dir);
    for (int i = 0; i < 3; ++i)
      frame[i] = frame[i].normalize();
    frame = apf::transpose(frame);
    apf::setMatrix(frames, v, 0, frame);
  }
  m->end(it);
  return frames;
}
 
static apf::Field *getERMSizeFrames(apf::Field *hessians, apf::Field *gradphi, apf::Field *frame_comps[3])
//Function that gets size frames, or the metric tensor, in the interface size field configuration
{
  apf::Mesh *m = apf::getMesh(gradphi);
  apf::Field *frames;
  frames = m->findField("proteus_size_frame");
  apf::MeshIterator *it = m->begin(0);
  apf::MeshEntity *v;
  while ((v = m->iterate(it)))
  {
    apf::Matrix3x3 frame(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0);
    apf::Vector3 gphi;
    apf::getVector(gradphi, v, 0, gphi);
/*
    apf::Vector3 dir;
    if (gphi.getLength() > 1e-16)
      dir = gphi.normalize();
    else
      dir = apf::Vector3(1, 0, 0);
*/
 
    //get eigen values and eigenvectors from hessian
    apf::Matrix3x3 hessian;
    apf::getMatrix(hessians, v, 0, hessian);
    apf::Vector3 eigenVectors[3];
    double eigenValues[3];
    apf::eigen(hessian, eigenVectors, eigenValues);
    SortingStruct ssa[3];
    for (int i = 0; i < 3; ++i)
    {
      ssa[i].v = eigenVectors[i];
      ssa[i].wm = std::fabs(eigenValues[i]);
    }
 
    //sort eigenvalues and eigenvectors
    std::sort(ssa, ssa + 3);
    assert(ssa[2].wm >= ssa[1].wm);
    assert(ssa[1].wm >= ssa[0].wm);
    double firstEigenvalue = ssa[2].wm;
    assert(firstEigenvalue > 1e-12);
    //frame[0] = dir;
    //if (firstEigenvalue > 1e-16)
    //{
      frame[0] = ssa[2].v;
      frame[1] = ssa[1].v;
      frame[2] = ssa[0].v;
    //}
    //else
    //  frame = apf::getFrame(dir);
    
/*
    apf::Vector3 test(1.0,0.0,0.0);
    frame = apf::getFrame(test);
    apf::setMatrix(frames,v,0,frame);
*/
 
    //normalize eigenvectors
    for (int i = 0; i < 3; ++i)
      frame[i] = frame[i].normalize();
    frame = apf::transpose(frame);
    apf::setMatrix(frames, v, 0, frame);
    apf::setVector(frame_comps[0], v, 0, frame[0]);
    apf::setVector(frame_comps[1], v, 0, frame[1]);
    apf::setVector(frame_comps[2], v, 0, frame[2]);
  }
  m->end(it);
  return frames;
}
 
int MeshAdaptPUMIDrvr::calculateAnisoSizeField()
//High level function that obtains the size scales and the size frames for anistropic interface-based adapt
{
  apf::Field *phif = m->findField("phi");
  assert(phif);
  apf::Field *gradphi = apf::recoverGradientByVolume(phif);
  apf::Field *grad2phi = apf::recoverGradientByVolume(gradphi);
  apf::Field *hess = computeHessianField(grad2phi);
  apf::destroyField(grad2phi);
  apf::Field *curves = getCurves(hess, gradphi);
  freeField(size_scale);
 
  size_scale = getSizeScales(phif, curves, hmin, hmax, nAdapt);
  apf::destroyField(curves);
  freeField(size_frame);
  size_frame = getSizeFrames(hess, gradphi);
  apf::destroyField(hess);
 
  apf::destroyField(gradphi);
  for (int i = 0; i < 2; ++i)
    SmoothField(size_scale);
 
  return 0;
}
 
struct Smoother : public apf::CavityOp
{
  Smoother(apf::Field *f) : apf::CavityOp(apf::getMesh(f))
  {
    field = f;
    int nc = apf::countComponents(f);
    newField = apf::createPackedField(mesh, "proteus_smooth_new", nc);
    sum.setSize(nc);
    value.setSize(nc);
    nApplied = 0;
  }
  ~Smoother()
  {
    copyData(field, newField);
    apf::destroyField(newField);
  }
  virtual Outcome setEntity(apf::MeshEntity *e)
  {
    if (apf::hasEntity(newField, e))
      return SKIP;
    if (!this->requestLocality(&e, 1))
      return REQUEST;
    vertex = e;
    return OK;
  }
  virtual void apply()
  {
    /* the smoothing function used here is the average of the
   vertex value and neighboring vertex values, with the
   center vertex weighted equally as the neighbors */
    apf::Up edges;
    mesh->getUp(vertex, edges);
    apf::getComponents(field, vertex, 0, &sum[0]);
    for (int i = 0; i < edges.n; ++i)
    {
      apf::MeshEntity *ov = apf::getEdgeVertOppositeVert(
          mesh, edges.e[i], vertex);
      apf::getComponents(field, ov, 0, &value[0]);
      sum += value;
    }
    sum /= edges.n + 1;
    apf::setComponents(newField, vertex, 0, &sum[0]);
    ++nApplied;
  }
  apf::Field *field;
  apf::Field *newField;
  apf::MeshEntity *vertex;
  apf::DynamicVector sum;
  apf::DynamicVector value;
  apf::MeshTag *emptyTag;
  int nApplied;
};
 
static void SmoothField(apf::Field *f)
{
  Smoother op(f);
  op.applyToDimension(0);
}
 
void getTargetError(apf::Mesh* m, apf::Field* errField, double &target_error,double totalError, PCU_t PCUObj){
  //Implemented for 3D and for serial case only so far
  //Need to communicate target error in parallel
  logEvent("WARNING/ERROR:Parallel implementation is not completed yet",4);
  if(m->getDimension()==2){
    target_error = totalError/sqrt(m->count(m->getDimension()));
    if(PCU_Comm_Self(PCUObj)==0)
      std::cout<<"The estimated target error is "<<target_error<<std::endl;
    return;
  }
  apf::Field* interfaceField = m->findField("vof");
  apf::Field* targetField = apf::createField(m,"targetError",apf::SCALAR,apf::getVoronoiShape(m->getDimension(),1));
  apf::MeshEntity* ent;
  apf::MeshIterator* it = m->begin(m->getDimension());
  apf::MeshElement* element;
  apf::Element* vofElem;
  std::vector <double> errVect;
  while( (ent = m->iterate(it))){
    element = apf::createMeshElement(m, ent);
    vofElem = apf::createElement(interfaceField,element);
    double vofVal = apf::getScalar(vofElem,apf::Vector3(1./3.,1./3.,1./3.));
 
    if(vofVal < 0.9 && vofVal > 0.1){ //at the interface
      double errorValue = apf::getScalar(errField,ent,0);
      errVect.push_back(errorValue);    
      apf::setScalar(targetField,ent,0,errorValue);
    }
    else{
      apf::setScalar(targetField,ent,0,0.0);
    }
  }
  m->end(it);
  if(errVect.size()==0){
    target_error = totalError/sqrt(m->count(m->getDimension()));
  }
  else{
    std::ofstream myfile;
    myfile.open("interfaceErrors.txt", std::ios::app );
    for(int i=0;i<errVect.size();i++){
      myfile << errVect[i]<<std::endl;
    }
    myfile.close();
    std::sort(errVect.begin(),errVect.end());
    int vectorSize = errVect.size();
    if(vectorSize %2 ==0){
      int idx1 = vectorSize/2-1;
      target_error = (errVect[idx1]+errVect[idx1+1])/2; //get average
    }
    else
      target_error = errVect[(vectorSize-1)/2];
  }
  if(PCU_Comm_Self(PCUObj)==0)
    std::cout<<"The estimated target error is "<<target_error<<std::endl;
  //std::abort();
}
 
int MeshAdaptPUMIDrvr::getERMSizeField(double err_total)
//High level function that obtains the size scales and the size frames for ERM-based adapt and uses the computed total error
{
 
  freeField(size_frame);
  freeField(size_scale);
  freeField(size_iso);
 
  //Initialize fields and needed types/variables
  apf::Field* errField;
  //apf::Mesh* m;
  if(hasERM)
    errField = m->findField("ErrorRegion");
  else if(hasVMS)
    errField = m->findField("VMSH1");
  assert(errField); 
  //apf::Mesh *m = apf::getMesh(vmsErrH1);
  //apf::getMesh(errField);
  apf::MeshIterator *it;
  apf::MeshEntity *v;
  apf::MeshElement *element;
  apf::MeshEntity *reg;
  //size_iso = apf::createLagrangeField(m, "proteus_size", apf::SCALAR, 1);
  if(m->findField("errorSize"))
    apf::destroyField(m->findField("errorSize"));
  apf::Field *errorSize = apf::createLagrangeField(m, "errorSize", apf::SCALAR, 1);
 
  if (adapt_type_config == "anisotropic"){
    size_scale = apf::createLagrangeField(m, "proteus_size_scale", apf::VECTOR, 1);
    size_frame = apf::createLagrangeField(m, "proteus_size_frame", apf::MATRIX, 1);
  }
  apf::Field *errorSize_reg = apf::createField(m, "iso_size", apf::SCALAR, apf::getConstant(nsd));
  apf::Field *clipped_vtx = apf::createLagrangeField(m, "iso_clipped", apf::SCALAR, 1);
  
  //Get total number of elements
  int numel = 0;
  int nsd = m->getDimension();
  numel = m->count(nsd);
  PCU_Add_Ints(PCUObj, &numel, 1);
 
  //if target error is not specified, choose one based on interface or based on equidistribution assumption
  if(target_error==0){
    if(m->findField("vof")!=NULL)
      getTargetError(m,errField,target_error,err_total,PCUObj);
    else
      target_error = err_total/sqrt(m->count(nsd));
  }
   
  // Get domain volume
  // should only need to be computed once unless geometry is complex
  if (domainVolume == 0)
  {
    double volTotal = 0.0;
    it = m->begin(nsd);
    while (reg = m->iterate(it))
    {
      volTotal += apf::measure(m, reg);
    }
    PCU_Add_Doubles(PCUObj, &volTotal, 1);
    domainVolume = volTotal;
    assert(domainVolume > 0);
  }
  //compute the new size field over elements
  it = m->begin(nsd);
  double err_curr = 0.0;
  double errRho_curr = 0.0;
  double errRho_target = target_error / sqrt(domainVolume);
  apf::Vector3 err_vect;
  while (reg = m->iterate(it))
  {
    double h_old;
    double h_new;
    element = apf::createMeshElement(m, reg);
 
    if (m->getDimension() == 2)
      //h_old = sqrt(apf::measure(element) * 4 / sqrt(3));
      h_old = apf::computeLargestHeightInTri(m,reg);
    else
      //h_old = pow(apf::measure(element) * 6 * sqrt(2), 1.0 / 3.0); //edge of a regular tet
      h_old = apf::computeLargestHeightInTet(m,reg);
    //err_curr = apf::getScalar(vmsErrH1, reg, 0);
    err_curr = apf::getScalar(errField, reg, 0);
    //err_curr = err_vect[0];
    //errRho_curr = apf::getScalar(errRho_reg, reg, 0);
    //h_new = h_old*errRho_target/errRho_curr;
    //h_new = h_old*sqrt(apf::measure(element))/sqrt(domainVolume)*target_error/err_curr;
    //
    //error-to-size relationship should be different between anisotropic and isotropic cases
    //consider moving this to where size frames are computed to get aspect ratio info
    if (adapt_type_config == "anisotropic")
      if(target_error/err_curr <= 1)
        h_new = h_old * pow((target_error / err_curr),2.0/(2.0*(1.0)+1.0)); //refinement
      else
        h_new = h_old * pow((target_error / err_curr),2.0/(2.0*(1.0)+3.0)); //coarsening
    else //isotropic
    {
      if(target_error/err_curr <= 1)
        h_new = h_old * pow((target_error / err_curr),2.0/(2.0*(1.0)+nsd));
      else
        h_new = h_old * pow((target_error / err_curr),2.0/(2.0*(1.0)+nsd));
        //h_new = h_old * pow((target_error / err_curr),2.0/(2.0*(1.0)+nsd+1.0)); //extra 1.0 to slow down coarsening
    }
 
    apf::setScalar(errorSize_reg, reg, 0, h_new);
    apf::destroyMeshElement(element);
  }
  m->end(it);
 
  //Transfer size field from elements to vertices through averaging
  PCU_Comm_Begin(PCUObj);
  it = m->begin(0);
  while ((v = m->iterate(it)))
  {
    //averageToEntity(errorSize_reg, errorSize, v);
    //volumeAverageToEntity(errorSize_reg, errorSize, v);
    //errorAverageToEntity(errorSize_reg, errorSize,errField, v);
    minToEntity(errorSize_reg, errorSize, v);
    if(!m->isOwned(v))
    {
        apf::Copies remotes;
        m->getRemotes(v,remotes);
        int owningPart=m->getOwner(v);
        PCU_COMM_PACK(PCUObj, owningPart, remotes[owningPart]);
        double currentSize = apf::getScalar(errorSize,v,0);
        PCU_COMM_PACK(PCUObj, owningPart,currentSize);
    }
 
  }
  m->end(it);
  PCU_Comm_Send(PCUObj);
  while(PCU_Comm_Receive(PCUObj))
  {
    apf::MeshEntity* receivedEnt;
    double receivedSize;
    PCU_COMM_UNPACK(PCUObj, receivedEnt);
    PCU_COMM_UNPACK(PCUObj, receivedSize);
 
    double currentSize = apf::getScalar(errorSize,receivedEnt,0);
 
    if(receivedSize<currentSize)
    {
        apf::setScalar(errorSize,receivedEnt,0,receivedSize);
    }
 
  }
 
 
 
  //Get the anisotropic size frame
  if (adapt_type_config == "anisotropic")
  {
    if(comm_rank==0)
      std::cout<<"Entering anisotropic loop to compute size scales and frames\n";
    double eps_u = 0.002; //distance from the interface
/*
    apf::Field *phif = m->findField("phi");
    apf::Field *gradphi = apf::recoverGradientByVolume(phif);
    apf::Field *grad2phi = apf::recoverGradientByVolume(gradphi);
*/
    apf::Field *speedF = extractSpeed(m->findField("velocity"));
    apf::Field *gradSpeed = apf::recoverGradientByVolume(speedF);
    apf::Field *grad2Speed = apf::recoverGradientByVolume(gradSpeed);
    //apf::Field *hess = computeHessianField(grad2phi);
    //apf::Field *curves = getCurves(hess, gradphi);
    //apf::Field* metricf = computeMetricField(gradphi,grad2phi,errorSize,eps_u);
    apf::Field *metricf = computeMetricField(gradSpeed, grad2Speed, errorSize, eps_u);
    apf::Field *frame_comps[3] = {apf::createLagrangeField(m, "frame_0", apf::VECTOR, 1), apf::createLagrangeField(m, "frame_1", apf::VECTOR, 1), apf::createLagrangeField(m, "frame_2", apf::VECTOR, 1)};
    //getERMSizeFrames(metricf, gradSpeed, frame_comps);
 
    //Set the size scale for vertices
    it = m->begin(0);
    apf::Vector3 scale;
    while ((v = m->iterate(it)))
    {
      double tempScale = apf::getScalar(errorSize, v, 0);
      if (tempScale < hmin)
        apf::setScalar(clipped_vtx, v, 0, -1);
      else if (tempScale > hmax)
        apf::setScalar(clipped_vtx, v, 0, 1);
      else
        apf::setScalar(clipped_vtx, v, 0, 0);
      clamp(tempScale, hmin, hmax);
      apf::setScalar(errorSize,v,0,tempScale);
    }
    it = m->begin(0);
    while( (v = m->iterate(it)) ){
      double phi;// = apf::getScalar(phif, v, 0);
      apf::Vector3 curve;
      //apf::getVector(curves, v, 0, curve);
 
      //metricf is the hessian
      apf::Matrix3x3 metric;
      apf::getMatrix(metricf, v, 0, metric);
 
      apf::Vector3 eigenVectors[3];
      double eigenValues[3];
      apf::eigen(metric, eigenVectors, eigenValues);
      // Sort the eigenvalues and corresponding vectors
      // Larger eigenvalues means a need for a finer mesh
      SortingStruct ssa[3];
      for (int i = 0; i < 3; ++i)
      {
        ssa[i].v = eigenVectors[i];
        ssa[i].wm = std::fabs(eigenValues[i]);
      }
      std::sort(ssa, ssa + 3);
 
      assert(ssa[2].wm >= ssa[1].wm);
      assert(ssa[1].wm >= ssa[0].wm);
 
      double lambda[3] = {ssa[2].wm, ssa[1].wm, ssa[0].wm};
 
      scaleFormulaERM(phi, hmin, hmax, apf::getScalar(errorSize, v, 0), curve, lambda, eps_u, scale,nsd,maxAspect);
      apf::setVector(size_scale, v, 0, scale);
      //get frames
 
      apf::Matrix3x3 frame(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0);
 
      //get eigen values and eigenvectors from hessian
      double firstEigenvalue = ssa[2].wm;
      assert(firstEigenvalue > 1e-12);
      frame[0] = ssa[2].v;
      frame[1] = ssa[1].v;
      frame[2] = ssa[0].v;
    
      //normalize eigenvectors
      for (int i = 0; i < 3; ++i)
        frame[i] = frame[i].normalize();
      frame = apf::transpose(frame);
      apf::setMatrix(size_frame, v, 0, frame);
 
    }
    m->end(it);
 
    //Do simple size and aspect ratio grading
    gradeAnisoMesh(m,gradingFactor,PCUObj);
    if(comm_rank==0)
      std::cout<<"Finished grading size 0\n";
    gradeAspectRatio(m,1,gradingFactor,PCUObj);
    if(comm_rank==0)
      std::cout<<"Finished grading size 1\n";
    gradeAspectRatio(m,2,gradingFactor,PCUObj);
    if(comm_rank==0)
      std::cout<<"Finished grading size 2\n";
 
    apf::synchronize(size_scale);
 
    //apf::destroyField(gradphi);
    //apf::destroyField(grad2phi);
    //apf::destroyField(curves);
    //apf::destroyField(hess);
 
    if(logging_config=="on"){
      char namebuffer[20];
      sprintf(namebuffer,"pumi_preadapt_aniso_%i",nAdapt);
      apf::writeVtkFiles(namebuffer, m);
    }
 
    apf::destroyField(metricf);
    apf::destroyField(frame_comps[0]);
    apf::destroyField(frame_comps[1]);
    apf::destroyField(frame_comps[2]);
    apf::destroyField(speedF);
    apf::destroyField(gradSpeed);
    apf::destroyField(grad2Speed);
  }
  else
  {
    it = m->begin(0);
    while ((v = m->iterate(it)))
    {
      double tempScale = apf::getScalar(errorSize, v, 0);
      if (tempScale < hmin)
        apf::setScalar(clipped_vtx, v, 0, -1);
      else if (tempScale > hmax)
        apf::setScalar(clipped_vtx, v, 0, 1);
      else
        apf::setScalar(clipped_vtx, v, 0, 0);
      clamp(tempScale, hmin, hmax);
      apf::setScalar(errorSize, v, 0, tempScale);
    }
    //gradeMesh();
    apf::synchronize(errorSize);
    m->end(it);
    if (target_element_count != 0)
    {
      sam::scaleIsoSizeField(errorSize, target_element_count);
      clampField(errorSize, hmin, hmax);
      //gradeMesh();
      //SmoothField(errorSize);
    }
    sizeFieldList.push(errorSize);
  } 
 
  //Destroy locally required fields
  apf::destroyField(errorSize_reg);
  apf::destroyField(clipped_vtx);
  if(comm_rank==0)
    std::cout<<"Finished Size Field\n";
  return 0;
}
 
int MeshAdaptPUMIDrvr::testIsotropicSizeField()
//Function that tests MeshAdapt by generating an isotropic sizefield based on hmin
{ 
    freeField(size_iso);
    size_iso = apf::createLagrangeField(m, "proteus_size",apf::SCALAR,1);
    apf::MeshIterator* it = m->begin(0);
    apf::MeshEntity* v;
    while(v = m->iterate(it)){
      double phi = hmin;
      clamp(phi,hmin,hmax);
      apf::setScalar(size_iso,v,0,phi);
    }
    return 0;
}
 
int gradeSizeModify(apf::Mesh* m, double gradingFactor, 
    double size[2], apf::Adjacent edgAdjVert, 
    apf::Adjacent vertAdjEdg,
    std::queue<apf::MeshEntity*> &markedEdges,
    apf::MeshTag* isMarked,
    int fieldType,
    int vecPos, //which idx of sizeVec to modify
    int idxFlag,
    PCU_t PCUObj)
 
//General function to actually modify sizes
{
    //Determine a switching scheme depending on which vertex needs a modification
    int idx1,idx2;
    if(idxFlag == 0){
      idx1=0;
      idx2=1;
    } 
    else{
      idx1=1;
      idx2 = 0;
    } 
    
    int marker[3] = {0,1,0}; 
    double marginVal = 0.01;
    int needsParallel=0;
 
    if(fieldType == apf::SCALAR){
      apf::Field* size_iso = m->findField("proteus_size");
 
      if(size[idx1]>(gradingFactor*size[idx2])*(1+marginVal))
      {
        if(m->isOwned(edgAdjVert[idx1]))
        {
          size[idx1] = gradingFactor*size[idx2];
          apf::setScalar(size_iso,edgAdjVert[idx1],0,size[idx1]);
          m->getAdjacent(edgAdjVert[idx1], 1, vertAdjEdg);
          for (std::size_t i=0; i<vertAdjEdg.getSize();++i){
            m->getIntTag(vertAdjEdg[i],isMarked,&marker[2]);
            //if edge is not already marked
            if(!marker[2]){
              m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]);
              markedEdges.push(vertAdjEdg[i]);
            }
          }
        } //end isOwned
        else
        { //Pack information to owning processor
          needsParallel=1;
          apf::Copies remotes;
          m->getRemotes(edgAdjVert[idx1],remotes);
          double newSize = gradingFactor*size[idx2];
          int owningPart=m->getOwner(edgAdjVert[idx1]);
          PCU_COMM_PACK(PCUObj, owningPart, remotes[owningPart]);
          PCU_COMM_PACK(PCUObj, owningPart,newSize);
        }
      }
 
    }//end if apf::SCALAR
    else{
      apf::Field* size_scale = m->findField("proteus_size_scale");
      apf::Vector3 sizeVec;
      if(size[idx1]>(gradingFactor*size[idx2])*(1+marginVal)){
        size[idx1] = gradingFactor*size[idx2];
        apf::getVector(size_scale,edgAdjVert[idx1],0,sizeVec);
        if(vecPos > 0){
          sizeVec[vecPos] = size[idx1]*sizeVec[0]; //realize the new aspect ratio
        }
        else{
          sizeVec[0] = size[idx1];
        }
        apf::setVector(size_scale,edgAdjVert[idx1],0,sizeVec);
        m->getAdjacent(edgAdjVert[idx1], 1, vertAdjEdg);
        for (std::size_t i=0; i<vertAdjEdg.getSize();++i){
          m->getIntTag(vertAdjEdg[i],isMarked,&marker[2]);
          //if edge is not already marked
          if(!marker[2]){
            m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]);
            markedEdges.push(vertAdjEdg[i]);
          }
        }
      }
    }
  return needsParallel;
}
 
void markEdgesInitial(apf::Mesh* m, std::queue<apf::MeshEntity*> &markedEdges,double gradingFactor)
//Function used to initially determine which edges need to be considered for gradation
{
  //marker structure for 0) not marked 1) marked 2)storage
  int marker[3] = {0,1,0}; 
 
  double size[2];
  apf::MeshTag* isMarked = m->findTag("isMarked");
  apf::Field* size_iso = m->findField("proteus_size");
  apf::Adjacent edgAdjVert;
  apf::MeshEntity* edge;
  apf::MeshIterator* it = m->begin(1);
  while((edge=m->iterate(it))){
    m->getAdjacent(edge, 0, edgAdjVert);
    for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){
      size[i]=apf::getScalar(size_iso,edgAdjVert[i],0);
    }
    if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){
      //add edge to a queue 
      markedEdges.push(edge);
      //tag edge to indicate that it is part of queue 
      m->setIntTag(edge,isMarked,&marker[1]); 
    }
    else{
      m->setIntTag(edge,isMarked,&marker[0]); 
    }
  }
  m->end(it); 
}
 
int serialGradation(apf::Mesh* m, std::queue<apf::MeshEntity*> &markedEdges,double gradingFactor, PCU_t PCUObj)
//Function used loop over the mesh edge queue for gradation and modify the sizes
{
  double size[2];
  //marker structure for 0) not marked 1) marked 2)storage
  int marker[3] = {0,1,0}; 
  apf::MeshTag* isMarked = m->findTag("isMarked");
  apf::Field* size_iso = m->findField("proteus_size");
  apf::Adjacent edgAdjVert;
  apf::Adjacent vertAdjEdg;
  apf::MeshEntity* edge;
  apf::MeshIterator* it = m->begin(1);
  int needsParallel=0;
 
  //perform serial gradation while packing necessary info for parallel
  while(!markedEdges.empty()){ 
    edge = markedEdges.front();
    m->getAdjacent(edge, 0, edgAdjVert);
    for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){
      size[i] = apf::getScalar(size_iso,edgAdjVert[i],0);
    }
 
    needsParallel+=gradeSizeModify(m, gradingFactor, size, edgAdjVert, 
      vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 0, PCUObj);
    needsParallel+=gradeSizeModify(m, gradingFactor, size, edgAdjVert, 
      vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 1, PCUObj);
 
    m->setIntTag(edge,isMarked,&marker[0]);
    markedEdges.pop();
  }
  return needsParallel;
}
 
int MeshAdaptPUMIDrvr::gradeMesh(double gradationFactor)
//Function to grade isotropic mesh through comparison of edge vertex size ratios
//This implementation accounts for parallel meshes as well
//First do serial gradation. 
//If a shared entity has its size modified, then send new size to owning copy.
//After full loop over entities, have owning copy take minimum of all sizes received
//Flag adjacent entities to owning copy.
//Communicate to remote copies that a size was modified, and flag adjacent edges to remote copies for further gradation
{
  //
  logEvent("Starting grading",4);
  apf::MeshEntity* edge;
  apf::Adjacent edgAdjVert;
  apf::Adjacent vertAdjEdg;
  double size[2];
  std::queue<apf::MeshEntity*> markedEdges;
  apf::MeshTag* isMarked = m->createIntTag("isMarked",1);
 
  //marker structure for 0) not marked 1) marked 2)storage
  int marker[3] = {0,1,0}; 
 
  apf::MeshIterator* it;
  markEdgesInitial(m,markedEdges,gradingFactor);
 
  int needsParallel=1;
  int nCount=1;
  while(needsParallel)
  {
    PCU_Comm_Begin(PCUObj);
    needsParallel = serialGradation(m,markedEdges,gradingFactor,PCUObj);
 
    PCU_Add_Ints(PCUObj, &needsParallel,1);
    if(comm_rank==0)
      std::cerr<<"Sending size info for gradation"<<std::endl;
    PCU_Comm_Send(PCUObj); 
 
    apf::MeshEntity* ent;
    double receivedSize;
    double currentSize;
    double newSize;
 
    //Need a container to get all entitites that need to be updated on remotes
    std::queue<apf::MeshEntity*> updateRemoteVertices;
 
    apf::Copies remotes;
    //owning copies are receiving
    while(PCU_Comm_Receive(PCUObj))
    {
      PCU_COMM_UNPACK(PCUObj, ent);
      PCU_COMM_UNPACK(PCUObj, receivedSize);
 
      if(!m->isOwned(ent)){
        logEvent("THERE WAS AN ERROR",4);
        std::exit(1);
      }
 
      currentSize = apf::getScalar(size_iso,ent,0);
      newSize = std::min(receivedSize,currentSize);
      apf::setScalar(size_iso,ent,0,newSize);
      
      //add adjacent edges into Q
      m->getAdjacent(ent, 1, vertAdjEdg);
      for (std::size_t i=0; i<vertAdjEdg.getSize();++i)
      {
        edge = vertAdjEdg[i];
        m->getIntTag(vertAdjEdg[i],isMarked,&marker[2]);
        if(!marker[2])
        {
          markedEdges.push(edge);
          //tag edge to indicate that it is part of queue 
          m->setIntTag(edge,isMarked,&marker[1]);
        }
      }
      updateRemoteVertices.push(ent);
    }
 
    PCU_Comm_Begin(PCUObj);
 
    while(!updateRemoteVertices.empty())
    { 
      ent = updateRemoteVertices.front();
      //get remote copies and send updated mesh sizes
      m->getRemotes(ent,remotes);
      currentSize = apf::getScalar(size_iso,ent,0);
      for(apf::Copies::iterator iter=remotes.begin(); iter!=remotes.end();++iter)
      {
        PCU_COMM_PACK(PCUObj, iter->first, iter->second);
      }
      updateRemoteVertices.pop();
    }
 
    PCU_Comm_Send(PCUObj);
    //while remote copies are receiving
    while(PCU_Comm_Receive(PCUObj))
    {
      //unpack
      PCU_COMM_UNPACK(PCUObj, ent);
      //PCU_COMM_UNPACK(PCUObj, receivedSize);
      assert(!m->isOwned(ent));
 
      if(m->isOwned(ent)){
        logEvent("Problem occurred",4);
        std::exit(1);
      }
 
      //add adjacent edges into Q
      m->getAdjacent(ent, 1, vertAdjEdg);
      for (std::size_t i=0; i<vertAdjEdg.getSize();++i)
      {
        edge = vertAdjEdg[i];
        m->getIntTag(vertAdjEdg[i],isMarked,&marker[2]);
        if(!marker[2])
        {
          markedEdges.push(edge);
          //tag edge to indicate that it is part of queue 
          m->setIntTag(edge,isMarked,&marker[1]);
        }
      }
    }
    apf::synchronize(size_iso);
 
  } //end outer while
 
  //Cleanup of edge marker field
  it = m->begin(1);
  while((edge=m->iterate(it))){
    m->removeTag(edge,isMarked);
  }
  m->end(it); 
  m->destroyTag(isMarked);
 
  //apf::synchronize(size_iso);
  logEvent("Completed grading",4);
  return needsParallel;
}
 
void gradeAnisoMesh(apf::Mesh* m, PCU_t PCUObj)
//Function to grade anisotropic mesh through comparison of edge vertex aspect ratios and minimum sizes
//For simplicity, we do not bother with accounting for entities across partitions
{
  //
  //if(comm_rank==0)
  //  std::cout<<"Starting anisotropic grading\n";
  apf::MeshIterator* it = m->begin(1);
  apf::MeshEntity* edge;
  apf::Adjacent edgAdjVert;
  apf::Adjacent vertAdjEdg;
  double gradingFactor = 1.3;
  double size[2];
  apf::Vector3 sizeVec;
  std::queue<apf::MeshEntity*> markedEdges;
  apf::MeshTag* isMarked = m->createIntTag("isMarked",1);
  apf::Field* size_scale = m->findField("proteus_size_scale");
 
  //marker structure for 0) not marked 1) marked 2)storage
  int marker[3] = {0,1,0}; 
 
  while((edge=m->iterate(it))){
    m->getAdjacent(edge, 0, edgAdjVert);
    for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){
      apf::getVector(size_scale,edgAdjVert[i],0,sizeVec);
      size[i]=sizeVec[0];
    }
    if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){
      //add edge to a queue 
      markedEdges.push(edge);
      //tag edge to indicate that it is part of queue 
      m->setIntTag(edge,isMarked,&marker[1]); 
 
    }
    else{
      m->setIntTag(edge,isMarked,&marker[0]); 
    }
  }
  m->end(it); 
  while(!markedEdges.empty()){
    edge = markedEdges.front();
    m->getAdjacent(edge, 0, edgAdjVert);
    for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){
      apf::getVector(size_scale,edgAdjVert[i],0,sizeVec);
      size[i]=sizeVec[0];
    }
    gradeSizeModify(m, gradingFactor, size, edgAdjVert, 
      vertAdjEdg, markedEdges, isMarked, apf::VECTOR,0, 0, PCUObj);
    gradeSizeModify(m, gradingFactor, size, edgAdjVert, 
      vertAdjEdg, markedEdges, isMarked, apf::VECTOR,0, 1, PCUObj);
 
/*
    if(size[0]>gradingFactor*size[1]){
      size[0] = gradingFactor*size[1];
      apf::getVector(size_scale,edgAdjVert[0],0,sizeVec);
      sizeVec[0] = size[0];
      apf::setVector(size_scale,edgAdjVert[0],0,sizeVec);
      m->getAdjacent(edgAdjVert[0], 1, vertAdjEdg);
      for (std::size_t i=0; i<vertAdjEdg.getSize();++i){
        m->getIntTag(vertAdjEdg[i],isMarked,&marker[2]);
        //if edge is not already marked
        if(!marker[2]){
          m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]);
          markedEdges.push(vertAdjEdg[i]);
        }
      }
    }
    if(size[1]>gradingFactor*size[0]){
      size[1] = gradingFactor*size[0];
      apf::getVector(size_scale,edgAdjVert[1],0,sizeVec);
      sizeVec[0] = size[1];
      apf::setVector(size_scale,edgAdjVert[1],0,sizeVec);
      m->getAdjacent(edgAdjVert[1], 1, vertAdjEdg);
      for (std::size_t i=0; i<vertAdjEdg.getSize();++i){
        m->getIntTag(vertAdjEdg[i],isMarked,&marker[2]);
        //if edge is not already marked
        if(!marker[2]){
          m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]);
          markedEdges.push(vertAdjEdg[i]);
        }
      }
    }
*/
    m->setIntTag(edge,isMarked,&marker[0]);
    markedEdges.pop();
  }
  it = m->begin(1);
  while((edge=m->iterate(it))){
    m->removeTag(edge,isMarked);
  }
  m->end(it); 
  m->destroyTag(isMarked);
  apf::synchronize(size_scale);
  //if(comm_rank==0)
  //  std::cout<<"Completed minimum size grading\n";
}
 
void gradeAspectRatio(apf::Mesh* m,int idx, PCU_t PCUObj)
//Function to grade anisotropic mesh through comparison of edge vertex aspect ratios and minimum sizes
//For simplicity, we do not bother with accounting for entities across partitions
{
  std::cout<<"Entered function\n"; 
  apf::MeshIterator* it = m->begin(1);
  apf::MeshEntity* edge;
  apf::Adjacent edgAdjVert;
  apf::Adjacent vertAdjEdg;
  double gradingFactor = 1.3;
  double size[2];
  apf::Vector3 sizeVec;
  std::queue<apf::MeshEntity*> markedEdges;
  apf::MeshTag* isMarked = m->createIntTag("isMarked",1);
  apf::Field* size_scale = m->findField("proteus_size_scale");
 
  //marker structure for 0) not marked 1) marked 2)storage
  int marker[3] = {0,1,0}; 
 
  while((edge=m->iterate(it))){
    m->getAdjacent(edge, 0, edgAdjVert);
    for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){
      apf::getVector(size_scale,edgAdjVert[i],0,sizeVec);
      size[i]=sizeVec[idx]/sizeVec[0];
    }
    if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){
      //add edge to a queue 
      markedEdges.push(edge);
      //tag edge to indicate that it is part of queue 
      m->setIntTag(edge,isMarked,&marker[1]); 
 
    }
    else{
      m->setIntTag(edge,isMarked,&marker[0]); 
    }
  }
  m->end(it); 
 
  std::cout<<"Got queue of size "<<markedEdges.size()<<std::endl; 
  while(!markedEdges.empty()){
    edge = markedEdges.front();
    m->getAdjacent(edge, 0, edgAdjVert);
    for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){
      apf::getVector(size_scale,edgAdjVert[i],0,sizeVec);
      size[i]=sizeVec[idx]/sizeVec[0];
    }
    gradeSizeModify(m, gradingFactor, size, edgAdjVert, 
      vertAdjEdg, markedEdges, isMarked, apf::VECTOR, idx, 0, PCUObj);
    gradeSizeModify(m, gradingFactor, size, edgAdjVert, 
      vertAdjEdg, markedEdges, isMarked, apf::VECTOR, idx, 1, PCUObj);
 
    m->setIntTag(edge,isMarked,&marker[0]);
    markedEdges.pop();
  }
  it = m->begin(1);
  while((edge=m->iterate(it))){
    m->removeTag(edge,isMarked);
  }
  m->end(it); 
  m->destroyTag(isMarked);
  apf::synchronize(size_scale);
}
 
void gradeAnisoMesh(apf::Mesh* m,double gradingFactor,PCU_t PCUObj)
//Function to grade anisotropic mesh through comparison of edge vertex aspect ratios and minimum sizes
//For simplicity, we do not bother with accounting for entities across partitions
{
  //
  //if(comm_rank==0)
  //  std::cout<<"Starting anisotropic grading\n";
  apf::MeshIterator* it = m->begin(1);
  apf::MeshEntity* edge;
  apf::Adjacent edgAdjVert;
  apf::Adjacent vertAdjEdg;
  //double gradingFactor = 1.3;
  double size[2];
  apf::Vector3 sizeVec;
  std::queue<apf::MeshEntity*> markedEdges;
  apf::MeshTag* isMarked = m->createIntTag("isMarked",1);
  apf::Field* size_scale = m->findField("proteus_size_scale");
 
  //marker structure for 0) not marked 1) marked 2)storage
  int marker[3] = {0,1,0}; 
 
  while((edge=m->iterate(it))){
    m->getAdjacent(edge, 0, edgAdjVert);
    for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){
      apf::getVector(size_scale,edgAdjVert[i],0,sizeVec);
      size[i]=sizeVec[0];
    }
    if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){
      //add edge to a queue 
      markedEdges.push(edge);
      //tag edge to indicate that it is part of queue 
      m->setIntTag(edge,isMarked,&marker[1]); 
 
    }
    else{
      m->setIntTag(edge,isMarked,&marker[0]); 
    }
  }
  m->end(it); 
  while(!markedEdges.empty()){
    edge = markedEdges.front();
    m->getAdjacent(edge, 0, edgAdjVert);
    for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){
      apf::getVector(size_scale,edgAdjVert[i],0,sizeVec);
      size[i]=sizeVec[0];
    }
    gradeSizeModify(m, gradingFactor, size, edgAdjVert, 
      vertAdjEdg, markedEdges, isMarked, apf::VECTOR,0, 0, PCUObj);
    gradeSizeModify(m, gradingFactor, size, edgAdjVert, 
      vertAdjEdg, markedEdges, isMarked, apf::VECTOR,0, 1, PCUObj);
 
/*
    if(size[0]>gradingFactor*size[1]){
      size[0] = gradingFactor*size[1];
      apf::getVector(size_scale,edgAdjVert[0],0,sizeVec);
      sizeVec[0] = size[0];
      apf::setVector(size_scale,edgAdjVert[0],0,sizeVec);
      m->getAdjacent(edgAdjVert[0], 1, vertAdjEdg);
      for (std::size_t i=0; i<vertAdjEdg.getSize();++i){
        m->getIntTag(vertAdjEdg[i],isMarked,&marker[2]);
        //if edge is not already marked
        if(!marker[2]){
          m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]);
          markedEdges.push(vertAdjEdg[i]);
        }
      }
    }
    if(size[1]>gradingFactor*size[0]){
      size[1] = gradingFactor*size[0];
      apf::getVector(size_scale,edgAdjVert[1],0,sizeVec);
      sizeVec[0] = size[1];
      apf::setVector(size_scale,edgAdjVert[1],0,sizeVec);
      m->getAdjacent(edgAdjVert[1], 1, vertAdjEdg);
      for (std::size_t i=0; i<vertAdjEdg.getSize();++i){
        m->getIntTag(vertAdjEdg[i],isMarked,&marker[2]);
        //if edge is not already marked
        if(!marker[2]){
          m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]);
          markedEdges.push(vertAdjEdg[i]);
        }
      }
    }
*/
    m->setIntTag(edge,isMarked,&marker[0]);
    markedEdges.pop();
  }
  it = m->begin(1);
  while((edge=m->iterate(it))){
    m->removeTag(edge,isMarked);
  }
  m->end(it); 
  m->destroyTag(isMarked);
  apf::synchronize(size_scale);
  //if(comm_rank==0)
  //  std::cout<<"Completed minimum size grading\n";
}
 
void gradeAspectRatio(apf::Mesh* m,int idx,double gradingFactor,PCU_t PCUObj)
//Function to grade anisotropic mesh through comparison of edge vertex aspect ratios and minimum sizes
//For simplicity, we do not bother with accounting for entities across partitions
{
  if(PCU_Comm_Self(PCUObj)==0)
    std::cout<<"Entered function\n"; 
  apf::MeshIterator* it = m->begin(1);
  apf::MeshEntity* edge;
  apf::Adjacent edgAdjVert;
  apf::Adjacent vertAdjEdg;
  //double gradingFactor = 1.3;
  double size[2];
  apf::Vector3 sizeVec;
  std::queue<apf::MeshEntity*> markedEdges;
  apf::MeshTag* isMarked = m->createIntTag("isMarked",1);
  apf::Field* size_scale = m->findField("proteus_size_scale");
 
  //marker structure for 0) not marked 1) marked 2)storage
  int marker[3] = {0,1,0}; 
 
  while((edge=m->iterate(it))){
    m->getAdjacent(edge, 0, edgAdjVert);
    for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){
      apf::getVector(size_scale,edgAdjVert[i],0,sizeVec);
      size[i]=sizeVec[idx]/sizeVec[0];
    }
    if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){
      //add edge to a queue 
      markedEdges.push(edge);
      //tag edge to indicate that it is part of queue 
      m->setIntTag(edge,isMarked,&marker[1]); 
 
    }
    else{
      m->setIntTag(edge,isMarked,&marker[0]); 
    }
  }
  m->end(it); 
 
  if(PCU_Comm_Self(PCUObj)==0)
    std::cout<<"Got queue of size "<<markedEdges.size()<<std::endl; 
  while(!markedEdges.empty()){
    edge = markedEdges.front();
    m->getAdjacent(edge, 0, edgAdjVert);
    for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){
      apf::getVector(size_scale,edgAdjVert[i],0,sizeVec);
      size[i]=sizeVec[idx]/sizeVec[0];
    }
    gradeSizeModify(m, gradingFactor, size, edgAdjVert, 
      vertAdjEdg, markedEdges, isMarked, apf::VECTOR, idx, 0, PCUObj);
    gradeSizeModify(m, gradingFactor, size, edgAdjVert, 
      vertAdjEdg, markedEdges, isMarked, apf::VECTOR, idx, 1, PCUObj);
 
    m->setIntTag(edge,isMarked,&marker[0]);
    markedEdges.pop();
  }
  it = m->begin(1);
  while((edge=m->iterate(it))){
    m->removeTag(edge,isMarked);
  }
  m->end(it); 
  m->destroyTag(isMarked);
  apf::synchronize(size_scale);
}