MoveMeshToCriticalLayer.cpp File Reference

#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <tinyxml.h>
#include <SpatialDomains/MeshGraph.h>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList2D.h>
#include <LibUtilities/BasicUtils/ParseUtils.h>
#include "ExtractCriticalLayerFunctions.h"

Go to the source code of this file.

Classes
struct	MoveVerts

Enumerations
enum	SolveType { eSolveX, eSolveY, eSolveXY, eNoSolve }

Functions
void	GetInterfaceVerts (const int compositeID, SpatialDomains::MeshGraphSharedPtr &mesh, vector< int > &InterfaceVerts)
void	GetStreakLocation (LibUtilities::SessionReaderSharedPtr &vSession, SpatialDomains::MeshGraphSharedPtr &mesh, string &fieldfile, Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc)
void	GetNewVertexLocation (TiXmlElement *doc, SpatialDomains::MeshGraphSharedPtr &mesh, vector< int > &InterfaceVerts, Array< OneD, NekDouble > &xstreak, Array< OneD, NekDouble > &ystreak, Array< OneD, NekDouble > &vertx, Array< OneD, NekDouble > &verty, int maxiter)
void	TurnOffEdges (TiXmlElement *doc, SpatialDomains::SegGeomMap &meshedges, Array< OneD, MoveVerts > &verts)
void	RedefineVertices (TiXmlElement *doc, Array< OneD, NekDouble > &dvertx, Array< OneD, NekDouble > &dverty)
void	EnforceRotationalSymmetry (SpatialDomains::MeshGraphSharedPtr &mesh, Array< OneD, NekDouble > &dvertx, Array< OneD, NekDouble > &dverty)
int	main (int argc, char *argv[])

Enumeration Type Documentation

SolveType

enum SolveType

Enumerator
eSolveX
eSolveY
eSolveXY
eNoSolve

Definition at line 50 of file MoveMeshToCriticalLayer.cpp.

{
    eSolveX,
    eSolveY,
    eSolveXY,
    eNoSolve
};

Function Documentation

EnforceRotationalSymmetry()

void EnforceRotationalSymmetry	(	SpatialDomains::MeshGraphSharedPtr &	mesh,
		Array< OneD, NekDouble > &	dvertx,
		Array< OneD, NekDouble > &	dverty
	)

Definition at line 678 of file MoveMeshToCriticalLayer.cpp.

References Vmath::Vmax().

Referenced by main().

{
    int i,j;
    int nverts = mesh->GetNvertices();
    Array<OneD, NekDouble> x(nverts),y(nverts);
    NekDouble xval,yval,zval;

    for(i = 0; i < nverts; ++i)
    {
        mesh->GetVertex(i)->GetCoords(xval,yval,zval);
        x[i] = xval + dvertx[i];
        y[i] = yval + dverty[i];
    }
    
    NekDouble xmax = Vmath::Vmax(nverts,x,1);
    NekDouble tol = 1e-5, dist2,xrot,yrot;
    Array<OneD,int> index(nverts);
    // find nearest 
    for(i = 0; i < nverts; ++i)
    {
        xrot = -x[i] + xmax;
        yrot = -y[i];
        tol  = 1.0;

        for(j = 0; j < nverts; ++j)
        {
            dist2 = (x[j]-xrot)*(x[j]-xrot) + (y[j]-yrot)*(y[j]-yrot);
            if(dist2 < tol)
            {
                index[i]  = j;
                tol = dist2;
            }
        }
    }

    //average points and recalcualte dvertx, dverty
    for(i = 0; i < nverts; ++i)
    {
        mesh->GetVertex(i)->GetCoords(xval,yval,zval);

        xrot = 0.5*(-x[index[i]] + xmax + x[i]);
        yrot = 0.5*(-y[index[i]] + y[i]);
        
        dvertx[i] = xrot - xval;
        dverty[i] = yrot - yval;
    }
}

GetInterfaceVerts()

void GetInterfaceVerts	(	const int	compositeID,
		SpatialDomains::MeshGraphSharedPtr &	mesh,
		vector< int > &	InterfaceVerts
	)

