Nektar::SpatialDomains::Geometry2D Class Reference

2D geometry information More...

#include <Geometry2D.h>

Inheritance diagram for Nektar::SpatialDomains::Geometry2D:

Public Member Functions
	Geometry2D ()
	Geometry2D (const int coordim, Curve *curve)
	~Geometry2D () override=default
Curve *	GetCurve ()
void	SetCurve (Curve *curvePtr)
Public Member Functions inherited from Nektar::SpatialDomains::Geometry
	Geometry ()
	Default constructor.
	Geometry (int coordim)
	Constructor when supplied a coordinate dimension.
virtual	~Geometry ()=default
int	GetCoordim () const
	Return the coordinate dimension of this object (i.e. the dimension of the space in which this object is embedded).
void	SetCoordim (int coordim)
	Sets the coordinate dimension of this object (i.e. the dimension of the space in which this object is embedded).
GeomFactorsUniquePtr	GenGeomFactors (LibUtilities::PointsKeyVector &keyTgt)
	Used by Expansion to generate associated GeomFactors.
LibUtilities::ShapeType	GetShapeType (void)
	Get the geometric shape type of this object.
int	GetGlobalID (void) const
	Get the ID of this object.
void	SetGlobalID (int globalid)
	Set the ID of this object.
int	GetVid (int i) const
	Returns global id of vertex i of this object.
int	GetEid (int i) const
	Get the ID of edge i of this object.
int	GetFid (int i) const
	Get the ID of face i of this object.
int	GetTid (int i) const
	Get the ID of trace i of this object.
PointGeom *	GetVertex (int i) const
	Returns vertex i of this object.
Geometry1D *	GetEdge (int i) const
	Returns edge i of this object.
Geometry2D *	GetFace (int i) const
	Returns face i of this object.
StdRegions::Orientation	GetEorient (const int i) const
	Returns the orientation of edge i with respect to the ordering of edges in the standard element.
StdRegions::Orientation	GetForient (const int i) const
	Returns the orientation of face i with respect to the ordering of faces in the standard element.
int	GetNumVerts () const
	Get the number of vertices of this object.
int	GetNumEdges () const
	Get the number of edges of this object.
int	GetNumFaces () const
	Get the number of faces of this object.
int	GetShapeDim () const
	Get the object's shape dimension.
StdRegions::StdExpansionSharedPtr	GetXmap () const
	Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standard element to physical element.
const Array< OneD, const NekDouble > &	GetCoeffs (const int i) const
	Return the coefficients of the transformation Geometry::m_xmap in coordinate direction i.
void	FillGeom ()
	Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements.
std::array< NekDouble, 6 >	GetBoundingBox ()
	Generates the bounding box for the element.
void	ClearBoundingBox ()
bool	ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0)
	Determine whether an element contains a particular Cartesian coordinate \((x,y,z)\).
bool	ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol)
	Determine whether an element contains a particular Cartesian coordinate \((x,y,z)\).
bool	ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol, NekDouble &dist)
	Determine whether an element contains a particular Cartesian coordinate \(\vec{x} = (x,y,z)\).
