ProcessWallNormalData.cpp

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////////
//
//  File: ProcessWallNormalData.cpp
//
//  For more information, please see: http://www.nektar.info/
//
//  The MIT License
//
//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
//  Department of Aeronautics, Imperial College London (UK), and Scientific
//  Computing and Imaging Institute, University of Utah (USA).
//
//  Permission is hereby granted, free of charge, to any person obtaining a
//  copy of this software and associated documentation files (the "Software"),
//  to deal in the Software without restriction, including without limitation
//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
//  and/or sell copies of the Software, and to permit persons to whom the
//  Software is furnished to do so, subject to the following conditions:
//
//  The above copyright notice and this permission notice shall be included
//  in all copies or substantial portions of the Software.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
//  DEALINGS IN THE SOFTWARE.
//
//  Description: Get the wall-normal data at a given origin.
//
////////////////////////////////////////////////////////////////////////////////
 
#include <iostream>
#include <string>
 
#include <FieldUtils/Interpolator.h>
#include <LibUtilities/BasicUtils/ParseUtils.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <MultiRegions/ExpList.h>
 
#include "ProcessWallNormalData.h"
 
using namespace std;
 
namespace Nektar::FieldUtils
{
 
ModuleKey ProcessWallNormalData::className =
    GetModuleFactory().RegisterCreatorFunction(
        ModuleKey(eProcessModule, "wallNormalData"),
        ProcessWallNormalData::create,
        "Export data in wall-normal direction from a single point on the "
        "wall.");
 
ProcessWallNormalData::ProcessWallNormalData(FieldSharedPtr f)
    : ProcessBoundaryExtract(f)
{
    f->m_writeBndFld = false; // turned on in upstream ProcessBoundaryExtract
 
    m_config["xorig"] =
        ConfigOption(false, "0.5,0.0,0.0",
                     "The point to be projected onto the wall to get the \
                          sampling origin. default=[0.5,0,0]");
    m_config["projDir"] =
        ConfigOption(false, "0.0,1.0,0.0",
                     "The direction to project the point onto the wall to \
                          get the sampling origin. default=[0,1,0]");
    m_config["maxDist"] = ConfigOption(
        false, "1.0", "Distance to limit projection distance to find the \
                          desired sampling origin. defalut=1.0");
    m_config["distH"] = ConfigOption(
        false, "0.01", "Sampling distance along the wall normal at the \
                          sampling origin. default=0.1");
    m_config["nptsH"] = ConfigOption(
        false, "5", "Number of sampling points along the wall normal. \
                          default=5");
    m_config["d"] = ConfigOption(
        false, "0.1", "The parameter that controls the sampling points' \
                          distribution. d should be in the range (0,inf). d \
                          in (0,0.95] gives controled points; d in (0.95,inf) \
                          gives evenly spaced points");
 
    // To allow some functions called somewhere else outside Process(),
    // m_spacedim should be defined while it cannot be defined here in the
    // constructor because f->m_exp[0]->GetCoordim(0) and
    // f->m_graph->GetMeshDimension() has not been defeind yet
    // *It is defeind in m_f->SetUpExp(vm), where vm is necessary.
}
 
ProcessWallNormalData::~ProcessWallNormalData()
{
}
 
/*Note
 * This module is used to get field data in the wall-normal direction.
 * The input cases can be 2D, 2.5D and 3D.
 * The data will be exported with .pts extension.
 *
 * The user defined parameters are: bnd, xorig, searchDir, maxDist, distH,
 * nptsH, and d.
 *  - bnd=0 is the boundary id. This boundary should contain the desired origin.
 *  - xorig="x,y,z" are the coordinates of the input sampling origin. They are
 *    used as the references to get exact origin on the wall.
 *  - projDir="0,1,0" is the projection direction to find the point on the wall
 *  - maxDist=1.0 is the distance that constrains the projection on the wall.
 *    The projection distance should be smaller than this distance. This
 *    parameter is set in case the boundary is wavey so that multiple
 *    projections exist at the same time.
 *  - distH=0.01 is the sampling depth in the wall-normal direction.
 *  - nptsH=11 is the number of sampling points along wall-normal direction.
 *  - d=0.1 is a destribution control parameter of the sampling points. It
 *    should be in the range (0,inf). d in range (0,0.95] for controlled array.
 *    d>0.95 gives evenly spaced array.
 */
void ProcessWallNormalData::v_Process(po::variables_map &vm)
{
    ProcessBoundaryExtract::v_Process(vm);
 
    // Get dim to store data
    const int nfields       = m_f->m_variables.size();
    const int nCoordDim     = m_f->m_exp[0]->GetCoordim(0);
    const int nBndLcoordDim = nCoordDim - 1;
    m_spacedim        = m_f->m_exp[0]->GetCoordim(0) + m_f->m_numHomogeneousDir;
    const int totVars = m_spacedim + nfields;
 
    // Initialize the sampling parameters
    std::vector<NekDouble> xorig, searchDir;
    ASSERTL0(ParseUtils::GenerateVector(m_config["xorig"].as<string>(), xorig),
             "Failed to interpret origin coordinates string");
    ASSERTL0(
        ParseUtils::GenerateVector(m_config["projDir"].as<string>(), searchDir),
        "Failed to interpret search direction string");
    const NekDouble maxDist = m_config["maxDist"].as<NekDouble>();
    const NekDouble distH   = m_config["distH"].as<NekDouble>();
    const int nptsH         = m_config["nptsH"].as<int>();
    const NekDouble delta   = m_config["d"].as<NekDouble>();
 
    Array<OneD, NekDouble> orig(3), projDir(3); // gloCoord of the origin
    for (int i = 0; i < 3; ++i)
    {
        orig[i]    = (xorig.size() >= (i + 1)) ? xorig[i] : 0.0;
        projDir[i] = (searchDir.size() >= (i + 1)) ? searchDir[i] : 0.0;
    }
    if (nCoordDim == 2)
    {
        projDir[2] = 0.0;
    }
    Vmath::Smul(3, 1.0 / sqrt(Vmath::Dot(3, projDir, 1, projDir, 1)), projDir,
                1, projDir, 1);
 
    // Update z according to the closest plane in 2.5D cases
    if (m_f->m_numHomogeneousDir == 1)
    {
        int nPlanes    = m_f->m_exp[0]->GetHomogeneousBasis()->GetZ().size();
        NekDouble lHom = m_f->m_exp[0]->GetHomoLen();
        if (orig[2] < 0.0 || orig[2] > lHom)
        {
            orig[2] = 0.0;
        }
        else
        {
            NekDouble dZ = lHom / nPlanes;
            NekDouble zTmp, zCur = 0.0, distTmp, distCur = 999.0;
            for (int i = 0; i <= nPlanes; ++i)
            {
                zTmp    = dZ * i;
                distTmp = fabs(orig[2] - zTmp);
                if (distTmp < distCur)
                {
                    distCur = distTmp;
                    zCur    = zTmp;
                    if (i == nPlanes)
                    {
                        zCur = 0.0;
                    }
                }
            }
            orig[2] = zCur;
        }
    }
 
    // Get the bnd id
    SpatialDomains::BoundaryConditions bcs(m_f->m_session,
                                           m_f->m_exp[0]->GetGraph());
    const SpatialDomains::BoundaryRegionCollection bregions =
        bcs.GetBoundaryRegions();
    map<int, int> BndRegionMap;
    int cnt = 0;
    for (auto &breg_it : bregions)
    {
        BndRegionMap[breg_it.first] = cnt++;
    }
    int bnd = BndRegionMap[m_f->m_bndRegionsToWrite[0]];
 
    // Get expansion list for boundary
    Array<OneD, MultiRegions::ExpListSharedPtr> BndExp(nfields);
    for (int i = 0; i < nfields; ++i)
    {
        BndExp[i] = m_f->m_exp[i]->UpdateBndCondExpansion(bnd);
    }
 
    //-------------------------------------------------------------------------
    // Find the element that contains the origin, and give the precise origin
    SpatialDomains::GeometrySharedPtr bndGeom;
    Array<OneD, NekDouble> locCoord(nBndLcoordDim, -999.0);
    NekDouble projDist;
    int elmtid;
    bool isInside = false;
 
    // Search and get precise locCoord
    const int nElmts = BndExp[0]->GetNumElmts();
    for (elmtid = 0; elmtid < nElmts; ++elmtid) // nElmts
    {
        bndGeom  = BndExp[0]->GetExp(elmtid)->GetGeom();
        isInside = BndElmtContainsPoint(bndGeom, orig, projDir, locCoord,
                                        projDist, maxDist);
        if (isInside)
        {
            break;
        }
    }
    ASSERTL0(isInside, "Failed to find the sampling origin on the boundary.");
 
    // Then Update the precise sampling position
    StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
    const int npts                            = bndXmap->GetTotPoints();
    Array<OneD, Array<OneD, NekDouble>> pts(nCoordDim);
    Array<OneD, Array<OneD, const NekDouble>> bndCoeffs(nCoordDim);
 
    for (int i = 0; i < nCoordDim; ++i)
    {
        pts[i]       = Array<OneD, NekDouble>(npts);
        bndCoeffs[i] = bndGeom->GetCoeffs(i); // 0/1/2 for x/y/z
        bndXmap->BwdTrans(bndCoeffs[i], pts[i]);
        orig[i] = bndXmap->PhysEvaluate(locCoord, pts[i]); // Update wall point
    }
 
    // Get outward-pointing normal vectors for all quadrature points on bnd
    // Use these normals as diretion references
    Array<OneD, Array<OneD, NekDouble>> normalsQ;
    m_f->m_exp[0]->GetBoundaryNormals(bnd, normalsQ);
 
    // Get the normals in the element
    // Set point key to get point id
    Array<OneD, LibUtilities::PointsKey> from_key(nBndLcoordDim);
    int from_nPtsPerElmt = 1;
    for (int i = 0; i < nBndLcoordDim; ++i)
    {
        from_key[i] = BndExp[0]->GetExp(elmtid)->GetBasis(i)->GetPointsKey();
        from_nPtsPerElmt *= from_key[i].GetNumPoints();
    }
 
    // Get ref normal
    Array<OneD, NekDouble> normalsRef(3, 0.0);
    int refId = elmtid * from_nPtsPerElmt + 0; // the 1st normal in the bnd elmt
    for (int i = 0; i < m_spacedim; ++i)
    {
        normalsRef[i] = normalsQ[i][refId];
    }
 
    // Get the precise normals and correct the direction to be inward
    Array<OneD, NekDouble> normals(3, 0.0);
    GetNormals(bndGeom, locCoord, normals);
 
    if (Vmath::Dot(3, normals, normalsRef) > 0.0)
    {
        Vmath::Neg(3, normals, 1);
    }
 
    // Output the info if -v is set
    if (m_f->m_verbose)
    {
        cout << "------ wallNormalData module ------\n";
        cout << "Input point:\n";
        cout << "  - [Px,Py,Pz] = [" << xorig[0] << ", " << xorig[1];
        if (xorig.size() >= 3)
        {
            cout << ", " << xorig[2] << "]\n";
        }
        else
        {
            cout << ", " << 0.0 << "]\n";
        }
        cout << "Projection direction:\n";
        cout << "  - [vx,vy,vz] = [" << projDir[0] << ", " << projDir[1] << ", "
             << projDir[2] << "]\n";
        cout << "Sampling origin on the wall:\n";
        cout << "  - [Ox,Oy,Oz] = [" << orig[0] << ", " << orig[1] << ", "
             << orig[2] << "]\n";
        cout << "Normals at the origin:\n";
        cout << "  - [nx,ny,nz] = [" << normals[0] << ", " << normals[1] << ", "
             << normals[2] << "]\n";
        cout
            << "Ref normals (at quadrature points in the projected element):\n";
        for (int i = 0; i < from_nPtsPerElmt; ++i)
        {
            cout << "  - " << i << ": [nx,ny,nz] = ["
                 << -normalsQ[0][elmtid * from_nPtsPerElmt + i] << ", "
                 << -normalsQ[1][elmtid * from_nPtsPerElmt + i];
            if (m_spacedim == 3)
            {
                cout << ", " << -normalsQ[2][elmtid * from_nPtsPerElmt + i]
                     << "]\n";
            }
            else
            {
                cout << "]\n";
            }
        }
        cout << endl;
    }
 
    //-------------------------------------------------------------------------
    // Set the depth of sampling array
    Array<OneD, NekDouble> h(nptsH, 0.0);
    if (delta > 0 && delta <= 0.95)
    {
        // Use expression in Agrawal's paper:
        // h = 1-tanh((1-ksi)*atanh(sqrt(1-delta)))/sqrt(1-delta), ksi in [0,1]
        NekDouble tmp1;
        const NekDouble tmp2 = 1.0 / (static_cast<NekDouble>(nptsH) - 1.0);
        const NekDouble tmp3 = sqrt(1.0 - delta);
        const NekDouble tmp4 = atanh(tmp3);
        const NekDouble tmp5 = 1.0 / tmp3;
        for (int i = 1; i < nptsH; ++i)
        {
            tmp1 = 1.0 - i * tmp2; // tmp1 = 1-ksi
            h[i] = 1 - tanh(tmp1 * tmp4) * tmp5;
            h[i] *= distH; // physical distance in normal direction
        }
    }
    else if (delta > 0.95)
    {
        // Evenly spaced array
        const NekDouble tmp1 = 1.0 / (static_cast<NekDouble>(nptsH) - 1.0);
        for (int i = 1; i < nptsH; ++i)
        {
            h[i] = i * tmp1;
            h[i] *= distH; // physical distance in normal direction
        }
    }
    else
    {
        ASSERTL0(false, "Input error. Delta needs to be greater than 0.0.");
    }
 
    // Set pts coordinates and interpoate the data
    Array<OneD, Array<OneD, NekDouble>> ptsH(totVars);
    for (int i = 0; i < totVars; ++i)
    {
        ptsH[i] = Array<OneD, NekDouble>(nptsH, 0.0);
    }
 
    for (int i = 0; i < m_spacedim; ++i)
    {
        for (int j = 0; j < nptsH; ++j)
        {
            ptsH[i][j] = orig[i] + h[j] * normals[i]; // x0+dist*nx
        }
    }
 
    m_f->m_fieldPts = MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(
        m_spacedim, m_f->m_variables, ptsH, LibUtilities::NullPtsInfoMap);
 
    Interpolator<std::vector<MultiRegions::ExpListSharedPtr>> interp;
    interp.Interpolate(m_f->m_exp, m_f->m_fieldPts,
                       NekConstants::kNekUnsetDouble);
}
 
/**
 * @brief Project a single point along the given direction to a plane
 * @param gloCoord     Global coordinate of the point. size=3.
 * @param projDir      Projection direction, which is also the normal vector of
 *                     the target plane. size=3, norm=1.
 * @param distToOrig   The distance from the origin (0,0,0) to the target plane
 * @param projGloCoord The global coordinate of the projecion result.
 */
void ProcessWallNormalData::ProjectPoint(
    const Array<OneD, const NekDouble> &gloCoord,
    const Array<OneD, const NekDouble> &projDir, const NekDouble distToOrig,
    Array<OneD, NekDouble> &projGloCoord)
{
    NekDouble tmp1;
    Array<OneD, NekDouble> tmp2(3);
 
    tmp1 = Vmath::Dot(3, gloCoord, 1, projDir, 1); // |xn| = x dot n
    Vmath::Smul(3, distToOrig - tmp1, projDir, 1, tmp2,
                1); // d_xn = (dist-|xn|) n
    Vmath::Vadd(3, gloCoord, 1, tmp2, 1, projGloCoord, 1); // x' = x + d_xn
}
 
/**
 * @brief Project the vertices of the elmt along the given direction to a plane
 * @param pts          Global coordinate of the vertices of the elmt. size=2/3.
 * @param projDir      Projection direction, which is also the normal vector of
 *                     the target plane. size=3, norm=1.
 * @param distToOrig   The distance from the origin (0,0,0) to the target plane.
 * @param projPts      The global coordinate of the projecion result.
 */
void ProcessWallNormalData::ProjectVertices(
    const Array<OneD, const Array<OneD, NekDouble>> &pts,
    const Array<OneD, const NekDouble> &projDir, const NekDouble distToOrig,
    Array<OneD, Array<OneD, NekDouble>> &projPts)
{
    const int nCoordDim = pts.size(); // 2 for 2.5D cases
    const int npts      = pts[0].size();
 
    NekDouble tmp1;
    Array<OneD, NekDouble> singlePnt(nCoordDim), tmp2(nCoordDim);
 
    for (int i = 0; i < npts; ++i)
    {
        // Get a point
        for (int j = 0; j < nCoordDim; ++j)
        {
            singlePnt[j] = pts[j][i];
        }
 
        // Projection
        tmp1 = Vmath::Dot(nCoordDim, singlePnt, 1, projDir, 1);
        Vmath::Smul(nCoordDim, distToOrig - tmp1, projDir, 1, tmp2, 1);
        Vmath::Vadd(nCoordDim, singlePnt, 1, tmp2, 1, singlePnt, 1);
 
        // Save to projPts
        for (int j = 0; j < nCoordDim; ++j)
        {
            projPts[j][i] = singlePnt[j];
        }
    }
}
 
/**
 * @brief Determine if the projected point is inside the projected element.
 * @param projGloCoord The global coordinate of the projected single point.
 * @param projPts      The global coordinate of the projected vertices,size=2/3
 * @param projDir      Projection direction, which is used as the reference
 *                     direction in the 3D routine. size=3, norm=1.
 * @param paralTol     Tolerence to check if two vectors are parallel.
 * @param angleTol     Tolerence to check if the total angle is 2*pi.
 * @return             Inside (true) or not (false)
 */
bool ProcessWallNormalData::isInProjectedArea2D(
    const Array<OneD, const NekDouble> &projGloCoord,
    const Array<OneD, const Array<OneD, NekDouble>> &projPts,
    const NekDouble paralTol)
{
    const NekDouble npts = projPts[0].size();
 
    Array<OneD, NekDouble> vec1(2, 0.0), vec2(2, 0.0);
 
    vec1[0] = projPts[0][0] - projGloCoord[0];
    vec1[1] = projPts[1][0] - projGloCoord[1];
    vec2[0] = projPts[0][npts - 1] - projGloCoord[0];
    vec2[1] = projPts[1][npts - 1] - projGloCoord[1];
 
    Vmath::Smul(2, 1.0 / sqrt(Vmath::Dot(2, vec1, 1, vec1, 1)), vec1, 1, vec1,
                1);
    Vmath::Smul(2, 1.0 / sqrt(Vmath::Dot(2, vec2, 1, vec2, 1)), vec2, 1, vec2,
                1);
 
    if ((fabs(vec1[0] + vec2[0]) + fabs(vec1[1] + vec2[1])) < paralTol)
    {
        // On the line segement, true
        return true;
    }
    else
    {
        return false;
    }
}
 
bool ProcessWallNormalData::isInProjectedArea3D(
    const Array<OneD, const NekDouble> &projGloCoord,
    const Array<OneD, const Array<OneD, NekDouble>> &projPts,
    const Array<OneD, const NekDouble> &projDir, const NekDouble paralTol,
    const NekDouble angleTol)
{
 
    // Generate a polygen (quad), re-order the points on the edge
    const NekDouble npts = projPts[0].size();
    const int nptsEdge =
        sqrt(npts); // num of points on edge for geo representation
    const int nptsPolygon = 4 * nptsEdge - 4;
    Array<OneD, Array<OneD, NekDouble>> ptsPolygon(3);
    for (int i = 0; i < 3; ++i)
    {
        ptsPolygon[i] = Array<OneD, NekDouble>(nptsPolygon);
 
        for (int j = 0; j < nptsEdge; ++j)
        {
            ptsPolygon[i][j] = projPts[i][j];
        }
        for (int j = 0; j < nptsEdge - 2; ++j)
        {
            ptsPolygon[i][nptsEdge + j] = projPts[i][(j + 2) * nptsEdge - 1];
        }
        for (int j = 0; j < nptsEdge; ++j)
        {
            ptsPolygon[i][2 * nptsEdge - 2 + j] = projPts[i][npts - 1 - j];
        }
        for (int j = 0; j < nptsEdge - 2; ++j)
        {
            ptsPolygon[i][3 * nptsEdge - 2 + j] =
                projPts[i][nptsEdge * (nptsEdge - j - 2)];
        }
    }
 
    // Determine relation using angle method (sum=2*pi)
    NekDouble angleCos, angleAbs, angleSign, angleSum = 0.0;
    Array<OneD, NekDouble> vec1(3, 0.0), vec2(3, 0.0), vec3(3, 0.0);
    int id1, id2;
 
    for (int i = 0; i < nptsPolygon; ++i)
    {
        id1 = i;
        id2 = (id1 == (nptsPolygon - 1)) ? 0 : (id1 + 1);
 
        for (int j = 0; j < 3; ++j)
        {
            vec1[j] = ptsPolygon[j][id1] - projGloCoord[j];
            vec2[j] = ptsPolygon[j][id2] - projGloCoord[j];
        }
 
        Vmath::Smul(3, 1.0 / sqrt(Vmath::Dot(3, vec1, 1, vec1, 1)), vec1, 1,
                    vec1, 1);
        Vmath::Smul(3, 1.0 / sqrt(Vmath::Dot(3, vec2, 1, vec2, 1)), vec2, 1,
                    vec2, 1);
 
        if ((fabs(vec1[0] + vec2[0]) + fabs(vec1[1] + vec2[1]) +
             fabs(vec1[2] + vec2[2])) < paralTol)
        {
            // On the line segement, true
            // This branch is used to aviod angle=pi but angleSign=-1 (not 1)
            return true;
        }
        else
        {
            // Off the line segement, compute angle
            // vec3 = vec1 x vec2, not scaled
            vec3[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
            vec3[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
            vec3[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
 
            // Use the projDir as reference direction (positive)
            angleSign = Vmath::Dot(3, vec3, 1, projDir, 1) > 0.0 ? 1.0 : -1.0;
            angleCos  = Vmath::Dot(3, vec1, 1, vec2, 1);
 
            // Avoid error in using acos()
            if (angleCos > 1.0)
            {
                angleCos = 1.0;
            }
            else if (angleCos < -1.0)
            {
                angleCos = -1.0;
            }
 
            angleAbs = acos(angleCos);
            angleSum += angleSign * angleAbs;
        }
    }
 
    // Check
    angleSum = fabs(angleSum);
    if (fabs(angleSum - 2.0 * M_PI) < angleTol)
    {
        return true;
    }
    else
    {
        return false;
    }
}
 
/**
 * @brief Use iteration to get the locCoord. This routine should be used after
 *        we have checked the projected point is inside the projected element.
 * @param bndGeom      Geometry to get the xmap.
 * @param gloCoord     Global coordinate of the point. size=3.
 * @param pts          Global coordinate of the vertices of the elmt. size=2/3.
 * @param dieUse       The main direction(s) used to compute local coordinate
 * @param locCoord     Iteration results for local coordinate(s)
 * @param dist         Returned distance in physical space if the collapsed
 *                     locCoord is out of range [-1,1].
 * @param iterTol      Tolerence for iteration.
 * @param iterMax      Maximum iteration steps
 * @return             Converged (true) or not (false)
 */
bool ProcessWallNormalData::BisectionForLocCoordOnBndElmt(
    SpatialDomains::GeometrySharedPtr bndGeom,
    const Array<OneD, const NekDouble> &gloCoord,
    const Array<OneD, const Array<OneD, NekDouble>> &pts,
    const Array<OneD, const int> &dirUse, Array<OneD, NekDouble> &locCoord,
    const NekDouble iterTol, const int iterMax)
{
    // Initial settings
    Array<OneD, NekDouble> etaLR(2); // range [-1,1]
    etaLR[0] = -1.0;                 // left
    etaLR[1] = 1.0;                  // right
    NekDouble tmpL, tmpR;            // tmp values for L/R
 
    StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
    locCoord[0]                               = -2.0;
    int cnt                                   = 0;
    bool isConverge                           = false;
    while (cnt < iterMax)
    {
        tmpL = bndXmap->PhysEvaluate(etaLR, pts[dirUse[0]]);
        tmpR = bndXmap->PhysEvaluate(etaLR + 1, pts[dirUse[0]]);
 
        if (fabs(gloCoord[dirUse[0]] - tmpL) >=
            fabs(gloCoord[dirUse[0]] - tmpR))
        {
            etaLR[0] = 0.5 * (etaLR[0] + etaLR[1]);
        }
        else
        {
            etaLR[1] = 0.5 * (etaLR[0] + etaLR[1]);
        }
 
        if ((etaLR[1] - etaLR[0]) < iterTol)
        {
            locCoord[0] = 0.5 * (etaLR[0] + etaLR[1]);
            isConverge  = true;
            break;
        }
 
        ++cnt;
    }
 
    // Warning if failed
    if (cnt >= iterMax)
    {
        WARNINGL1(false, "Bisection iteration is not converged");
    }
 
    return isConverge;
}
 
bool ProcessWallNormalData::NewtonIterForLocCoordOnBndElmt(
    SpatialDomains::GeometrySharedPtr bndGeom,
    const Array<OneD, const NekDouble> &gloCoord,
    const Array<OneD, const Array<OneD, NekDouble>> &pts,
    const Array<OneD, const int> &dirUse, Array<OneD, NekDouble> &locCoord,
    NekDouble &dist, const NekDouble iterTol, const int iterMax)
{
 
    const NekDouble LcoordDiv = 15.0;
 
    StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
 
    Array<OneD, const NekDouble> Jac =
        bndGeom->GetMetricInfo()->GetJac(bndXmap->GetPointsKeys());
    NekDouble scaledTol =
        Vmath::Vsum(Jac.size(), Jac, 1) / ((NekDouble)Jac.size());
    scaledTol *= iterTol;
 
    // Set the gloCoord used to compute locCoord
    const int dir1 = dirUse[0];
    const int dir2 = dirUse[1];
 
    Array<OneD, NekDouble> Dx1D1(pts[dir1].size());
    Array<OneD, NekDouble> Dx1D2(pts[dir1].size());
    Array<OneD, NekDouble> Dx2D1(pts[dir1].size());
    Array<OneD, NekDouble> Dx2D2(pts[dir1].size());
 
    // Ideally this will be stored in m_geomfactors
    bndXmap->PhysDeriv(pts[dir1], Dx1D1, Dx1D2);
    bndXmap->PhysDeriv(pts[dir2], Dx2D1, Dx2D2);
 
    // Initial the locCoord, in [-1,1]
    locCoord[0] = 0.0;
    locCoord[1] = 0.0;
 
    NekDouble x1map, x2map, F1, F2;
    NekDouble derx1_1, derx1_2, derx2_1, derx2_2, jac;
    NekDouble resid;
    int cnt         = 0;
    bool isConverge = false;
 
    F1 = F2 = 2000; // Starting value of Function
    resid   = sqrt(F1 * F1 + F2 * F2);
    while (cnt++ < iterMax)
    {
        x1map = bndXmap->PhysEvaluate(locCoord, pts[dir1]);
        x2map = bndXmap->PhysEvaluate(locCoord, pts[dir2]);
 
        F1 = gloCoord[dir1] - x1map;
        F2 = gloCoord[dir2] - x2map;
 
        if (F1 * F1 + F2 * F2 < scaledTol)
        {
            resid      = sqrt(F1 * F1 + F2 * F2);
            isConverge = true;
            break;
        }
 
        // Interpolate derivative metric at locCoord
        derx1_1 = bndXmap->PhysEvaluate(locCoord, Dx1D1);
        derx1_2 = bndXmap->PhysEvaluate(locCoord, Dx1D2);
        derx2_1 = bndXmap->PhysEvaluate(locCoord, Dx2D1);
        derx2_2 = bndXmap->PhysEvaluate(locCoord, Dx2D2);
 
        jac = derx2_2 * derx1_1 - derx2_1 * derx1_2;
 
        // Use analytical inverse of derivitives which are
        // also similar to those of metric factors.
        locCoord[0] = locCoord[0] + (derx2_2 * (gloCoord[dir1] - x1map) -
                                     derx1_2 * (gloCoord[dir2] - x2map)) /
                                        jac;
 
        locCoord[1] = locCoord[1] + (-derx2_1 * (gloCoord[dir1] - x1map) +
                                     derx1_1 * (gloCoord[dir2] - x2map)) /
                                        jac;
 
        // locCoord have diverged so stop iteration
        if (!(std::isfinite(locCoord[0]) && std::isfinite(locCoord[1])))
        {
            dist = 1e16;
            std::ostringstream ss;
            ss << "nan or inf found in NewtonIterForLocCoordOnProjBndElmt in "
                  "element "
               << bndGeom->GetGlobalID();
            WARNINGL1(false, ss.str());
            return false;
        }
        if (fabs(locCoord[0]) > LcoordDiv || fabs(locCoord[1]) > LcoordDiv)
        {
            break;
        }
    }
 
    // Check distance for collapsed coordinate
    Array<OneD, NekDouble> eta(2);
    bndXmap->LocCoordToLocCollapsed(locCoord, eta);
 
    if (bndGeom->ClampLocCoords(eta, 0.0))
    {
        // calculate the global point corresponding to locCoord
        x1map = bndXmap->PhysEvaluate(eta, pts[dir1]);
        x2map = bndXmap->PhysEvaluate(eta, pts[dir2]);
 
        F1 = gloCoord[dir1] - x1map;
        F2 = gloCoord[dir2] - x2map;
 
        dist = sqrt(F1 * F1 + F2 * F2);
    }
    else
    {
        dist = 0.0;
    }
 
    // Warning if failed
    if (cnt >= iterMax)
    {
        Array<OneD, NekDouble> collCoords(2);
        bndXmap->LocCoordToLocCollapsed(locCoord, collCoords);
 
        // if coordinate is inside element dump error!
        if ((collCoords[0] >= -1.0 && collCoords[0] <= 1.0) &&
            (collCoords[1] >= -1.0 && collCoords[1] <= 1.0))
        {
            std::ostringstream ss;
 
            ss << "Reached MaxIterations (" << iterMax
               << ") in Newton iteration ";
            ss << "Init value (" << setprecision(4) << 0 << "," << 0 << ","
               << ") ";
            ss << "Fin  value (" << locCoord[0] << "," << locCoord[1] << ","
               << ") ";
            ss << "Resid = " << resid << " Tolerance = " << sqrt(scaledTol);
 
            WARNINGL1(cnt < iterMax, ss.str());
        }
    }
 
    return isConverge;
}
 
/**
 * @brief Check if a point can be projected onto an oundary element in a given
 *        direction. If yes, give the local coordinates of the projected point.
 *        we have checked the projected point is inside the projected element.
 * @param bndGeom      Pointer to the geometry of the boundary element.
 * @param gloCoord     Global coordinate of the point. size=3.
 * @param projDir      Projection direction, which is used as the reference
 *                     direction in the 3D routine. size=3, norm=1.
 * @param locCoord     Iteration results for local coordinates (if inside).
 * @param projDist     Projection distance betweem the point to the wall point.
 * @param maxDist      Disntance to check if the wall point is desired.
 * @param iterTol      Tolerence for iteration.
 * @return             Inside (true) or not (false)
 */
bool ProcessWallNormalData::BndElmtContainsPoint(
    SpatialDomains::GeometrySharedPtr bndGeom,
    const Array<OneD, const NekDouble> &gloCoord,
    const Array<OneD, const NekDouble> &projDir,
    Array<OneD, NekDouble> &locCoord, NekDouble &projDist,
    const NekDouble maxDist, const NekDouble iterTol)
{
    // Get variables
    StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
    const int npts                            = bndXmap->GetTotPoints();
    const int nCoordDim = m_f->m_exp[0]->GetCoordim(0); // =2 for 2.5D cases
 
    Array<OneD, Array<OneD, const NekDouble>> bndCoeffs(nCoordDim);
    Array<OneD, Array<OneD, NekDouble>> pts(nCoordDim);
    Array<OneD, Array<OneD, NekDouble>> projPts(nCoordDim);
    Array<OneD, NekDouble> projGloCoord(3, 0.0);
 
    for (int i = 0; i < nCoordDim; ++i)
    {
        pts[i]       = Array<OneD, NekDouble>(npts);
        projPts[i]   = Array<OneD, NekDouble>(npts);
        bndCoeffs[i] = bndGeom->GetCoeffs(i); // 0/1/2 for x/y/z
        bndXmap->BwdTrans(bndCoeffs[i], pts[i]);
    }
 
    // Project the point and vertices of the element in the input direction
    ProjectPoint(gloCoord, projDir, 0.0, projGloCoord);
    ProjectVertices(pts, projDir, 0.0, projPts);
 
    // Set the main direction(s) and the minor direction
    // The gloCoord for minor direction will not be used for locCoord iteration
    // dirUse[0] is the main dir for 2D/2.5D cases, dirUse[1]/[2] is the minor
    // one dirUse[0] and dirUse[1] are the main dir for 3D cases, dirUse[2] is
    // hte minor one
    int dirMaxId = 0; // id to get the dir with largest projDir component
    for (int i = 1; i < nCoordDim; ++i)
    {
        if (fabs(projDir[i]) > fabs(projDir[dirMaxId]))
        {
            dirMaxId = i;
        }
    }
 
    Array<OneD, int> dirUse(3, 0);
    if (nCoordDim == 2)
    {
        // 2D or 2.5D cases
        if (dirMaxId == 0)
        {
            dirUse[0] = 1;
            dirUse[1] = 0;
            dirUse[2] = 2;
        }
        else
        {
            dirUse[0] = 0;
            dirUse[1] = 1;
            dirUse[2] = 2;
        }
    }
    else
    {
        // 3D cases
        if (dirMaxId == 0)
        {
            dirUse[0] = 1;
            dirUse[1] = 2;
            dirUse[2] = 0;
        }
        else if (dirMaxId == 1)
        {
            dirUse[0] = 2;
            dirUse[1] = 0;
            dirUse[2] = 1;
        }
        else
        {
            dirUse[0] = 0;
            dirUse[1] = 1;
            dirUse[2] = 2;
        }
    }
 
    // Check if the projected point is in the projected elmt
    // If yes, then compute the locCoord and check if desired point is found
    if (nCoordDim == 2)
    {
        if (isInProjectedArea2D(projGloCoord, projPts, 1.0e-12))
        {
            bool isConverge, isDesired;
 
            isConverge = BisectionForLocCoordOnBndElmt(
                bndGeom, projGloCoord, projPts, dirUse, locCoord, iterTol);
 
            Array<OneD, NekDouble> tmp(2, 0.0);
            tmp[0]   = bndXmap->PhysEvaluate(locCoord, pts[0]) - gloCoord[0];
            tmp[1]   = bndXmap->PhysEvaluate(locCoord, pts[1]) - gloCoord[1];
            projDist = Vmath::Dot(2, tmp, 1, projDir, 1); // can be negative
 
            isDesired = (projDist > 0.0) && (projDist < maxDist);
 
            return isConverge && isDesired;
        }
        else
        {
            return false;
        }
    }
    else
    {
        if (isInProjectedArea3D(projGloCoord, projPts, projDir, 1.0e-12,
                                1.0e-6))
        {
            NekDouble dist;
            bool isConverge, isDesired;
 
            isConverge =
                NewtonIterForLocCoordOnBndElmt(bndGeom, projGloCoord, projPts,
                                               dirUse, locCoord, dist, iterTol);
 
            if (dist > iterTol)
            {
                std::ostringstream ss;
                ss << "Collapsed locCoord out of range.\n"
                   << "Newton iteration gives the distance: " << dist;
                WARNINGL1(false, ss.str());
            }
 
            Array<OneD, NekDouble> tmp(3, 0.0);
            tmp[0]   = bndXmap->PhysEvaluate(locCoord, pts[0]) - gloCoord[0];
            tmp[1]   = bndXmap->PhysEvaluate(locCoord, pts[1]) - gloCoord[1];
            tmp[2]   = bndXmap->PhysEvaluate(locCoord, pts[2]) - gloCoord[2];
            projDist = Vmath::Dot(3, tmp, 1, projDir, 1); // can be negative
 
            isDesired = (projDist > 0.0) && (projDist < maxDist);
 
            return isConverge && isDesired;
        }
        else
        {
            return false;
        }
    }
}
 
/**
 * @brief Get the normals for a given locCoord
 * @param bndGeom      Pointer to the geometry of the boundary element.
 * @param locCoord     Iteration results for local coordinates (if inside).
 * @param normals      Wall normal as the result
 */
void ProcessWallNormalData::GetNormals(
    SpatialDomains::GeometrySharedPtr bndGeom,
    const Array<OneD, const NekDouble> &locCoord,
    Array<OneD, NekDouble> &normals)
{
    const int nCoordDim = m_f->m_exp[0]->GetCoordim(0); // =2 for 2.5D cases
    StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
    Array<OneD, Array<OneD, NekDouble>> pts(m_spacedim);
    Array<OneD, Array<OneD, const NekDouble>> bndCoeffs(m_spacedim);
    int npts = bndXmap->GetTotPoints();
 
    // Get pts in the element
    for (int i = 0; i < nCoordDim; ++i)
    {
        pts[i]       = Array<OneD, NekDouble>(npts);
        bndCoeffs[i] = bndGeom->GetCoeffs(i); // 0/1/2 for x/y/z
        bndXmap->BwdTrans(bndCoeffs[i], pts[i]);
    }
 
    // Get the outward-pointing normals according to the given locCoord
    if (nCoordDim == 2)
    {
        Array<OneD, NekDouble> DxD1(pts[0].size());
        Array<OneD, NekDouble> DyD1(pts[0].size());
 
        bndXmap->PhysDeriv(pts[0], DxD1);
        bndXmap->PhysDeriv(pts[1], DyD1);
 
        NekDouble dxd1, dyd1, norm;
        dxd1 = bndXmap->PhysEvaluate(locCoord, DxD1);
        dyd1 = bndXmap->PhysEvaluate(locCoord, DyD1);
        norm = sqrt(dxd1 * dxd1 + dyd1 * dyd1);
 
        normals[0] = dyd1 / norm;
        normals[1] = -dxd1 / norm;
        normals[2] = 0.0;
    }
    else
    {
        Array<OneD, NekDouble> DxD1(pts[0].size());
        Array<OneD, NekDouble> DxD2(pts[0].size());
        Array<OneD, NekDouble> DyD1(pts[0].size());
        Array<OneD, NekDouble> DyD2(pts[0].size());
        Array<OneD, NekDouble> DzD1(pts[0].size());
        Array<OneD, NekDouble> DzD2(pts[0].size());
 
        bndXmap->PhysDeriv(pts[0], DxD1, DxD2);
        bndXmap->PhysDeriv(pts[1], DyD1, DyD2);
        bndXmap->PhysDeriv(pts[2], DzD1, DzD2);
 
        NekDouble dxd1, dyd1, dzd1, dxd2, dyd2, dzd2;
        dxd1 = bndXmap->PhysEvaluate(locCoord, DxD1);
        dxd2 = bndXmap->PhysEvaluate(locCoord, DxD2);
        dyd1 = bndXmap->PhysEvaluate(locCoord, DyD1);
        dyd2 = bndXmap->PhysEvaluate(locCoord, DyD2);
        dzd1 = bndXmap->PhysEvaluate(locCoord, DzD1);
        dzd2 = bndXmap->PhysEvaluate(locCoord, DzD2);
 
        NekDouble n1, n2, n3, norm;
        n1   = dyd1 * dzd2 - dyd2 * dzd1;
        n2   = dzd1 * dxd2 - dzd2 * dxd1;
        n3   = dxd1 * dyd2 - dxd2 * dyd1;
        norm = sqrt(n1 * n1 + n2 * n2 + n3 * n3);
 
        normals[0] = n1 / norm;
        normals[1] = n2 / norm;
        normals[2] = n3 / norm;
    }
}
 
} // namespace Nektar::FieldUtils