Definition at line 160 of file MoveMeshToCriticalLayer.cpp.

Referenced by main().

{
    SpatialDomains::CompositeSharedPtr composite;
    composite = mesh->GetComposite(compositeID);
    int compsize = composite->m_geomVec.size();
    map<int,int> vertmap;

    for(int i = 0; i < compsize; ++i)
    {
        if(vertmap.count(composite->m_geomVec[i]->GetVid(0)) == 0)
        {
            InterfaceVerts.push_back(composite->m_geomVec[i]->GetVid(0)); 
            vertmap[composite->m_geomVec[i]->GetVid(0)]  = 1;
        }

        if(vertmap.count(composite->m_geomVec[i]->GetVid(1)) == 0)
        {
            InterfaceVerts.push_back(composite->m_geomVec[i]->GetVid(1)); 
            vertmap[composite->m_geomVec[i]->GetVid(1)]  = 1;
        }
    }
}

GetNewVertexLocation()

void GetNewVertexLocation	(	TiXmlElement *	doc,
		SpatialDomains::MeshGraphSharedPtr &	mesh,
		vector< int > &	InterfaceVerts,
		Array< OneD, NekDouble > &	xstreak,
		Array< OneD, NekDouble > &	ystreak,
		Array< OneD, NekDouble > &	vertx,
		Array< OneD, NekDouble > &	verty,
		int	maxiter
	)

Definition at line 221 of file MoveMeshToCriticalLayer.cpp.

References ASSERTL0, Nektar::SpatialDomains::PointGeom::dist(), eNoSolve, eSolveX, eSolveXY, eSolveY, and TurnOffEdges().

Referenced by main().