NekDouble	GetLocCoords (const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
	Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this geometry object.
NekDouble	GetCoord (const int i, const Array< OneD, const NekDouble > &Lcoord)
	Given local collapsed coordinate Lcoord, return the value of physical coordinate in direction i.
int	PreliminaryCheck (const Array< OneD, const NekDouble > &gloCoord)
	A fast and robust check if a given global coord is outside of a deformed element. For regular elements, this check is unnecessary.
bool	MinMaxCheck (const Array< OneD, const NekDouble > &gloCoord)
	Check if given global coord is within the BoundingBox of the element.
bool	ClampLocCoords (Array< OneD, NekDouble > &locCoord, NekDouble tol=std::numeric_limits< NekDouble >::epsilon())
	Clamp local coords to be within standard regions [-1, 1]^dim.
NekDouble	FindDistance (const Array< OneD, const NekDouble > &xs, Array< OneD, NekDouble > &xi)
int	GetVertexEdgeMap (int i, int j) const
	Returns the standard element edge IDs that are connected to a given vertex.
int	GetVertexFaceMap (int i, int j) const
	Returns the standard element face IDs that are connected to a given vertex.
int	GetEdgeFaceMap (int i, int j) const
	Returns the standard element edge IDs that are connected to a given face.
int	GetEdgeNormalToFaceVert (int i, int j) const
	Returns the standard lement edge IDs that are normal to a given face vertex.
int	GetDir (const int i, const int j=0) const
	Returns the element coordinate direction corresponding to a given face coordinate direction.
GeomType	CalcGeomType ()
void	Reset (CurveMap &curvedEdges, CurveMap &curvedFaces)
	Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated GeomFactors.
void	ResetNonRecursive (CurveMap &curvedEdges, CurveMap &curvedFaces)
	Reset this geometry object non-recursively: unset the current state, zero Geometry::m_coeffs and remove allocated GeomFactors.
void	Setup ()

Static Public Attributes
static const int	kDim = 2

Protected Member Functions
NekDouble	v_GetLocCoords (const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords) override
	Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this geometry object.
void	NewtonIterationForLocCoord (const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &ptsx, const Array< OneD, const NekDouble > &ptsy, Array< OneD, NekDouble > &Lcoords, NekDouble &dist)
void	SolveStraightEdgeQuad (const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
void	v_CalculateInverseIsoParam () override
int	v_GetShapeDim () const override
	Get the object's shape dimension.
NekDouble	v_FindDistance (const Array< OneD, const NekDouble > &xs, Array< OneD, NekDouble > &xi) override
Protected Member Functions inherited from Nektar::SpatialDomains::Geometry
virtual int	v_GetVid (int i) const
	Get the ID of vertex i of this object.
virtual PointGeom *	v_GetVertex (const int i) const
	Returns vertex i of this object.
virtual Geometry1D *	v_GetEdge (const int i) const
	Returns edge i of this object.
virtual Geometry2D *	v_GetFace (const int i) const
	Returns face i of this object.
virtual StdRegions::Orientation	v_GetEorient (const int i) const
	Returns the orientation of edge i with respect to the ordering of edges in the standard element.
virtual StdRegions::Orientation	v_GetForient (const int i) const
	Returns the orientation of face i with respect to the ordering of faces in the standard element.
virtual int	v_GetNumVerts () const
	Get the number of vertices of this object.
virtual int	v_GetNumEdges () const
	Get the number of edges of this object.
virtual int	v_GetNumFaces () const
	Get the number of faces of this object.
virtual GeomFactorsUniquePtr	v_GenGeomFactors (LibUtilities::PointsKeyVector &keyTgt)
	Used by Expansion to generate associated GeomFactors.
virtual StdRegions::StdExpansionSharedPtr	v_GetXmap () const
	Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standard element to physical element.
virtual void	v_FillGeom ()
	Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements.
virtual bool	v_ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol, NekDouble &dist)
	Determine whether an element contains a particular Cartesian coordinate \(\vec{x} = (x,y,z)\).
virtual int	v_AllLeftCheck (const Array< OneD, const NekDouble > &gloCoord)
virtual NekDouble	v_GetCoord (const int i, const Array< OneD, const NekDouble > &Lcoord)
	Given local collapsed coordinate Lcoord, return the value of physical coordinate in direction i.
virtual int	v_GetVertexEdgeMap (int i, int j) const
	Returns the standard element edge IDs that are connected to a given vertex.
virtual int	v_GetVertexFaceMap (int i, int j) const
	Returns the standard element face IDs that are connected to a given vertex.
virtual int	v_GetEdgeFaceMap (int i, int j) const
	Returns the standard element edge IDs that are connected to a given face.
virtual int	v_GetEdgeNormalToFaceVert (const int i, const int j) const
	Returns the standard lement edge IDs that are normal to a given face vertex.
virtual int	v_GetDir (const int faceidx, const int facedir) const
	Returns the element coordinate direction corresponding to a given face coordinate direction.
virtual GeomType	v_CalcGeomType ()
virtual void	v_Reset (CurveMap &curvedEdges, CurveMap &curvedFaces)
	Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated GeomFactors.
virtual void	v_Setup ()
void	SetUpCoeffs (const int nCoeffs)
	Initialise the Geometry::m_coeffs array.

Protected Attributes
Curve *	m_curve
Array< OneD, int >	m_manifold
Array< OneD, Array< OneD, NekDouble > >	m_edgeNormal
Protected Attributes inherited from Nektar::SpatialDomains::Geometry
int	m_coordim
	Coordinate dimension of this geometry object.
StdRegions::StdExpansionSharedPtr	m_xmap
	\(\chi\) mapping containing isoparametric transformation.
GeomState	m_state
	Enumeration to dictate whether coefficients are filled.
bool	m_setupState
	Wether or not the setup routines have been run.
LibUtilities::ShapeType	m_shapeType
	Type of shape.
int	m_globalID
	Global ID.
std::vector< Array< OneD, NekDouble > >	m_coeffs
	Array containing expansion coefficients of m_xmap.
Array< OneD, NekDouble >	m_boundingBox
	Array containing bounding box.
Array< OneD, Array< OneD, NekDouble > >	m_isoParameter
Array< OneD, Array< OneD, NekDouble > >	m_invIsoParam
int	m_straightEdge

Detailed Description

2D geometry information

Definition at line 49 of file Geometry2D.h.

Constructor & Destructor Documentation

Geometry2D() [1/2]

Nektar::SpatialDomains::Geometry2D::Geometry2D

(

)

Definition at line 46 of file Geometry2D.cpp.

{
}

Geometry2D() [2/2]

Nektar::SpatialDomains::Geometry2D::Geometry2D	(	const int	coordim,
		Curve *	curve
	)

Definition at line 50 of file Geometry2D.cpp.

    : Geometry(coordim), m_curve(curve)
{
    ASSERTL0(m_coordim > 1,
             "Coordinate dimension should be at least 2 for a 2D geometry");
}

References ASSERTL0, and Nektar::SpatialDomains::Geometry::m_coordim.

~Geometry2D()

Nektar::SpatialDomains::Geometry2D::~Geometry2D ( )

overridedefault

Member Function Documentation

GetCurve()

Curve * Nektar::SpatialDomains::Geometry2D::GetCurve ( )

inline

Definition at line 57 of file Geometry2D.h.

    {
        return m_curve;
    }

References m_curve.

Referenced by export_GeomElements(), and Nektar::SpatialDomains::ZoneBase::ZoneBase().

NewtonIterationForLocCoord()

void Nektar::SpatialDomains::Geometry2D::NewtonIterationForLocCoord ( const Array< OneD, const NekDouble > &  coords, const Array< OneD, const NekDouble > &  ptsx, const Array< OneD, const NekDouble > &  ptsy, Array< OneD, NekDouble > &  Lcoords, NekDouble &  dist  )

protected

Definition at line 84 of file Geometry2D.cpp.

{
    // Maximum iterations for convergence
    const int MaxIterations = NekConstants::kNewtonIterations;
    // |x-xp|^2 < EPSILON  error    tolerance
    const NekDouble Tol = 1.e-8;
    // |r,s|    > LcoordDIV stop   the search
    const NekDouble LcoordDiv = 15.0;
 
    LibUtilities::PointsKeyVector ptsKeys = m_xmap->GetPointsKeys();
    Array<OneD, const NekDouble> Jac      = v_GenGeomFactors(ptsKeys)->GetJac();
 
    NekDouble ScaledTol =
        Vmath::Vsum(Jac.size(), Jac, 1) / ((NekDouble)Jac.size());
    ScaledTol *= Tol;
 
    NekDouble xmap, ymap, F1, F2;
    NekDouble derx_1, derx_2, dery_1, dery_2, jac;
 
    // save intiial guess for later reference if required.
    NekDouble init0 = Lcoords[0], init1 = Lcoords[1];
 
    Array<OneD, NekDouble> DxD1(ptsx.size());
    Array<OneD, NekDouble> DxD2(ptsx.size());
    Array<OneD, NekDouble> DyD1(ptsx.size());
    Array<OneD, NekDouble> DyD2(ptsx.size());
 
    // Ideally this will be stored in m_geomfactors
    m_xmap->PhysDeriv(ptsx, DxD1, DxD2);
    m_xmap->PhysDeriv(ptsy, DyD1, DyD2);
 
    int cnt = 0;
    Array<OneD, DNekMatSharedPtr> I(2);
    Array<OneD, NekDouble> eta(2);
 
    F1 = F2         = 2000; // Starting value of Function
    NekDouble resid = sqrt(F1 * F1 + F2 * F2);
    while (cnt++ < MaxIterations)
    {
        //  evaluate lagrange interpolant at Lcoords
        m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
        I[0] = m_xmap->GetBasis(0)->GetI(eta);
        I[1] = m_xmap->GetBasis(1)->GetI(eta + 1);
 
        // calculate the global point `corresponding to Lcoords
        xmap = m_xmap->PhysEvaluate(I, ptsx);
        ymap = m_xmap->PhysEvaluate(I, ptsy);
 
        F1 = coords[0] - xmap;
        F2 = coords[1] - ymap;
 
        if (F1 * F1 + F2 * F2 < ScaledTol)
        {
            resid = sqrt(F1 * F1 + F2 * F2);
            break;
        }
 
        // Interpolate derivative metric at Lcoords
        derx_1 = m_xmap->PhysEvaluate(I, DxD1);
        derx_2 = m_xmap->PhysEvaluate(I, DxD2);
        dery_1 = m_xmap->PhysEvaluate(I, DyD1);
        dery_2 = m_xmap->PhysEvaluate(I, DyD2);
 
        jac = dery_2 * derx_1 - dery_1 * derx_2;
 
        // use analytical inverse of derivitives which are
        // also similar to those of metric factors.
        Lcoords[0] =
            Lcoords[0] +
            (dery_2 * (coords[0] - xmap) - derx_2 * (coords[1] - ymap)) / jac;
 
        Lcoords[1] =
            Lcoords[1] +
            (-dery_1 * (coords[0] - xmap) + derx_1 * (coords[1] - ymap)) / jac;
 
        if (!(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1])))
        {
            dist = 1e16;
            std::ostringstream ss;
            ss << "nan or inf found in NewtonIterationForLocCoord in element "
               << GetGlobalID();
            WARNINGL1(false, ss.str());
            return;
        }
        if (fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv)
        {
            break; // lcoords have diverged so stop iteration
        }
    }
 
    m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
    if (ClampLocCoords(eta, 0.))
    {
        I[0] = m_xmap->GetBasis(0)->GetI(eta);
        I[1] = m_xmap->GetBasis(1)->GetI(eta + 1);
        // calculate the global point corresponding to Lcoords
        xmap = m_xmap->PhysEvaluate(I, ptsx);
        ymap = m_xmap->PhysEvaluate(I, ptsy);
        F1   = coords[0] - xmap;
        F2   = coords[1] - ymap;
        dist = sqrt(F1 * F1 + F2 * F2);
    }
    else
    {
        dist = 0.;
    }
 
    if (cnt >= MaxIterations)
    {
        Array<OneD, NekDouble> collCoords(2);
        m_xmap->LocCoordToLocCollapsed(Lcoords, 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 (" << MaxIterations
               << ") in Newton iteration ";
            ss << "Init value (" << std::setprecision(4) << init0 << ","
               << init1 << ","
               << ") ";
            ss << "Fin  value (" << Lcoords[0] << "," << Lcoords[1] << ","
               << ") ";
            ss << "Resid = " << resid
               << " Tolerance = " << std::sqrt(ScaledTol);
 
            WARNINGL1(cnt < MaxIterations, ss.str());
        }
    }
}

