3D geometry information More...

#include <Geometry3D.h>

Inheritance diagram for Nektar::SpatialDomains::Geometry3D:

Public Member Functions
	Geometry3D ()

	Geometry3D (const int coordim)

	~Geometry3D () override

Public Member Functions inherited from Nektar::SpatialDomains::Geometry
	Geometry ()
	Default constructor. More...

	Geometry (int coordim)
	Constructor when supplied a coordinate dimension. More...

virtual	~Geometry ()
	Default destructor. More...

int	GetCoordim () const
	Return the coordinate dimension of this object (i.e. the dimension of the space in which this object is embedded). More...

void	SetCoordim (int coordim)
	Sets the coordinate dimension of this object (i.e. the dimension of the space in which this object is embedded). More...

GeomFactorsSharedPtr	GetGeomFactors ()
	Get the geometric factors for this object, generating them if required. More...

GeomFactorsSharedPtr	GetRefGeomFactors (const Array< OneD, const LibUtilities::BasisSharedPtr > &tbasis)

GeomFactorsSharedPtr	GetMetricInfo ()
	Get the geometric factors for this object. More...

LibUtilities::ShapeType	GetShapeType (void)
	Get the geometric shape type of this object. More...

int	GetGlobalID (void) const
	Get the ID of this object. More...

void	SetGlobalID (int globalid)
	Set the ID of this object. More...

int	GetVid (int i) const
	Get the ID of vertex `i` of this object. More...

int	GetEid (int i) const
	Get the ID of edge `i` of this object. More...

int	GetFid (int i) const
	Get the ID of face `i` of this object. More...

int	GetTid (int i) const
	Get the ID of trace `i` of this object. More...

PointGeomSharedPtr	GetVertex (int i) const
	Returns vertex `i` of this object. More...

Geometry1DSharedPtr	GetEdge (int i) const
	Returns edge `i` of this object. More...

Geometry2DSharedPtr	GetFace (int i) const
	Returns face `i` of this object. More...

StdRegions::Orientation	GetEorient (const int i) const
	Returns the orientation of edge `i` with respect to the ordering of edges in the standard element. More...

StdRegions::Orientation	GetForient (const int i) const
	Returns the orientation of face `i` with respect to the ordering of faces in the standard element. More...

int	GetNumVerts () const
	Get the number of vertices of this object. More...

int	GetNumEdges () const
	Get the number of edges of this object. More...

int	GetNumFaces () const
	Get the number of faces of this object. More...

int	GetShapeDim () const
	Get the object's shape dimension. More...

StdRegions::StdExpansionSharedPtr	GetXmap () const
	Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standard element to physical element. More...

const Array< OneD, const NekDouble > &	GetCoeffs (const int i) const
	Return the coefficients of the transformation Geometry::m_xmap in coordinate direction `i`. More...

void	FillGeom ()
	Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements. More...