{
    int i,j,k;
    int nverts = mesh->GetNvertices(); 

    SpatialDomains::SegGeomMap meshedges = mesh->GetAllSegGeoms();

    Array<OneD,MoveVerts> Verts(nverts);
    
    // loop mesh edges and fill in verts info
    SpatialDomains::PointGeomSharedPtr v0,v1;
    SpatialDomains::PointGeom dist;

    int vid0,vid1;
    NekDouble kspring;
    NekDouble x,y,x1,y1,z1,x2,y2,z2;

    // Setup intiial spring and verts
    for(auto &segIter : meshedges)
    {
        vid0 = (segIter.second)->GetVid(0);
        vid1 = (segIter.second)->GetVid(1);

        v0 = (segIter.second)->GetVertex(0);
        v1 = (segIter.second)->GetVertex(1);
        
        kspring = 1.0/v0->dist(*v1); 
        
        Verts[vid0].kspring.push_back(kspring);
        Verts[vid0].springVid.push_back(vid1);
        Verts[vid1].kspring.push_back(kspring);
        Verts[vid1].springVid.push_back(vid0);

    }
    
    // Scale spring by total around vertex and turn on all vertices. 
    for(i = 0; i < nverts; ++i)
    {
        NekDouble invktot = 0.0;
        for(j = 0; j < Verts[i].kspring.size(); ++j)
        {
            invktot += Verts[i].kspring[j];
        }
        invktot = 1.0/invktot;
        for(j = 0; j < Verts[i].kspring.size(); ++j)
        {
            Verts[i].kspring[j] *= invktot;
        }

        Verts[i].solve = eSolveXY;
    }        


    // Turn off all edges defined by composite lists of correct dimension
    TurnOffEdges(doc,meshedges,Verts);
    
    NekDouble z,h0,h1,h2;
    // Set interface vertices to lie on critical layer
    for(i = 0; i < InterfaceVerts.size(); ++i)
    {
        Verts[InterfaceVerts[i]].solve = eNoSolve; 
        mesh->GetVertex(InterfaceVerts[i])->GetCoords(x,y,z);
        
        for(j = 0; j < xstreak.num_elements()-1; ++j)
        {
            if((x >= xstreak[j])&&(x <= xstreak[j+1]))
            {
                break;
            }
        }
        
        ASSERTL0(j != xstreak.num_elements(),"Did not find x location along critical layer");
        
        k = (j==0)?1: j; // offset index at beginning
        
        // quadraticalling interpolate points
        h0 = (x-xstreak[k])*(x-xstreak[k+1])/
            ((xstreak[k-1]-xstreak[k])*(xstreak[k-1]-xstreak[k+1]));
        h1 = (x-xstreak[k-1])*(x-xstreak[k+1])/
            ((xstreak[k]-xstreak[k-1])*(xstreak[k]-xstreak[k+1]));
        h2 = (x-xstreak[k-1])*(x-xstreak[k])/
            ((xstreak[k+1]-xstreak[k-1])*(xstreak[k+1]-xstreak[k]));
        
        dvertx[InterfaceVerts[i]] =  (xstreak[k-1]*h0 + xstreak[k]*h1 + xstreak[k+1]*h2) - x; 
        dverty[InterfaceVerts[i]] =  (ystreak[k-1]*h0 + ystreak[k]*h1 + ystreak[k+1]*h2) - y; 
    }

    // shift quads in critical layer to move more or less rigidly
    SpatialDomains::QuadGeomMap quadgeom = mesh->GetAllQuadGeoms();
    for(auto &quadIter : quadgeom)
    {
        for(i = 0; i < 4; ++i)
        {
            vid0 = (quadIter.second)->GetVid(i);
            
            switch(Verts[vid0].solve)
            {
            case eSolveXY:
                {
                    mesh->GetVertex(vid0)->GetCoords(x,y,z);

                    // find nearest interface vert
                    mesh->GetVertex(InterfaceVerts[0])->GetCoords(x1,y1,z1);
                    for(j = 0; j < InterfaceVerts.size()-1; ++j)
                    {
                        mesh->GetVertex(InterfaceVerts[j+1])->GetCoords(x2,y2,z2);
                        if((x >= x1)&&(x < x2))
                        {
                            break;
                        }
                        x1 = x2;
                        y1 = y2;
                    }
                    
                    // currently just shift vert as average of two sides
                    dvertx[vid0] = (x2-x)/(x2-x1)*dvertx[InterfaceVerts[j]]+
                        (x-x1)/(x2-x1)*dvertx[InterfaceVerts[j+1]];
                    dverty[vid0] = (x2-x)/(x2-x1)*dverty[InterfaceVerts[j]]+
                        (x-x1)/(x2-x1)*dverty[InterfaceVerts[j+1]];
                }
                break;
            case eSolveY:
                {
                    mesh->GetVertex(vid0)->GetCoords(x,y,z);
                    mesh->GetVertex(InterfaceVerts[0])->GetCoords(x1,y1,z1);
                    
                    if(fabs(x-x1) < 1e-6)
                    {
                        dverty[vid0] = dverty[InterfaceVerts[0]];
                    }
                    else
                    {
                        dverty[vid0] = dverty[InterfaceVerts[InterfaceVerts.size()-1]];
                    }
                }
                break;
            default:
                break;
            }
            Verts[vid0].solve = eNoSolve;
        }
    }

    


            
    // Iterate internal vertices 
    bool ContinueToIterate = true;
    int cnt  = 0;
    int nsum = 0;
    NekDouble dsum,dx,dy,sum,prev_sum = 0.0;
    NekDouble tol = 1e-3,fac;
    int blend = 50;

    while (ContinueToIterate)
    {
        
        sum = 0.0;
        nsum = 0;

        // use a ramping function to help move interior slowly 
        fac = (cnt < blend)? 1.0/(blend+1.0)*(cnt+1): 1.0;
        
        for(i = 0; i < nverts; ++i)
        {
            if(Verts[i].solve != eNoSolve)
            {
                dx = dy = 0.0;

                if((Verts[i].solve == eSolveX)||(Verts[i].solve == eSolveXY))
                {
                    for(j = 0; j < Verts[i].kspring.size(); ++j)
                    {
                        dx += fac*Verts[i].kspring[j]*dvertx[Verts[i].springVid[j]];
                    }
                }

                if((Verts[i].solve == eSolveY)||(Verts[i].solve == eSolveXY))
                {
                    for(j = 0; j < Verts[i].kspring.size(); ++j)
                    {
                        dy += fac*Verts[i].kspring[j]*dverty[Verts[i].springVid[j]];
                    }
                }
                
                dsum = (dx*dx + dy*dy);
                
                dvertx[i] = dx;
                dverty[i] = dy;
                
                if(dsum > 1e-16)
                {
                    //sum  += dsum/(dvertx[i]*dvertx[i] + dverty[i]*dverty[i]);
                    sum  += dsum;
                    nsum += 1;
                }
            }
        }

        if(nsum)
        {
            sum = sqrt(sum/(NekDouble)nsum);
            
            NekDouble chg = sum-prev_sum;
            prev_sum = sum;

            cerr << "Iteration " << cnt << " : " << chg << endl;

            if((chg < tol)&&(cnt > blend))
            {
                ContinueToIterate = false;
            }
            
        }
        else if(cnt > blend)
        {
            ContinueToIterate = false;

        }            

        if(cnt++ > maxiter)
        {
            ContinueToIterate = false;
        }
    }
}

