ExtractCriticalLayerFunctions.cpp File Reference

#include <cstdio>
#include <MultiRegions/ExpList.h>

Go to the source code of this file.

Functions
void	Computestreakpositions (MultiRegions::ExpListSharedPtr &streak, Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, NekDouble cr, NekDouble trans)

Function Documentation

Computestreakpositions()

void Computestreakpositions	(	MultiRegions::ExpListSharedPtr &	streak,
		Array< OneD, NekDouble > &	xc,
		Array< OneD, NekDouble > &	yc,
		NekDouble	cr,
		NekDouble	trans
	)

Definition at line 42 of file ExtractCriticalLayerFunctions.cpp.

References ASSERTL0, Nektar::MultiRegions::eY, Nektar::NekConstants::kNekUnsetDouble, Vmath::Vmax(), and Vmath::Vmin().

{
    int i;
    int npts = xc.num_elements();

    int nq = streak->GetTotPoints();    
    Array<OneD, NekDouble> derstreak(nq);
    Array<OneD, NekDouble> x(nq);
    Array<OneD, NekDouble> y(nq);    
    streak->GetCoords(x,y);         

    streak->BwdTrans(streak->GetCoeffs(),streak->UpdatePhys());
    streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), derstreak);

    // set intiial xc to be equispaced over mesh and yc to be zero 
    NekDouble x_max = Vmath::Vmax(nq,x,1);
    NekDouble x_min = Vmath::Vmin(nq,x,1);

    for(i = 0; i < npts; ++i)
    {
        xc[i]  = x_min + (x_max - x_min)*i/((NekDouble)(npts-1));
        yc[i] = 0.0;
    }
    
    
    int elmtid, offset,cnt;    
    NekDouble U,dU;
    NekDouble F;
    NekDouble ConvTol  = 1e-9;
    NekDouble CoordTol = 1e-5;
    int maxiter = 100;
    Array<OneD, NekDouble> coord(2);

    // Do Newton iteration on y direction
    cerr << "[";
    for(int e=0; e<npts; e++)
    {
        coord[0] = xc[e];
        coord[1] = yc[e];

        if(!(e%10))
        {
            cerr << ".";
        }

        F = 1000;
        cnt = 0;
        while((abs(F)> ConvTol)&&(cnt < maxiter))
        {
            elmtid = streak->GetExpIndex(coord,CoordTol);
            offset = streak->GetPhys_Offset(elmtid);

            U  = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
            dU = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);

            coord[1] = coord[1] - (U-cr)/dU;   
            
            F = U-cr;   
            cnt++;
        }
        ASSERTL0(cnt < maxiter, "Failed to converge Newton iteration");

        yc[e] = coord[1];
    }
    cerr << "]" << endl;

    if(trans != NekConstants::kNekUnsetDouble)
    {
        // output to interface file
        FILE *fp = fopen("interfacedat.geo","w");
        
        NekDouble y_max = Vmath::Vmax(nq,y,1);
        NekDouble y_min = Vmath::Vmin(nq,y,1);
        
        cnt = 1;
        fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
                cnt++,x_min,y_min);
        fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
                cnt++,x_max,y_min);
        fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
                cnt++,x_max,y_max);
        fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
                cnt++,x_min,y_max);
        
        for(i = 0; i < npts; ++i)
        {
            fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n",
                    cnt++,xc[i],yc[i]);
        }
        
        fclose(fp);

        
        // output to interface_up file as bend of vertical shift and 45 degrees shift
        fp = fopen("interfacedat_up.geo","w");
        
        
        NekDouble nx,ny,norm;
        
        fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n",cnt++,xc[0],yc[0]+trans);
        
        for(i = 1; i < npts-1; ++i)
        {
            norm = sqrt((xc[i+1]-xc[i-1])*(xc[i+1]-xc[i-1])+(yc[i+1]-yc[i-1])*(yc[i+1]-yc[i-1]));
            nx = (yc[i-1]-yc[i+1])/norm;
            ny = (xc[i+1]-xc[i-1])/norm;
            
            fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n",
                    cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
        }
        
        fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n",cnt++,xc[npts-1],yc[npts-1]+trans);
        
        
        // output to interface_up file as bend of vertical shift and 45 degrees shift
        fp = fopen("interfacedat_dn.geo","w");
        
        trans = -trans;
        
        fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n",cnt++,xc[0],yc[0]+trans);
        
        for(i = 1; i < npts-1; ++i)
        {
            norm = sqrt((xc[i+1]-xc[i-1])*(xc[i+1]-xc[i-1])+(yc[i+1]-yc[i-1])*(yc[i+1]-yc[i-1]));
            nx = (yc[i-1]-yc[i+1])/norm;
            ny = (xc[i+1]-xc[i-1])/norm;
            
            fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n",
                    cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
        }
        
        fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n",cnt++,xc[npts-1],yc[npts-1]+trans);
    }

}