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.

{
    int i;
    int npts = xc.size();
 
    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);
    }
}

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