GetStreakLocation()

void GetStreakLocation	(	LibUtilities::SessionReaderSharedPtr &	vSession,
		SpatialDomains::MeshGraphSharedPtr &	mesh,
		string &	fieldfile,
		Array< OneD, NekDouble > &	xc,
		Array< OneD, NekDouble > &	yc
	)

Definition at line 183 of file MoveMeshToCriticalLayer.cpp.

References Computestreakpositions(), and Nektar::LibUtilities::Import().

Referenced by main().

{
    //-------------------------------------------------------------
    // Define Streak Expansion   
    MultiRegions::ExpListSharedPtr streak;   
    
    streak = MemoryManager<MultiRegions::ExpList2D>
        ::AllocateSharedPtr(vSession,mesh);
    //---------------------------------------------------------------

    //----------------------------------------------
    // Import field file.
    vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
    vector<vector<NekDouble> > fielddata;
    LibUtilities::Import(fieldfile,fielddef,fielddata);
    //----------------------------------------------

    //----------------------------------------------
    // Copy data from field file
    string  streak_field("w");
    for(unsigned int i = 0; i < fielddata.size(); ++i)
    {
        streak->ExtractDataToCoeffs(fielddef [i],
                                    fielddata[i],
                                    streak_field,
                                    streak->UpdateCoeffs());
    }
    //----------------------------------------------
    
    NekDouble cr = 0.0;

    cerr << "Extracting Critical Layer ";
    Computestreakpositions(streak,xc, yc,cr);    
}

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 94 of file MoveMeshToCriticalLayer.cpp.

References EnforceRotationalSymmetry(), GetInterfaceVerts(), GetNewVertexLocation(), GetStreakLocation(), and RedefineVertices().