References Nektar::SpatialDomains::Geometry::ClampLocCoords(), Nektar::SpatialDomains::Geometry::GetGlobalID(), Nektar::NekConstants::kNewtonIterations, Nektar::SpatialDomains::Geometry::m_xmap, tinysimd::sqrt(), Nektar::SpatialDomains::Geometry::v_GenGeomFactors(), Vmath::Vsum(), and WARNINGL1.

Referenced by v_GetLocCoords().

SetCurve()

void Nektar::SpatialDomains::Geometry2D::SetCurve ( Curve *  curvePtr )

inline

Definition at line 61 of file Geometry2D.h.

    {
        m_curve = curvePtr;
    }

References m_curve.

SolveStraightEdgeQuad()

void Nektar::SpatialDomains::Geometry2D::SolveStraightEdgeQuad ( const Array< OneD, const NekDouble > &  coords, Array< OneD, NekDouble > &  Lcoords  )

protected

Definition at line 57 of file Geometry2D.cpp.

{
    int i0 = 0, i1 = 1, j1 = 0, j2 = 1;
    if (m_straightEdge & 2)
    {
        i0 = 1;
        i1 = 0;
    }
    if (m_straightEdge & 4)
    {
        j1 = 1;
        j2 = 0;
    }
    NekDouble beta  = m_isoParameter[1][2];
    NekDouble gamma = m_isoParameter[1][3];
    NekDouble tty   = (coords[i1] - gamma * coords[i0] - m_isoParameter[1][0]) *
                    m_isoParameter[1][1];
    NekDouble denom = 1. / (m_isoParameter[0][2] + m_isoParameter[0][3] * tty);
    NekDouble epsilon = -m_isoParameter[0][3] * beta * denom;
    NekDouble h =
        (m_isoParameter[0][0] + m_isoParameter[0][1] * tty - coords[i0]) *
        denom;
    Lcoords[j2] = -h / (0.5 + sqrt(0.25 - epsilon * h));
    Lcoords[j1] = -beta * Lcoords[j2] + tty;
}