std::array< NekDouble, 6 >	GetBoundingBox ()
	Generates the bounding box for the element. More...

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)\). More...

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)\). More...

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)\). More...

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. More...

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`. More...

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. More...

bool	MinMaxCheck (const Array< OneD, const NekDouble > &gloCoord)
	Check if given global coord is within the BoundingBox of the element. More...

bool	ClampLocCoords (Array< OneD, NekDouble > &locCoord, NekDouble tol=std::numeric_limits< NekDouble >::epsilon())
	Clamp local coords to be within standard regions [-1, 1]^dim. More...

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. More...

int	GetVertexFaceMap (int i, int j) const
	Returns the standard element face IDs that are connected to a given vertex. More...

int	GetEdgeFaceMap (int i, int j) const
	Returns the standard element edge IDs that are connected to a given face. More...

int	GetEdgeNormalToFaceVert (int i, int j) const
	Returns the standard lement edge IDs that are normal to a given face vertex. More...

int	GetDir (const int i, const int j=0) const
	Returns the element coordinate direction corresponding to a given face coordinate direction. More...

void	Reset (CurveMap &curvedEdges, CurveMap &curvedFaces)
	Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated GeomFactors. More...

void	ResetNonRecursive (CurveMap &curvedEdges, CurveMap &curvedFaces)
	Reset this geometry object non-recursively: unset the current state, zero Geometry::m_coeffs and remove allocated GeomFactors. More...

void	Setup ()

void	GenGeomFactors ()
	Handles generation of geometry factors. More...

Static Public Attributes
static const int	kDim = 3

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. More...

void	NewtonIterationForLocCoord (const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &ptsx, const Array< OneD, const NekDouble > &ptsy, const Array< OneD, const NekDouble > &ptsz, Array< OneD, NekDouble > &Lcoords, NekDouble &dist)

void	NewtonIterationForLocCoord (const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)

void	v_FillGeom () override
	Put all quadrature information into face/edge structure and backward transform. More...

NekDouble	v_GetCoord (const int i, const Array< OneD, const NekDouble > &Lcoord) override
	Given local collapsed coordinate Lcoord return the value of physical coordinate in direction i. More...

void	v_CalculateInverseIsoParam () override

int	v_AllLeftCheck (const Array< OneD, const NekDouble > &gloCoord) override

int	v_GetShapeDim () const override
	Get the object's shape dimension. More...

int	v_GetNumVerts () const override
	Get the number of vertices of this object. More...

int	v_GetNumEdges () const override
	Get the number of edges of this object. More...

int	v_GetNumFaces () const override
	Get the number of faces of this object. More...

PointGeomSharedPtr	v_GetVertex (int i) const override

Geometry1DSharedPtr	v_GetEdge (int i) const override
	Returns edge `i` of this object. More...

Geometry2DSharedPtr	v_GetFace (int i) const override
	Returns face `i` of this object. More...

StdRegions::Orientation	v_GetEorient (const int i) const override
	Returns the orientation of edge `i` with respect to the ordering of edges in the standard element. More...

StdRegions::Orientation	v_GetForient (const int i) const override
	Returns the orientation of face `i` with respect to the ordering of faces in the standard element. More...

Protected Member Functions inherited from Nektar::SpatialDomains::Geometry
virtual PointGeomSharedPtr	v_GetVertex (int i) const =0

virtual Geometry1DSharedPtr	v_GetEdge (int i) const
	Returns edge `i` of this object. More...

virtual Geometry2DSharedPtr	v_GetFace (int i) const
	Returns face `i` of this object. More...

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. More...

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. More...

virtual int	v_GetNumVerts () const
	Get the number of vertices of this object. More...

virtual int	v_GetNumEdges () const
	Get the number of edges of this object. More...

virtual int	v_GetNumFaces () const
	Get the number of faces of this object. More...

virtual int	v_GetShapeDim () const
	Get the object's shape dimension. More...

virtual StdRegions::StdExpansionSharedPtr	v_GetXmap () const
	Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standard element to physical element. More...

virtual void	v_FillGeom ()
	Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements. More...

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)\). More...

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`. More...

virtual NekDouble	v_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. More...

virtual NekDouble	v_FindDistance (const Array< OneD, const NekDouble > &xs, Array< OneD, NekDouble > &xi)

virtual int	v_GetVertexEdgeMap (int i, int j) const
	Returns the standard element edge IDs that are connected to a given vertex. More...

virtual int	v_GetVertexFaceMap (int i, int j) const
	Returns the standard element face IDs that are connected to a given vertex. More...

virtual int	v_GetEdgeFaceMap (int i, int j) const
	Returns the standard element edge IDs that are connected to a given face. More...

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. More...

virtual int	v_GetDir (const int faceidx, const int facedir) const
	Returns the element coordinate direction corresponding to a given face coordinate direction. More...

virtual void	v_Reset (CurveMap &curvedEdges, CurveMap &curvedFaces)
	Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated GeomFactors. More...

virtual void	v_Setup ()

virtual void	v_GenGeomFactors ()=0

void	SetUpCoeffs (const int nCoeffs)
	Initialise the Geometry::m_coeffs array. More...

virtual void	v_CalculateInverseIsoParam ()

Protected Attributes
PointGeomVector	m_verts

SegGeomVector	m_edges

Geometry2DVector	m_faces

std::vector< StdRegions::Orientation >	m_eorient

std::vector< StdRegions::Orientation >	m_forient

int	m_eid

bool	m_ownverts

Protected Attributes inherited from Nektar::SpatialDomains::Geometry
int	m_coordim
	Coordinate dimension of this geometry object. More...

GeomFactorsSharedPtr	m_geomFactors
	Geometric factors. More...

GeomState	m_geomFactorsState
	State of the geometric factors. More...

StdRegions::StdExpansionSharedPtr	m_xmap
	\(\chi\) mapping containing isoparametric transformation. More...

GeomState	m_state
	Enumeration to dictate whether coefficients are filled. More...

bool	m_setupState
	Wether or not the setup routines have been run. More...

GeomType	m_geomType
	Type of geometry. More...

LibUtilities::ShapeType	m_shapeType
	Type of shape. More...

int	m_globalID
	Global ID. More...

Array< OneD, Array< OneD, NekDouble > >	m_coeffs
	Array containing expansion coefficients of `m_xmap`. More...

Array< OneD, NekDouble >	m_boundingBox
	Array containing bounding box. More...

Array< OneD, Array< OneD, NekDouble > >	m_isoParameter

Array< OneD, Array< OneD, NekDouble > >	m_invIsoParam

int	m_straightEdge

Additional Inherited Members
Static Protected Member Functions inherited from Nektar::SpatialDomains::Geometry
static GeomFactorsSharedPtr	ValidateRegGeomFactor (GeomFactorsSharedPtr geomFactor)
	Check to see if a geometric factor has already been created that contains the same regular information. More...