{
    if(argc !=4)
    {
        fprintf(stderr,"Usage: ./MoveMeshToCriticalLayer  meshfile streakfile  outfile\n");
        exit(1);
    }
    
    //------------------------------------------------------------
    // Create Session file by reading only first meshfile
    //-----------------------------------------------------------
    LibUtilities::SessionReaderSharedPtr vSession
        = LibUtilities::SessionReader::CreateInstance(2, argv);
    SpatialDomains::MeshGraphSharedPtr mesh =
        SpatialDomains::MeshGraph::Read(vSession);
    
    //-------------------------------------------------------------
    // Read in mesh from input file
    //------------------------------------------------------------
    TiXmlDocument& meshdoc = vSession->GetDocument();
    TiXmlHandle docHandle(&meshdoc);
    TiXmlElement* doc = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
    
    //------------------------------------------------------------
    // Get list of vertices along interface
    //------------------------------------------------------------
    vector<int> InterfaceVerts;
    GetInterfaceVerts(300,mesh,InterfaceVerts);

    //------------------------------------------------------------
    // Determine location of new interface vertices
    //------------------------------------------------------------
    Array<OneD,NekDouble> xstreak(50),ystreak(50);
    string fieldfile(argv[argc-2]);
    GetStreakLocation(vSession,mesh,fieldfile,xstreak,ystreak);

    //------------------------------------------------------------
    // Move internal mesh using critical layer info and under string analogy 
    //------------------------------------------------------------
    Array<OneD,NekDouble>  dvertx(mesh->GetNvertices(),0.0), dverty(mesh->GetNvertices(),0.0); 
    int maxiter;
    vSession->LoadParameter("MoveMeshMaxIterations",maxiter,100);

    GetNewVertexLocation(doc, mesh,InterfaceVerts,xstreak,ystreak,dvertx,dverty,maxiter);

    //------------------------------------------------------------
    // Enforce rotational symmetry on mesh 
    //------------------------------------------------------------
    if(vSession->DefinesSolverInfo("EnforceRotationalSymmetry"))
    {
        EnforceRotationalSymmetry(mesh,dvertx,dverty);
    }

    //------------------------------------------------------------
    // Redfine vertices in doc 
    //------------------------------------------------------------
    RedefineVertices(doc,dvertx,dverty);

    //------------------------------------------------------------
    // Write out moved mesh file 
    //------------------------------------------------------------
    std::string outfile(argv[argc-1]);
    meshdoc.SaveFile(outfile);
}

RedefineVertices()

void RedefineVertices	(	TiXmlElement *	doc,
		Array< OneD, NekDouble > &	dvertx,
		Array< OneD, NekDouble > &	dverty
	)

Error value returned by TinyXML.

Definition at line 601 of file MoveMeshToCriticalLayer.cpp.

References ASSERTL0.

Referenced by main().