References Nektar::LibUtilities::beta, Nektar::SpatialDomains::Geometry::m_isoParameter, Nektar::SpatialDomains::Geometry::m_straightEdge, and tinysimd::sqrt().

Referenced by v_GetLocCoords().

v_CalculateInverseIsoParam()

void Nektar::SpatialDomains::Geometry2D::v_CalculateInverseIsoParam ( )

overrideprotectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 517 of file Geometry2D.cpp.

{
    NekDouble Jac = m_isoParameter[0][1] * m_isoParameter[1][2] -
                    m_isoParameter[1][1] * m_isoParameter[0][2];
    Jac = 1. / Jac;
    // a12, -a02, -a11, a01
    m_invIsoParam       = Array<OneD, Array<OneD, NekDouble>>(2);
    m_invIsoParam[0]    = Array<OneD, NekDouble>(2);
    m_invIsoParam[1]    = Array<OneD, NekDouble>(2);
    m_invIsoParam[0][0] = m_isoParameter[1][2] * Jac;
    m_invIsoParam[0][1] = -m_isoParameter[0][2] * Jac;
    m_invIsoParam[1][0] = -m_isoParameter[1][1] * Jac;
    m_invIsoParam[1][1] = m_isoParameter[0][1] * Jac;
}

References Nektar::SpatialDomains::Geometry::m_invIsoParam, and Nektar::SpatialDomains::Geometry::m_isoParameter.