Static Protected Attributes inherited from Nektar::SpatialDomains::Geometry
static GeomFactorsVector	m_regGeomFactorsManager

Detailed Description

3D geometry information

Definition at line 64 of file Geometry3D.h.

Constructor & Destructor Documentation

◆ Geometry3D() [1/2]

Nektar::SpatialDomains::Geometry3D::Geometry3D ( )

Definition at line 49 of file Geometry3D.cpp.

50{

51}

◆ Geometry3D() [2/2]

Nektar::SpatialDomains::Geometry3D::Geometry3D ( const int coordim )

Definition at line 53 of file Geometry3D.cpp.

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

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

◆ ~Geometry3D()

Nektar::SpatialDomains::Geometry3D::~Geometry3D ( )

override

Definition at line 59 of file Geometry3D.cpp.

60{

61}

Member Function Documentation

◆ NewtonIterationForLocCoord() [1/2]

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

protected

Definition at line 70 of file Geometry3D.cpp.

{
    // maximum iterations for convergence
    const int MaxIterations = 51;
    // |x-xp|^2 < EPSILON  error tolerance
    const NekDouble Tol = 1.e-8;
    // |r,s|    > LcoordDIV stop the search
    const NekDouble LcoordDiv = 15.0;
 
    Array<OneD, const NekDouble> Jac =
        m_geomFactors->GetJac(m_xmap->GetPointsKeys());
 
    NekDouble ScaledTol =
        Vmath::Vsum(Jac.size(), Jac, 1) / ((NekDouble)Jac.size());
    ScaledTol *= Tol * Tol;
 
    NekDouble xmap, ymap, zmap, res;
    int cnt = 0;
    Array<OneD, NekDouble> var(8, 1.);
    Array<OneD, NekDouble> deriv(9);
    Array<OneD, NekDouble> tmp(3);
    while (cnt++ < MaxIterations)
    {
        var[1] = Lcoords[0];
        var[2] = Lcoords[1];
        var[3] = Lcoords[2];
        var[4] = Lcoords[0] * Lcoords[1];
        var[5] = Lcoords[1] * Lcoords[2];
        var[6] = Lcoords[0] * Lcoords[2];
        var[7] = var[4] * Lcoords[2];
        // calculate the global point corresponding to Lcoords
        xmap =
            Vmath::Dot(m_isoParameter[0].size(), var, 1, m_isoParameter[0], 1);
        ymap =
            Vmath::Dot(m_isoParameter[0].size(), var, 1, m_isoParameter[1], 1);
        zmap =
            Vmath::Dot(m_isoParameter[0].size(), var, 1, m_isoParameter[2], 1);
 
        tmp[0] = coords[0] - xmap;
        tmp[1] = coords[1] - ymap;
        tmp[2] = coords[2] - zmap;
 
        res = tmp[0] * tmp[0] + tmp[1] * tmp[1] + tmp[2] * tmp[2];
        if (res < ScaledTol)
        {
            break;
        }
 
        // Interpolate derivative metric at Lcoords (ddx1, ddx2, ddx3)
        deriv[0] = m_isoParameter[0][1] + m_isoParameter[0][4] * Lcoords[1];
        deriv[1] = m_isoParameter[0][2] + m_isoParameter[0][4] * Lcoords[0];
        deriv[2] = m_isoParameter[0][3];
        deriv[3] = m_isoParameter[1][1] + m_isoParameter[1][4] * Lcoords[1];
        deriv[4] = m_isoParameter[1][2] + m_isoParameter[1][4] * Lcoords[0];
        deriv[5] = m_isoParameter[1][3];
        deriv[6] = m_isoParameter[2][1] + m_isoParameter[2][4] * Lcoords[1];
        deriv[7] = m_isoParameter[2][2] + m_isoParameter[2][4] * Lcoords[0];
        deriv[8] = m_isoParameter[2][3];
        if (m_isoParameter[0].size() >= 6)
        {
            deriv[1] += m_isoParameter[0][5] * Lcoords[2];
            deriv[2] += m_isoParameter[0][5] * Lcoords[1];
            deriv[4] += m_isoParameter[1][5] * Lcoords[2];
            deriv[5] += m_isoParameter[1][5] * Lcoords[1];
            deriv[7] += m_isoParameter[2][5] * Lcoords[2];
            deriv[8] += m_isoParameter[2][5] * Lcoords[1];
        }
        if (m_isoParameter[0].size() >= 8)
        {
            deriv[0] += m_isoParameter[0][6] * Lcoords[2] +
                        m_isoParameter[0][7] * var[5];
            deriv[1] += m_isoParameter[0][7] * var[6];
            deriv[2] += m_isoParameter[0][6] * Lcoords[0] +
                        m_isoParameter[0][7] * var[4];
            deriv[3] += m_isoParameter[1][6] * Lcoords[2] +
                        m_isoParameter[1][7] * var[5];
            deriv[4] += m_isoParameter[1][7] * var[6];
            deriv[5] += m_isoParameter[1][6] * Lcoords[0] +
                        m_isoParameter[1][7] * var[4];
            deriv[6] += m_isoParameter[2][6] * Lcoords[2] +
                        m_isoParameter[2][7] * var[5];
            deriv[7] += m_isoParameter[2][7] * var[6];
            deriv[8] += m_isoParameter[2][6] * Lcoords[0] +
                        m_isoParameter[2][7] * var[4];
        }
        DNekMatSharedPtr mat =
            MemoryManager<DNekMat>::AllocateSharedPtr(3, 3, 0, eFULL);
        Vmath::Vcopy(9, deriv, 1, mat->GetPtr(), 1);
        mat->Invert();
        Lcoords[0] += Vmath::Dot(3, mat->GetPtr(), 1, tmp, 1);
        Lcoords[1] += Vmath::Dot(3, mat->GetPtr() + 3, 1, tmp, 1);
        Lcoords[2] += Vmath::Dot(3, mat->GetPtr() + 6, 1, tmp, 1);
 
        if (!(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1]) &&
              std::isfinite(Lcoords[2])) ||
            fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv ||
            fabs(Lcoords[2]) > LcoordDiv)
        {
            std::ostringstream ss;
            ss << "Iteration has diverged in NewtonIterationForLocCoord in "
                  "element "
               << GetGlobalID();
            WARNINGL1(false, ss.str());
            return;
        }
    }
 
    if (cnt >= MaxIterations)
    {
        std::ostringstream ss;
 
        ss << "Reached MaxIterations (" << MaxIterations
           << ") in Newton iteration ";
 
        WARNINGL1(cnt < MaxIterations, ss.str());
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Vmath::Dot(), Nektar::eFULL, Nektar::SpatialDomains::Geometry::GetGlobalID(), Nektar::SpatialDomains::Geometry::m_geomFactors, Nektar::SpatialDomains::Geometry::m_isoParameter, Nektar::SpatialDomains::Geometry::m_xmap, Vmath::Vcopy(), Vmath::Vsum(), and WARNINGL1.

◆ NewtonIterationForLocCoord() [2/2]

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

protected

Definition at line 196 of file Geometry3D.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;
 
    Array<OneD, const NekDouble> Jac =
        m_geomFactors->GetJac(m_xmap->GetPointsKeys());
 
    NekDouble ScaledTol =
        Vmath::Vsum(Jac.size(), Jac, 1) / ((NekDouble)Jac.size());
    ScaledTol *= Tol;
 
    NekDouble xmap, ymap, zmap, F1, F2, F3;
 
    NekDouble derx_1, derx_2, derx_3, dery_1, dery_2, dery_3, derz_1, derz_2,
        derz_3, jac;
 
    // save intiial guess for later reference if required.
    NekDouble init0 = Lcoords[0], init1 = Lcoords[1], init2 = Lcoords[2];
 
    Array<OneD, NekDouble> DxD1(ptsx.size());
    Array<OneD, NekDouble> DxD2(ptsx.size());
    Array<OneD, NekDouble> DxD3(ptsx.size());
    Array<OneD, NekDouble> DyD1(ptsx.size());
    Array<OneD, NekDouble> DyD2(ptsx.size());
    Array<OneD, NekDouble> DyD3(ptsx.size());
    Array<OneD, NekDouble> DzD1(ptsx.size());
    Array<OneD, NekDouble> DzD2(ptsx.size());
    Array<OneD, NekDouble> DzD3(ptsx.size());
 
    // Ideally this will be stored in m_geomfactors
    m_xmap->PhysDeriv(ptsx, DxD1, DxD2, DxD3);
    m_xmap->PhysDeriv(ptsy, DyD1, DyD2, DyD3);
    m_xmap->PhysDeriv(ptsz, DzD1, DzD2, DzD3);
 
    int cnt = 0;
    Array<OneD, DNekMatSharedPtr> I(3);
    Array<OneD, NekDouble> eta(3);
 
    F1 = F2 = F3    = 2000; // Starting value of Function
    NekDouble resid = sqrt(F1 * F1 + F2 * F2 + F3 * F3);
    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);
        I[2] = m_xmap->GetBasis(2)->GetI(eta + 2);
 
        // calculate the global point `corresponding to Lcoords
        xmap = m_xmap->PhysEvaluate(I, ptsx);
        ymap = m_xmap->PhysEvaluate(I, ptsy);
        zmap = m_xmap->PhysEvaluate(I, ptsz);
 
        F1 = coords[0] - xmap;
        F2 = coords[1] - ymap;
        F3 = coords[2] - zmap;
 
        if (F1 * F1 + F2 * F2 + F3 * F3 < ScaledTol)
        {
            resid = sqrt(F1 * F1 + F2 * F2 + F3 * F3);
            break;
        }
 
        // Interpolate derivative metric at Lcoords
        derx_1 = m_xmap->PhysEvaluate(I, DxD1);
        derx_2 = m_xmap->PhysEvaluate(I, DxD2);
        derx_3 = m_xmap->PhysEvaluate(I, DxD3);
        dery_1 = m_xmap->PhysEvaluate(I, DyD1);
        dery_2 = m_xmap->PhysEvaluate(I, DyD2);
        dery_3 = m_xmap->PhysEvaluate(I, DyD3);
        derz_1 = m_xmap->PhysEvaluate(I, DzD1);
        derz_2 = m_xmap->PhysEvaluate(I, DzD2);
        derz_3 = m_xmap->PhysEvaluate(I, DzD3);
 
        jac = derx_1 * (dery_2 * derz_3 - dery_3 * derz_2) -
              derx_2 * (dery_1 * derz_3 - dery_3 * derz_1) +
              derx_3 * (dery_1 * derz_2 - dery_2 * derz_1);
 
        // use analytical inverse of derivitives which are also similar to
        // those of metric factors.
        Lcoords[0] =
            Lcoords[0] +
            ((dery_2 * derz_3 - dery_3 * derz_2) * (coords[0] - xmap) -
             (derx_2 * derz_3 - derx_3 * derz_2) * (coords[1] - ymap) +
             (derx_2 * dery_3 - derx_3 * dery_2) * (coords[2] - zmap)) /
                jac;
 
        Lcoords[1] =
            Lcoords[1] -
            ((dery_1 * derz_3 - dery_3 * derz_1) * (coords[0] - xmap) -
             (derx_1 * derz_3 - derx_3 * derz_1) * (coords[1] - ymap) +
             (derx_1 * dery_3 - derx_3 * dery_1) * (coords[2] - zmap)) /
                jac;
 
        Lcoords[2] =
            Lcoords[2] +
            ((dery_1 * derz_2 - dery_2 * derz_1) * (coords[0] - xmap) -
             (derx_1 * derz_2 - derx_2 * derz_1) * (coords[1] - ymap) +
             (derx_1 * dery_2 - derx_2 * dery_1) * (coords[2] - zmap)) /
                jac;
 
        if (!(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1]) &&
              std::isfinite(Lcoords[2])))
        {
            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 ||
            fabs(Lcoords[2]) > 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);
        I[2] = m_xmap->GetBasis(2)->GetI(eta + 2);
        // calculate the global point corresponding to Lcoords
        xmap = m_xmap->PhysEvaluate(I, ptsx);
        ymap = m_xmap->PhysEvaluate(I, ptsy);
        zmap = m_xmap->PhysEvaluate(I, ptsz);
        F1   = coords[0] - xmap;
        F2   = coords[1] - ymap;
        F3   = coords[2] - zmap;
        dist = sqrt(F1 * F1 + F2 * F2 + F3 * F3);
    }
    else
    {
        dist = 0.;
    }
 
    if (cnt >= MaxIterations)
    {
        Array<OneD, NekDouble> collCoords(3);
        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) &&
            (collCoords[2] >= -1.0 && collCoords[2] <= 1.0))
        {
            std::ostringstream ss;
 
            ss << "Reached MaxIterations (" << MaxIterations
               << ") in Newton iteration ";
            ss << "Init value (" << setprecision(4) << init0 << "," << init1
               << "," << init2 << ") ";
            ss << "Fin  value (" << Lcoords[0] << "," << Lcoords[1] << ","
               << Lcoords[2] << ") ";
            ss << "Resid = " << resid << " Tolerance = " << sqrt(ScaledTol);
 
            WARNINGL1(cnt < MaxIterations, ss.str());
        }
    }
}

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