{


    TiXmlElement* element = doc->FirstChildElement("VERTEX");
    ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
    
    TiXmlElement *vertex = element->FirstChildElement("V");
    
    int indx;
    int nextVertexNumber = -1;
    int err;    /// Error value returned by TinyXML.
    
    vector<NekDouble> xpts,ypts,zpts;
    NekDouble xval, yval, zval;
        
    while (vertex)
    {
        nextVertexNumber++;
        
        TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
        std::string attrName(vertexAttr->Name());

        ASSERTL0(attrName == "ID", (std::string("Unknown attribute name: ") + attrName).c_str());
        
        err = vertexAttr->QueryIntValue(&indx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
        
        // Now read body of vertex
        std::string vertexBodyStr;
        
        TiXmlNode *vertexBody = vertex->FirstChild();

        while (vertexBody)
        {
            // Accumulate all non-comment body data.
            if (vertexBody->Type() == TiXmlNode::TINYXML_TEXT)
            {
                vertexBodyStr += vertexBody->ToText()->Value();
                vertexBodyStr += " ";
            }
            
            vertexBody = vertexBody->NextSibling();
        }

        ASSERTL0(!vertexBodyStr.empty(), "Vertex definitions must contain vertex data.");

        // Get vertex data from the data string.
        std::istringstream vertexDataStrm(vertexBodyStr.c_str());

        try
        {
            while(!vertexDataStrm.fail())
            {
                vertexDataStrm >> xval >> yval >> zval;                
            }

            xval += dvertx[indx];
            yval += dverty[indx];
            
            stringstream s;
            s << scientific << setprecision(8) 
              << xval << " " << yval << " " << zval;

            vertex->ReplaceChild(vertex->FirstChild(), TiXmlText(s.str()));
        }
        catch(...)
        {
            ASSERTL0(false, "Unable to read VERTEX data.");
        }   
        vertex = vertex->NextSiblingElement("V");
    }

}

TurnOffEdges()

void TurnOffEdges	(	TiXmlElement *	doc,
		SpatialDomains::SegGeomMap &	meshedges,
		Array< OneD, MoveVerts > &	verts
	)

All elements are of the form: "<C ID = "N"> ... </C>". Read the ID field first.

Parse out the element components corresponding to type of element.

Keep looking

Definition at line 458 of file MoveMeshToCriticalLayer.cpp.

References ASSERTL0, eNoSolve, eSolveX, eSolveY, and NEKERROR.

Referenced by GetNewVertexLocation().

{
    TiXmlElement* field = doc->FirstChildElement("COMPOSITE");
    ASSERTL0(field, "Unable to find COMPOSITE tag in file.");

    int nextCompositeNumber = -1;

    /// All elements are of the form: "<C ID = "N"> ... </C>".
    /// Read the ID field first.
    TiXmlElement *composite = field->FirstChildElement("C");

    while (composite)
    {
        nextCompositeNumber++;
        
        int indx;
        int err = composite->QueryIntAttribute("ID", &indx);
        ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
        
        TiXmlNode* compositeChild = composite->FirstChild();
        // This is primarily to skip comments that may be present.
        // Comments appear as nodes just like elements.
        // We are specifically looking for text in the body
        // of the definition.
        while(compositeChild && compositeChild->Type() != TiXmlNode::TINYXML_TEXT)
        {
            compositeChild = compositeChild->NextSibling();
        }
        
        ASSERTL0(compositeChild, "Unable to read composite definition body.");
        std::string compositeStr = compositeChild->ToText()->ValueStr();
        
        /// Parse out the element components corresponding to type of element.
        std::istringstream compositeDataStrm(compositeStr.c_str());

        try
        {
            std::string compositeElementStr;
            compositeDataStrm >> compositeElementStr;
            
            std::istringstream tokenStream(compositeElementStr);
            char type;
            
            tokenStream >> type;
            
            // in what follows we are assuming there is only one block of data 
            std::string::size_type indxBeg = compositeElementStr.find_first_of('[') + 1;
            std::string::size_type indxEnd = compositeElementStr.find_last_of(']') - 1;
        
            ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") +  compositeElementStr).c_str());
            
            std::string indxStr = compositeElementStr.substr(indxBeg, indxEnd - indxBeg + 1);
            std::vector<unsigned int> seqVector;
            
            bool err = ParseUtils::GenerateSeqVector(indxStr, seqVector);
            
            ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());
            
            switch(type)
            {
            case 'E':  // Turn off vertices along composite edges 
                {
                    int seqlen = seqVector.size();
                    NekDouble x0,y0,z0,x1,y1,z1;
                    int vid0,vid1;

                    for(int i = 0; i < seqlen; ++i)
                    {
                        meshedges[seqVector[i]]->GetVertex(0)->GetCoords(x0,y0,z0);
                        meshedges[seqVector[i]]->GetVertex(1)->GetCoords(x1,y1,z1);
                        vid0 = meshedges[seqVector[i]]->GetVid(0);
                        vid1 = meshedges[seqVector[i]]->GetVid(1);

                        if(fabs(x0-x1) < 1e-8)
                        {
                            //identify corners by double visit
                            if(Verts[vid0].solve == eSolveX)
                            {
                                Verts[vid0].solve = eNoSolve;
                            }
                            else
                            {
                                Verts[vid0].solve = eSolveY;
                            }
                            
                            if(Verts[vid1].solve == eSolveX)
                            {
                                Verts[vid1].solve = eNoSolve;
                            }
                            else
                            {
                                Verts[vid1].solve = eSolveY;
                            }
                        }

                        if(fabs(y0-y1) < 1e-8)
                        {
                            //identify corners by double visit
                            if(Verts[vid0].solve == eSolveY)
                            {
                                Verts[vid0].solve = eNoSolve;
                            }
                            else
                            {
                                Verts[vid0].solve = eSolveX;
                            }
                            
                            if(Verts[vid1].solve == eSolveY)
                            {
                                Verts[vid1].solve = eNoSolve;
                            }
                            else
                            {
                                Verts[vid1].solve = eSolveX;
                            }
                        }
                    }

                }
                break;
                
                case 'T':  case 'Q':  // do nothing
                {
                    break;
                }
            default:
                NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + compositeElementStr).c_str());
            }
            
        }
        catch(...)
        {
            NEKERROR(ErrorUtil::efatal,
                     (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
        }
        
        /// Keep looking
        composite = composite->NextSiblingElement("C");
    }
}