Referenced by Nektar::SpatialDomains::QuadGeom::v_CalcGeomType(), and Nektar::SpatialDomains::TriGeom::v_CalcGeomType().

v_FindDistance()

NekDouble Nektar::SpatialDomains::Geometry2D::v_FindDistance ( const Array< OneD, const NekDouble > &  xs, Array< OneD, NekDouble > &  xi  )

overrideprotectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 290 of file Geometry2D.cpp.

{
    GeomType Gtype = CalcGeomType();
    if (Gtype == eRegular)
    {
        xiOut = Array<OneD, NekDouble>(2, 0.0);
 
        GetLocCoords(xs, xiOut);
        ClampLocCoords(xiOut);
 
        Array<OneD, NekDouble> gloCoord(3);
        gloCoord[0] = GetCoord(0, xiOut);
        gloCoord[1] = GetCoord(1, xiOut);
        gloCoord[2] = GetCoord(2, xiOut);
 
        return sqrt((xs[0] - gloCoord[0]) * (xs[0] - gloCoord[0]) +
                    (xs[1] - gloCoord[1]) * (xs[1] - gloCoord[1]) +
                    (xs[2] - gloCoord[2]) * (xs[2] - gloCoord[2]));
    }
    // If deformed edge then the inverse mapping is non-linear so need to
    // numerically solve for the local coordinate
    else if (Gtype == eDeformed)
    {
        // Choose starting based on closest quad
        Array<OneD, NekDouble> xi(2, 0.0), eta(2, 0.0);
        m_xmap->LocCollapsedToLocCoord(eta, xi);
 
        // Armijo constants:
        // https://en.wikipedia.org/wiki/Backtracking_line_search
        const NekDouble c1 = 1e-4, c2 = 0.9;
 
        int nq = m_xmap->GetTotPoints();
 
        Array<OneD, NekDouble> x(nq), y(nq), z(nq);
        m_xmap->BwdTrans(m_coeffs[0], x);
        m_xmap->BwdTrans(m_coeffs[1], y);
        m_xmap->BwdTrans(m_coeffs[2], z);
 
        Array<OneD, NekDouble> xderxi1(nq, 0.0), yderxi1(nq, 0.0),
            zderxi1(nq, 0.0), xderxi2(nq, 0.0), yderxi2(nq, 0.0),
            zderxi2(nq, 0.0), xderxi1xi1(nq, 0.0), yderxi1xi1(nq, 0.0),
            zderxi1xi1(nq, 0.0), xderxi1xi2(nq, 0.0), yderxi1xi2(nq, 0.0),
            zderxi1xi2(nq, 0.0), xderxi2xi1(nq, 0.0), yderxi2xi1(nq, 0.0),
            zderxi2xi1(nq, 0.0), xderxi2xi2(nq, 0.0), yderxi2xi2(nq, 0.0),
            zderxi2xi2(nq, 0.0);
 
        // Get first & second derivatives & partial derivatives of x,y,z values
        std::array<NekDouble, 3> xc_derxi, yc_derxi, zc_derxi;
 
        m_xmap->PhysDeriv(x, xderxi1, xderxi2);
        m_xmap->PhysDeriv(y, yderxi1, yderxi2);
        m_xmap->PhysDeriv(z, zderxi1, zderxi2);
 
        m_xmap->PhysDeriv(xderxi1, xderxi1xi1, xderxi1xi2);
        m_xmap->PhysDeriv(yderxi1, yderxi1xi1, yderxi1xi2);
        m_xmap->PhysDeriv(zderxi1, zderxi1xi1, zderxi1xi2);
 
        m_xmap->PhysDeriv(yderxi2, yderxi2xi1, yderxi2xi2);
        m_xmap->PhysDeriv(xderxi2, xderxi2xi1, xderxi2xi2);
        m_xmap->PhysDeriv(zderxi2, zderxi2xi1, zderxi2xi2);
 
        // Minimisation loop (Quasi-newton method)
        NekDouble fx_prev = std::numeric_limits<NekDouble>::max();
        for (int i = 0; i < NekConstants::kNewtonIterations; ++i)
        {
            // Compute the objective function, f(x_k) and its derivatives
            NekDouble xc = m_xmap->PhysEvaluate(xi, x, xc_derxi);
            NekDouble yc = m_xmap->PhysEvaluate(xi, y, yc_derxi);
            NekDouble zc = m_xmap->PhysEvaluate(xi, z, zc_derxi);
 
            NekDouble xc_derxi1xi1 = m_xmap->PhysEvaluate(xi, xderxi1xi1);
            NekDouble yc_derxi1xi1 = m_xmap->PhysEvaluate(xi, yderxi1xi1);
            NekDouble zc_derxi1xi1 = m_xmap->PhysEvaluate(xi, zderxi1xi1);
 
            NekDouble xc_derxi1xi2 = m_xmap->PhysEvaluate(xi, xderxi1xi2);
            NekDouble yc_derxi1xi2 = m_xmap->PhysEvaluate(xi, yderxi1xi2);
            NekDouble zc_derxi1xi2 = m_xmap->PhysEvaluate(xi, zderxi1xi2);
 
            NekDouble xc_derxi2xi2 = m_xmap->PhysEvaluate(xi, xderxi2xi2);
            NekDouble yc_derxi2xi2 = m_xmap->PhysEvaluate(xi, yderxi2xi2);
            NekDouble zc_derxi2xi2 = m_xmap->PhysEvaluate(xi, zderxi2xi2);
 
            // Objective function is the distance to the search point
            NekDouble xdiff = xc - xs[0];
            NekDouble ydiff = yc - xs[1];
            NekDouble zdiff = zc - xs[2];
 
            NekDouble fx = xdiff * xdiff + ydiff * ydiff + zdiff * zdiff;
 
            NekDouble fx_derxi1 = 2.0 * xdiff * xc_derxi[0] +
                                  2.0 * ydiff * yc_derxi[0] +
                                  2.0 * zdiff * zc_derxi[0];
 
            NekDouble fx_derxi2 = 2.0 * xdiff * xc_derxi[1] +
                                  2.0 * ydiff * yc_derxi[1] +
                                  2.0 * zdiff * zc_derxi[1];
 
            NekDouble fx_derxi1xi1 =
                2.0 * xdiff * xc_derxi1xi1 + 2.0 * xc_derxi[0] * xc_derxi[0] +
                2.0 * ydiff * yc_derxi1xi1 + 2.0 * yc_derxi[0] * yc_derxi[0] +
                2.0 * zdiff * zc_derxi1xi1 + 2.0 * zc_derxi[0] * zc_derxi[0];
 
            NekDouble fx_derxi1xi2 =
                2.0 * xdiff * xc_derxi1xi2 + 2.0 * xc_derxi[1] * xc_derxi[0] +
                2.0 * ydiff * yc_derxi1xi2 + 2.0 * yc_derxi[1] * yc_derxi[0] +
                2.0 * zdiff * zc_derxi1xi2 + 2.0 * zc_derxi[1] * zc_derxi[0];
 
            NekDouble fx_derxi2xi2 =
                2.0 * xdiff * xc_derxi2xi2 + 2.0 * xc_derxi[1] * xc_derxi[1] +
                2.0 * ydiff * yc_derxi2xi2 + 2.0 * yc_derxi[1] * yc_derxi[1] +
                2.0 * zdiff * zc_derxi2xi2 + 2.0 * zc_derxi[1] * zc_derxi[1];
 
            // Jacobian
            NekDouble jac[2];
            jac[0] = fx_derxi1;
            jac[1] = fx_derxi2;
 
            // Inverse of 2x2 hessian
            NekDouble hessInv[2][2];
 
            NekDouble det =
                1 / (fx_derxi1xi1 * fx_derxi2xi2 - fx_derxi1xi2 * fx_derxi1xi2);
            hessInv[0][0] = det * fx_derxi2xi2;
            hessInv[0][1] = det * -fx_derxi1xi2;
            hessInv[1][0] = det * -fx_derxi1xi2;
            hessInv[1][1] = det * fx_derxi1xi1;
 
            // Check for convergence
            if (abs(fx - fx_prev) < 1e-12)
            {
                fx_prev = fx;
                break;
            }
            else
            {
                fx_prev = fx;
            }
 
            NekDouble gamma = 1.0;
            bool conv       = false;
 
            // Search direction: Newton's method
            NekDouble pk[2];
            pk[0] = -(hessInv[0][0] * jac[0] + hessInv[1][0] * jac[1]);
            pk[1] = -(hessInv[0][1] * jac[0] + hessInv[1][1] * jac[1]);
 
            // Backtracking line search
            while (gamma > 1e-10)
            {
                Array<OneD, NekDouble> xi_pk(2);
                xi_pk[0] = xi[0] + pk[0] * gamma;
                xi_pk[1] = xi[1] + pk[1] * gamma;
 
                Array<OneD, NekDouble> eta_pk(2, 0.0);
                m_xmap->LocCoordToLocCollapsed(xi_pk, eta_pk);
 
                if (eta_pk[0] <
                        (-1 - std::numeric_limits<NekDouble>::epsilon()) ||
                    eta_pk[0] >
                        (1 + std::numeric_limits<NekDouble>::epsilon()) ||
                    eta_pk[1] <
                        (-1 - std::numeric_limits<NekDouble>::epsilon()) ||
                    eta_pk[1] > (1 + std::numeric_limits<NekDouble>::epsilon()))
                {
                    gamma /= 2.0;
                    continue;
                }
 
                std::array<NekDouble, 3> xc_pk_derxi, yc_pk_derxi, zc_pk_derxi;
 
                NekDouble xc_pk = m_xmap->PhysEvaluate(xi_pk, x, xc_pk_derxi);
                NekDouble yc_pk = m_xmap->PhysEvaluate(xi_pk, y, yc_pk_derxi);
                NekDouble zc_pk = m_xmap->PhysEvaluate(xi_pk, z, zc_pk_derxi);
 
                NekDouble xc_pk_diff = xc_pk - xs[0];
                NekDouble yc_pk_diff = yc_pk - xs[1];
                NekDouble zc_pk_diff = zc_pk - xs[2];
 
                NekDouble fx_pk = xc_pk_diff * xc_pk_diff +
                                  yc_pk_diff * yc_pk_diff +
                                  zc_pk_diff * zc_pk_diff;
 
                NekDouble fx_pk_derxi1 = 2.0 * xc_pk_diff * xc_pk_derxi[0] +
                                         2.0 * yc_pk_diff * yc_pk_derxi[0] +
                                         2.0 * zc_pk_diff * zc_pk_derxi[0];
 
                NekDouble fx_pk_derxi2 = 2.0 * xc_pk_diff * xc_pk_derxi[1] +
                                         2.0 * yc_pk_diff * yc_pk_derxi[1] +
                                         2.0 * zc_pk_diff * zc_pk_derxi[1];
 
                // Check Wolfe conditions using Armijo constants
                // https://en.wikipedia.org/wiki/Wolfe_conditions
                NekDouble tmp  = pk[0] * fx_derxi1 + pk[1] * fx_derxi2;
                NekDouble tmp2 = pk[0] * fx_pk_derxi1 + pk[1] * fx_pk_derxi2;
                if ((fx_pk - (fx + c1 * gamma * tmp)) <
                        std::numeric_limits<NekDouble>::epsilon() &&
                    (-tmp2 - (-c2 * tmp)) <
                        std::numeric_limits<NekDouble>::epsilon())
                {
                    conv = true;
                    break;
                }
 
                gamma /= 2.0;
            }
 
            if (!conv)
            {
                break;
            }
 
            xi[0] += gamma * pk[0];
            xi[1] += gamma * pk[1];
        }
 
        xiOut = xi;
        return sqrt(fx_prev);
    }
    else
    {
        NEKERROR(ErrorUtil::efatal, "Geometry type unknown");
    }
 
    return -1.0;
}