Referenced by v_GetLocCoords().

◆ v_AllLeftCheck()

int Nektar::SpatialDomains::Geometry3D::v_AllLeftCheck ( const Array< OneD, const NekDouble > & gloCoord )

overrideprotectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 189 of file Geometry3D.cpp.

{
    // todo: check only plane surfaces
    return 0;
}

◆ v_CalculateInverseIsoParam()

void Nektar::SpatialDomains::Geometry3D::v_CalculateInverseIsoParam ( )

overrideprotectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 570 of file Geometry3D.cpp.

{
    DNekMatSharedPtr mat =
        MemoryManager<DNekMat>::AllocateSharedPtr(3, 3, 0, eFULL);
    for (int i = 0; i < 3; ++i)
    {
        for (int j = 0; j < 3; ++j)
        {
 
            mat->SetValue(i, j, m_isoParameter[i][j + 1]);
        }
    }
    mat->Invert();
    m_invIsoParam = Array<OneD, Array<OneD, NekDouble>>(3);
    for (int i = 0; i < 3; ++i)
    {
        m_invIsoParam[i] = Array<OneD, NekDouble>(3);
        for (int j = 0; j < 3; ++j)
        {
            m_invIsoParam[i][j] = mat->GetValue(i, j);
        }
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::SpatialDomains::Geometry::m_invIsoParam, and Nektar::SpatialDomains::Geometry::m_isoParameter.

Referenced by Nektar::SpatialDomains::HexGeom::v_GenGeomFactors(), Nektar::SpatialDomains::PrismGeom::v_GenGeomFactors(), Nektar::SpatialDomains::PyrGeom::v_GenGeomFactors(), and Nektar::SpatialDomains::TetGeom::v_GenGeomFactors().

◆ v_FillGeom()

void Nektar::SpatialDomains::Geometry3D::v_FillGeom ( )

overrideprotectedvirtual

Put all quadrature information into face/edge structure and backward transform.

Note verts, edges, and faces are listed according to anticlockwise convention but points in _coeffs have to be in array format from left to right.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 447 of file Geometry3D.cpp.

{
    if (m_state == ePtsFilled)
    {
        return;
    }
 
    int i, j, k;
 
    for (i = 0; i < m_forient.size(); i++)
    {
        m_faces[i]->FillGeom();
 
        int nFaceCoeffs = m_faces[i]->GetXmap()->GetNcoeffs();
 
        Array<OneD, unsigned int> mapArray(nFaceCoeffs);
        Array<OneD, int> signArray(nFaceCoeffs);
 
        if (m_forient[i] < 9)
        {
            m_xmap->GetTraceToElementMap(
                i, mapArray, signArray, m_forient[i],
                m_faces[i]->GetXmap()->GetTraceNcoeffs(0),
                m_faces[i]->GetXmap()->GetTraceNcoeffs(1));
        }
        else
        {
            m_xmap->GetTraceToElementMap(
                i, mapArray, signArray, m_forient[i],
                m_faces[i]->GetXmap()->GetTraceNcoeffs(1),
                m_faces[i]->GetXmap()->GetTraceNcoeffs(0));
        }
 
        for (j = 0; j < m_coordim; j++)
        {
            const Array<OneD, const NekDouble> &coeffs =
                m_faces[i]->GetCoeffs(j);
 
            for (k = 0; k < nFaceCoeffs; k++)
            {
                NekDouble v              = signArray[k] * coeffs[k];
                m_coeffs[j][mapArray[k]] = v;
            }
        }
    }
 
    m_state = ePtsFilled;
}

References Nektar::SpatialDomains::ePtsFilled, Nektar::SpatialDomains::Geometry::GetXmap(), Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_coordim, m_faces, m_forient, Nektar::SpatialDomains::Geometry::m_state, and Nektar::SpatialDomains::Geometry::m_xmap.

Referenced by Nektar::SpatialDomains::HexGeom::v_GenGeomFactors(), Nektar::SpatialDomains::PrismGeom::v_GenGeomFactors(), Nektar::SpatialDomains::PyrGeom::v_GenGeomFactors(), Nektar::SpatialDomains::TetGeom::v_GenGeomFactors(), and v_GetLocCoords().

◆ v_GetCoord()

NekDouble Nektar::SpatialDomains::Geometry3D::v_GetCoord	(	const int	i,
		const Array< OneD, const NekDouble > &	Lcoord
	)

overrideprotectedvirtual

Given local collapsed coordinate Lcoord return the value of physical coordinate in direction i.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 500 of file Geometry3D.cpp.

{
    ASSERTL1(m_state == ePtsFilled, "Geometry is not in physical space");
 
    Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
    m_xmap->BwdTrans(m_coeffs[i], tmp);
 
    return m_xmap->PhysEvaluate(Lcoord, tmp);
}

References ASSERTL1, Nektar::SpatialDomains::ePtsFilled, Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_state, and Nektar::SpatialDomains::Geometry::m_xmap.

◆ v_GetEdge()

Geometry1DSharedPtr Nektar::SpatialDomains::Geometry3D::v_GetEdge ( int i ) const

overrideprotectedvirtual

Returns edge i of this object.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 540 of file Geometry3D.cpp.

{
    ASSERTL2(i >= 0 && i <= m_edges.size() - 1,
             "Edge ID must be between 0 and " +
                 boost::lexical_cast<string>(m_edges.size() - 1));
    return m_edges[i];
}

References ASSERTL2, and m_edges.

◆ v_GetEorient()

StdRegions::Orientation Nektar::SpatialDomains::Geometry3D::v_GetEorient ( const int i ) const

inlineoverrideprotectedvirtual

Returns the orientation of edge i with respect to the ordering of edges in the standard element.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 554 of file Geometry3D.cpp.

{
    ASSERTL2(i >= 0 && i <= m_edges.size() - 1,
             "Edge ID must be between 0 and " +
                 boost::lexical_cast<string>(m_edges.size() - 1));
    return m_eorient[i];
}

References ASSERTL2, m_edges, and m_eorient.

◆ v_GetFace()

Geometry2DSharedPtr Nektar::SpatialDomains::Geometry3D::v_GetFace ( int i ) const

overrideprotectedvirtual

Returns face i of this object.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 548 of file Geometry3D.cpp.

{
    ASSERTL2((i >= 0) && (i <= 5), "Edge id must be between 0 and 4");
    return m_faces[i];
}

References ASSERTL2, and m_faces.

◆ v_GetForient()

StdRegions::Orientation Nektar::SpatialDomains::Geometry3D::v_GetForient ( const int i ) const

overrideprotectedvirtual

Returns the orientation of face i with respect to the ordering of faces in the standard element.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 562 of file Geometry3D.cpp.

{
    ASSERTL2(i >= 0 && i <= m_faces.size() - 1,
             "Face ID must be between 0 and " +
                 boost::lexical_cast<string>(m_faces.size() - 1));
    return m_forient[i];
}

References ASSERTL2, m_faces, and m_forient.

◆ v_GetLocCoords()

NekDouble Nektar::SpatialDomains::Geometry3D::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 369 of file Geometry3D.cpp.

{
    NekDouble dist = std::numeric_limits<double>::max();
    Array<OneD, NekDouble> tmpcoords(3);
    tmpcoords[0] = coords[0];
    tmpcoords[1] = coords[1];
    tmpcoords[2] = coords[2];
    if (GetMetricInfo()->GetGtype() == eRegular)
    {
        tmpcoords[0] -= m_isoParameter[0][0];
        tmpcoords[1] -= m_isoParameter[1][0];
        tmpcoords[2] -= m_isoParameter[2][0];
        Lcoords[0] = m_invIsoParam[0][0] * tmpcoords[0] +
                     m_invIsoParam[0][1] * tmpcoords[1] +
                     m_invIsoParam[0][2] * tmpcoords[2];
        Lcoords[1] = m_invIsoParam[1][0] * tmpcoords[0] +
                     m_invIsoParam[1][1] * tmpcoords[1] +
                     m_invIsoParam[1][2] * tmpcoords[2];
        Lcoords[2] = m_invIsoParam[2][0] * tmpcoords[0] +
                     m_invIsoParam[2][1] * tmpcoords[1] +
                     m_invIsoParam[2][2] * tmpcoords[2];
    }
    else if (m_straightEdge)
    {
        ClampLocCoords(Lcoords, 0.);
        NewtonIterationForLocCoord(tmpcoords, Lcoords);
    }
    else
    {
        v_FillGeom();
        // Determine nearest point of coords  to values in m_xmap
        int npts = m_xmap->GetTotPoints();
        Array<OneD, NekDouble> ptsx(npts), ptsy(npts), ptsz(npts);
        Array<OneD, NekDouble> tmp1(npts), tmp2(npts);
 
        m_xmap->BwdTrans(m_coeffs[0], ptsx);
        m_xmap->BwdTrans(m_coeffs[1], ptsy);
        m_xmap->BwdTrans(m_coeffs[2], ptsz);
 
        const Array<OneD, const NekDouble> za = m_xmap->GetPoints(0);
        const Array<OneD, const NekDouble> zb = m_xmap->GetPoints(1);
        const Array<OneD, const NekDouble> zc = m_xmap->GetPoints(2);
 
        // guess the first local coords based on nearest point
        Vmath::Sadd(npts, -coords[0], ptsx, 1, tmp1, 1);
        Vmath::Vmul(npts, tmp1, 1, tmp1, 1, tmp1, 1);
        Vmath::Sadd(npts, -coords[1], ptsy, 1, tmp2, 1);
        Vmath::Vvtvp(npts, tmp2, 1, tmp2, 1, tmp1, 1, tmp1, 1);
        Vmath::Sadd(npts, -coords[2], ptsz, 1, tmp2, 1);
        Vmath::Vvtvp(npts, tmp2, 1, tmp2, 1, tmp1, 1, tmp1, 1);
 
        int min_i = Vmath::Imin(npts, tmp1, 1);
 
        // Get Local coordinates
        int qa = za.size(), qb = zb.size();
 
        Array<OneD, NekDouble> eta(3, 0.);
        eta[2] = zc[min_i / (qa * qb)];
        min_i  = min_i % (qa * qb);
        eta[1] = zb[min_i / qa];
        eta[0] = za[min_i % qa];
        m_xmap->LocCollapsedToLocCoord(eta, Lcoords);
 
        // Perform newton iteration to find local coordinates
        NewtonIterationForLocCoord(coords, ptsx, ptsy, ptsz, Lcoords, dist);
    }
    return dist;
}

References Nektar::SpatialDomains::Geometry::ClampLocCoords(), Nektar::SpatialDomains::eRegular, Nektar::SpatialDomains::Geometry::GetMetricInfo(), Vmath::Imin(), Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_invIsoParam, Nektar::SpatialDomains::Geometry::m_isoParameter, Nektar::SpatialDomains::Geometry::m_straightEdge, Nektar::SpatialDomains::Geometry::m_xmap, NewtonIterationForLocCoord(), Vmath::Sadd(), v_FillGeom(), Vmath::Vmul(), and Vmath::Vvtvp().

◆ v_GetNumEdges()

int Nektar::SpatialDomains::Geometry3D::v_GetNumEdges ( ) const

overrideprotectedvirtual

Get the number of edges of this object.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 525 of file Geometry3D.cpp.

{
    return m_edges.size();
}

References m_edges.

◆ v_GetNumFaces()

int Nektar::SpatialDomains::Geometry3D::v_GetNumFaces ( ) const

overrideprotectedvirtual

Get the number of faces of this object.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 530 of file Geometry3D.cpp.

{
    return m_faces.size();
}

References m_faces.

◆ v_GetNumVerts()

int Nektar::SpatialDomains::Geometry3D::v_GetNumVerts ( ) const

overrideprotectedvirtual

Get the number of vertices of this object.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 520 of file Geometry3D.cpp.

{
    return m_verts.size();
}

References m_verts.

◆ v_GetShapeDim()

int Nektar::SpatialDomains::Geometry3D::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 515 of file Geometry3D.cpp.

{
    return 3;
}

◆ v_GetVertex()

PointGeomSharedPtr Nektar::SpatialDomains::Geometry3D::v_GetVertex ( int i ) const

overrideprotectedvirtual

Implements Nektar::SpatialDomains::Geometry.

Definition at line 535 of file Geometry3D.cpp.

{
    return m_verts[i];
}

References m_verts.

Member Data Documentation

◆ kDim

const int Nektar::SpatialDomains::Geometry3D::kDim = 3

static

Definition at line 74 of file Geometry3D.h.

◆ m_edges

SegGeomVector Nektar::SpatialDomains::Geometry3D::m_edges

protected

Definition at line 78 of file Geometry3D.h.

◆ m_eid

int Nektar::SpatialDomains::Geometry3D::m_eid

protected

Definition at line 82 of file Geometry3D.h.

◆ m_eorient

std::vector<StdRegions::Orientation> Nektar::SpatialDomains::Geometry3D::m_eorient

protected

Definition at line 80 of file Geometry3D.h.

Referenced by Nektar::SpatialDomains::HexGeom::HexGeom(), Nektar::SpatialDomains::PrismGeom::PrismGeom(), Nektar::SpatialDomains::PyrGeom::PyrGeom(), Nektar::SpatialDomains::HexGeom::SetUpEdgeOrientation(), Nektar::SpatialDomains::PrismGeom::SetUpEdgeOrientation(), Nektar::SpatialDomains::PyrGeom::SetUpEdgeOrientation(), Nektar::SpatialDomains::TetGeom::SetUpEdgeOrientation(), Nektar::SpatialDomains::TetGeom::TetGeom(), and v_GetEorient().

◆ m_faces

Geometry2DVector Nektar::SpatialDomains::Geometry3D::m_faces

protected

Definition at line 79 of file Geometry3D.h.

◆ m_forient

std::vector<StdRegions::Orientation> Nektar::SpatialDomains::Geometry3D::m_forient

protected

Definition at line 81 of file Geometry3D.h.

◆ m_ownverts

bool Nektar::SpatialDomains::Geometry3D::m_ownverts

protected

Definition at line 83 of file Geometry3D.h.

◆ m_verts

PointGeomVector Nektar::SpatialDomains::Geometry3D::m_verts

protected

Definition at line 77 of file Geometry3D.h.

Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ Geometry3D() [1/2]

◆ Geometry3D() [2/2]

◆ ~Geometry3D()

Member Function Documentation

◆ NewtonIterationForLocCoord() [1/2]

◆ NewtonIterationForLocCoord() [2/2]

◆ v_AllLeftCheck()

◆ v_CalculateInverseIsoParam()

◆ v_FillGeom()

◆ v_GetCoord()

◆ v_GetEdge()

◆ v_GetEorient()

◆ v_GetFace()

◆ v_GetForient()

◆ v_GetLocCoords()

◆ v_GetNumEdges()

◆ v_GetNumFaces()

◆ v_GetNumVerts()

◆ v_GetShapeDim()

◆ v_GetVertex()

Member Data Documentation

◆ kDim

◆ m_edges

◆ m_eid

◆ m_eorient

◆ m_faces

◆ m_forient

◆ m_ownverts

◆ m_verts