References tinysimd::abs(), Nektar::SpatialDomains::Geometry::CalcGeomType(), Nektar::SpatialDomains::Geometry::ClampLocCoords(), Nektar::SpatialDomains::eDeformed, Nektar::ErrorUtil::efatal, Nektar::SpatialDomains::eRegular, Nektar::SpatialDomains::Geometry::GetCoord(), Nektar::SpatialDomains::Geometry::GetLocCoords(), Nektar::NekConstants::kNewtonIterations, Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_xmap, NEKERROR, and tinysimd::sqrt().

v_GetLocCoords()

NekDouble Nektar::SpatialDomains::Geometry2D::v_GetLocCoords ( const Array< OneD, const NekDouble > &  coords, Array< OneD, NekDouble > &  Lcoords  )

overrideprotectedvirtual

Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this geometry object.

For curvilinear and non-affine elements (i.e. where the Jacobian varies as a function of the standard element coordinates), this is a non-linear optimisation problem that requires the use of a Newton iteration. Note therefore that this can be an expensive operation.

Note that, clearly, the provided Cartesian coordinate lie outside the element. The function therefore returns the minimum distance from some position in the element to . Lcoords will also be constrained to fit within the range \([-1,1]^d\) where \( d \) is the dimension of the element.

Parameters

coords	Input Cartesian global coordinates
Lcoords	Corresponding local coordinates

Returns

Distance between obtained coordinates and provided ones.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 221 of file Geometry2D.cpp.

{
    NekDouble dist = std::numeric_limits<double>::max();
    Array<OneD, NekDouble> tmpcoords(2);
    tmpcoords[0]   = coords[m_manifold[0]];
    tmpcoords[1]   = coords[m_manifold[1]];
    GeomType Gtype = CalcGeomType();
    if (Gtype == eRegular)
    {
        tmpcoords[0] -= m_isoParameter[0][0];
        tmpcoords[1] -= m_isoParameter[1][0];
        Lcoords[0] = m_invIsoParam[0][0] * tmpcoords[0] +
                     m_invIsoParam[0][1] * tmpcoords[1];
        Lcoords[1] = m_invIsoParam[1][0] * tmpcoords[0] +
                     m_invIsoParam[1][1] * tmpcoords[1];
    }
    else if (m_straightEdge)
    {
        SolveStraightEdgeQuad(tmpcoords, Lcoords);
    }
    else if (Gtype == eDeformed)
    {
        v_FillGeom();
        // Determine nearest point of coords  to values in m_xmap
        int npts = m_xmap->GetTotPoints();
        Array<OneD, NekDouble> ptsx(npts), ptsy(npts);
        Array<OneD, NekDouble> tmpx(npts), tmpy(npts);
 
        // Determine 3D manifold orientation
        m_xmap->BwdTrans(m_coeffs[m_manifold[0]], ptsx);
        m_xmap->BwdTrans(m_coeffs[m_manifold[1]], ptsy);
 
        Array<OneD, NekDouble> eta(2, 0.);
        m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
        ClampLocCoords(eta, 0.);
 
        m_xmap->LocCollapsedToLocCoord(eta, Lcoords);
 
        // Perform newton iteration to find local coordinates
        NewtonIterationForLocCoord(tmpcoords, ptsx, ptsy, Lcoords, dist);
    }
    if (m_coordim == 3)
    {
        Array<OneD, NekDouble> eta(2, 0.), xi(2, 0.);
        m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
        ClampLocCoords(eta, 0.);
        m_xmap->LocCollapsedToLocCoord(eta, xi);
        int npts = m_xmap->GetTotPoints();
        Array<OneD, NekDouble> ptsz(npts);
        m_xmap->BwdTrans(m_coeffs[m_manifold[2]], ptsz);
        NekDouble z = m_xmap->PhysEvaluate(xi, ptsz) - coords[m_manifold[2]];
        if (Gtype == eDeformed)
        {
            dist = sqrt(z * z + dist * dist);
        }
        else
        {
            dist = fabs(z);
        }
    }
    return dist;
}

References Nektar::SpatialDomains::Geometry::CalcGeomType(), Nektar::SpatialDomains::Geometry::ClampLocCoords(), Nektar::SpatialDomains::eDeformed, Nektar::SpatialDomains::eRegular, Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_coordim, Nektar::SpatialDomains::Geometry::m_invIsoParam, Nektar::SpatialDomains::Geometry::m_isoParameter, m_manifold, Nektar::SpatialDomains::Geometry::m_straightEdge, Nektar::SpatialDomains::Geometry::m_xmap, NewtonIterationForLocCoord(), SolveStraightEdgeQuad(), tinysimd::sqrt(), and Nektar::SpatialDomains::Geometry::v_FillGeom().

v_GetShapeDim()

int Nektar::SpatialDomains::Geometry2D::v_GetShapeDim ( ) const

overrideprotectedvirtual

Get the object's shape dimension.

For example, a segment is one dimensional and quadrilateral is two dimensional.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 285 of file Geometry2D.cpp.

{
    return 2;
}

Member Data Documentation

kDim

const int Nektar::SpatialDomains::Geometry2D::kDim = 2

static

Definition at line 56 of file Geometry2D.h.

m_curve

Curve* Nektar::SpatialDomains::Geometry2D::m_curve

protected

Definition at line 67 of file Geometry2D.h.

Referenced by GetCurve(), SetCurve(), Nektar::SpatialDomains::QuadGeom::v_FillGeom(), Nektar::SpatialDomains::TriGeom::v_FillGeom(), Nektar::SpatialDomains::QuadGeom::v_Reset(), and Nektar::SpatialDomains::TriGeom::v_Reset().

m_edgeNormal

Array<OneD, Array<OneD, NekDouble> > Nektar::SpatialDomains::Geometry2D::m_edgeNormal

protected

Definition at line 69 of file Geometry2D.h.

Referenced by Nektar::SpatialDomains::QuadGeom::v_AllLeftCheck(), and Nektar::SpatialDomains::TriGeom::v_AllLeftCheck().

m_manifold

Array<OneD, int> Nektar::SpatialDomains::Geometry2D::m_manifold

protected

Definition at line 68 of file Geometry2D.h.

Referenced by Nektar::SpatialDomains::QuadGeom::v_AllLeftCheck(), Nektar::SpatialDomains::TriGeom::v_AllLeftCheck(), Nektar::SpatialDomains::QuadGeom::v_CalcGeomType(), Nektar::SpatialDomains::TriGeom::v_CalcGeomType(), and v_GetLocCoords().