#include <PrismGeom.h>

Inheritance diagram for Nektar::SpatialDomains::PrismGeom:

Collaboration diagram for Nektar::SpatialDomains::PrismGeom:

Public Member Functions
	PrismGeom ()

	PrismGeom (const Geometry2DSharedPtr faces[])

	~PrismGeom ()

Public Member Functions inherited from Nektar::LibUtilities::GraphVertexObject
	GraphVertexObject ()

	GraphVertexObject (const GraphVertexID id)

int	getid ()

void	setid (const GraphVertexID id)

virtual	~GraphVertexObject ()

Public Member Functions inherited from Nektar::SpatialDomains::Geometry3D
	Geometry3D ()

	Geometry3D (const int coordim)

virtual	~Geometry3D ()

int	GetEid (int i) const
	Return the ID of edge i in this element. More...

const Geometry1DSharedPtr	GetEdge (int i) const

Geometry2DSharedPtr	GetFace (int i)
	Return face i in this element. More...

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

Public Member Functions inherited from Nektar::SpatialDomains::Geometry
	Geometry ()

	Geometry (int coordim)

virtual	~Geometry ()

bool	IsElmtConnected (int gvo_id, int locid) const

void	AddElmtConnected (int gvo_id, int locid)

int	NumElmtConnected () const

int	GetCoordim () const

void	SetCoordim (int coordim)

GeomFactorsSharedPtr	GetGeomFactors ()

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

GeomFactorsSharedPtr	GetMetricInfo ()

LibUtilities::ShapeType	GetShapeType (void)

int	GetGlobalID (void)

void	SetGlobalID (int globalid)

int	GetVid (int i) const

int	GetEid (int i) const

int	GetFid (int i) const

int	GetTid (int i) const

int	GetNumVerts () const

PointGeomSharedPtr	GetVertex (int i) const

StdRegions::Orientation	GetEorient (const int i) const

StdRegions::Orientation	GetPorient (const int i) const

StdRegions::Orientation	GetForient (const int i) const

int	GetNumEdges () const

int	GetNumFaces () const

int	GetShapeDim () const

StdRegions::StdExpansionSharedPtr	GetXmap () const

const Array< OneD, const NekDouble > &	GetCoeffs (const int i) const

bool	ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0)

bool	ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol)

bool	ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol, NekDouble &resid)

int	GetVertexEdgeMap (int i, int j) const

int	GetVertexFaceMap (int i, int j) const
	return the id of the $j^{th}$ face attached to the $i^{th}$ vertex More...

int	GetEdgeFaceMap (int i, int j) const

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

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

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

void	SetOwnData ()

const LibUtilities::BasisSharedPtr	GetBasis (const int i)
	Return the j-th basis of the i-th co-ordinate dimension. More...

const LibUtilities::PointsKeyVector	GetPointsKeys ()

void	Reset (CurveMap &curvedEdges, CurveMap &curvedFaces)

Static Public Attributes
static const int	kNverts = 6

static const int	kNedges = 9

static const int	kNqfaces = 3

static const int	kNtfaces = 2

static const int	kNfaces

static const std::string	XMLElementType

Protected Member Functions
virtual void	v_GenGeomFactors ()

virtual NekDouble	v_GetLocCoords (const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)

virtual bool	v_ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, NekDouble tol)
	Determines if a point specified in global coordinates is located within this tetrahedral geometry. More...

virtual bool	v_ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol)

virtual bool	v_ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol, NekDouble &resid)
	Determines if a point specified in global coordinates is located within this tetrahedral geometry. More...

virtual int	v_GetNumVerts () const

virtual int	v_GetNumEdges () const

virtual int	v_GetNumFaces () const

virtual int	v_GetVertexEdgeMap (const int i, const int j) const

virtual int	v_GetVertexFaceMap (const int i, const int j) const

virtual int	v_GetEdgeFaceMap (const int i, const int j) const

virtual int	v_GetDir (const int faceidx, const int facedir) const

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

Protected Member Functions inherited from Nektar::SpatialDomains::Geometry3D
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 &resid)

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

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 int	v_GetShapeDim () const
	Return the dimension of this element. More...

virtual int	v_GetVid (int i) const
	Return the vertex ID of vertex i. More...

virtual PointGeomSharedPtr	v_GetVertex (int i) const
	Return vertex i in this element. More...

virtual const SegGeomSharedPtr	v_GetEdge (int i) const
	Return edge i of this element. More...

virtual StdRegions::Orientation	v_GetEorient (const int i) const
	Return the orientation of edge i in this element. More...

virtual int	v_GetEid (int i) const
	Return the ID of edge i in this element. More...

virtual const Geometry2DSharedPtr	v_GetFace (int i) const
	Return face i in this element. More...

virtual StdRegions::Orientation	v_GetForient (const int i) const
	Return the orientation of face i in this element. More...

virtual int	v_GetFid (int i) const
	Return the ID of face i in this element. More...

virtual int	v_GetEid () const
	Return the ID of this element. More...

virtual int	v_WhichEdge (SegGeomSharedPtr edge)
	Return the local ID of a given edge. More...

virtual int	v_WhichFace (Geometry2DSharedPtr face)
	Return the local ID of a given face. More...

virtual const LibUtilities::BasisSharedPtr	v_GetBasis (const int i)
	Return the j-th basis of the i-th co-ordinate dimension. More...

virtual void	v_AddElmtConnected (int gvo_id, int locid)

virtual bool	v_IsElmtConnected (int gvo_id, int locid) const

virtual int	v_NumElmtConnected () const

virtual void	v_SetOwnData ()

Protected Member Functions inherited from Nektar::SpatialDomains::Geometry
void	GenGeomFactors ()

virtual StdRegions::Orientation	v_GetPorient (const int i) const

virtual StdRegions::StdExpansionSharedPtr	v_GetXmap () const

virtual int	v_GetCoordim () const

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

Private Member Functions
void	SetUpLocalEdges ()

void	SetUpLocalVertices ()

void	SetUpEdgeOrientation ()

void	SetUpFaceOrientation ()

void	SetUpXmap ()
	Set up the m_xmap object by determining the order of each direction from derived faces. More...

Static Private Attributes
static const unsigned int	VertexEdgeConnectivity [6][3]

static const unsigned int	VertexFaceConnectivity [6][3]

static const unsigned int	EdgeFaceConnectivity [9][2]

Additional Inherited Members
Static Protected Member Functions inherited from Nektar::SpatialDomains::Geometry
static GeomFactorsSharedPtr	ValidateRegGeomFactor (GeomFactorsSharedPtr geomFactor)

Protected Attributes inherited from Nektar::LibUtilities::GraphVertexObject
GraphVertexID	m_id

Protected Attributes inherited from Nektar::SpatialDomains::Geometry3D
PointGeomVector	m_verts

SegGeomVector	m_edges

Geometry2DVector	m_faces

std::vector < StdRegions::Orientation >	m_eorient

std::vector < StdRegions::Orientation >	m_forient

std::list< CompToElmt >	m_elmtmap

bool	m_owndata

int	m_eid

bool	m_ownverts

Protected Attributes inherited from Nektar::SpatialDomains::Geometry
int	m_coordim
	coordinate dimension More...

GeomFactorsSharedPtr	m_geomFactors

GeomState	m_geomFactorsState

StdRegions::StdExpansionSharedPtr	m_xmap

GeomState	m_state
	enum identifier to determine if quad points are filled More...

GeomType	m_geomType

LibUtilities::ShapeType	m_shapeType

int	m_globalID

Array< OneD, Array< OneD, NekDouble > >	m_coeffs

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

Detailed Description

Definition at line 48 of file PrismGeom.h.

Constructor & Destructor Documentation

Nektar::SpatialDomains::PrismGeom::PrismGeom ( )

Definition at line 55 of file PrismGeom.cpp.

References Nektar::LibUtilities::ePrism, and Nektar::SpatialDomains::Geometry::m_shapeType.

         {
             m_shapeType = LibUtilities::ePrism;
         }

Nektar::SpatialDomains::PrismGeom::PrismGeom ( const Geometry2DSharedPtr faces[] )

Copy the face shared pointers.

Set up orientation vectors with correct amount of elements.

Set up local objects.

Determine necessary order for standard region.

Definition at line 60 of file PrismGeom.cpp.

References Nektar::LibUtilities::ePrism, kNedges, kNfaces, Nektar::SpatialDomains::Geometry3D::m_eorient, Nektar::SpatialDomains::Geometry3D::m_faces, Nektar::SpatialDomains::Geometry3D::m_forient, Nektar::SpatialDomains::Geometry::m_shapeType, Nektar::SpatialDomains::Geometry::m_xmap, Nektar::SpatialDomains::Geometry::SetUpCoeffs(), SetUpEdgeOrientation(), SetUpFaceOrientation(), SetUpLocalEdges(), SetUpLocalVertices(), and SetUpXmap().

                                                              :
             Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
         {
             m_shapeType = LibUtilities::ePrism;
 
             /// Copy the face shared pointers.
             m_faces.insert(m_faces.begin(), faces, faces+PrismGeom::kNfaces);
 
             /// Set up orientation vectors with correct amount of elements.
             m_eorient.resize(kNedges);
             m_forient.resize(kNfaces);
 
             /// Set up local objects.
             SetUpLocalEdges();
             SetUpLocalVertices();
             SetUpEdgeOrientation();
             SetUpFaceOrientation();
 
             /// Determine necessary order for standard region.
             SetUpXmap();
             SetUpCoeffs(m_xmap->GetNcoeffs());
         }

Nektar::SpatialDomains::PrismGeom::~PrismGeom ( )

Definition at line 83 of file PrismGeom.cpp.

84 {

85 }

Member Function Documentation

void Nektar::SpatialDomains::PrismGeom::SetUpEdgeOrientation ( )

private

Definition at line 582 of file PrismGeom.cpp.

References ASSERTL0, Nektar::StdRegions::eBackwards, Nektar::StdRegions::eForwards, Nektar::SpatialDomains::Geometry::GetVid(), kNedges, Nektar::SpatialDomains::Geometry3D::m_edges, Nektar::SpatialDomains::Geometry3D::m_eorient, and Nektar::SpatialDomains::Geometry3D::m_verts.

Referenced by PrismGeom().

                                             {
 
             // This 2D array holds the local id's of all the vertices
             // for every edge. For every edge, they are ordered to what we
             // define as being Forwards
             const unsigned int edgeVerts[kNedges][2] =
                 { {0,1} ,
                   {1,2} ,
                   {3,2} ,
                   {0,3} ,
                   {0,4} ,
                   {1,4} ,
                   {2,5} ,
                   {3,5} ,
                   {4,5} };
 
 
             int i;
             for(i = 0; i < kNedges; i++)
             {
                 if( m_edges[i]->GetVid(0) == m_verts[ edgeVerts[i][0] ]->GetVid() )
                 {
                     m_eorient[i] = StdRegions::eForwards;
                 }
                 else if( m_edges[i]->GetVid(0) == m_verts[ edgeVerts[i][1] ]->GetVid() )
                 {
                     m_eorient[i] = StdRegions::eBackwards;
                 }
                 else
                 {
                     ASSERTL0(false,"Could not find matching vertex for the edge");
                 }
             }
 
 
         };

void Nektar::SpatialDomains::PrismGeom::SetUpFaceOrientation ( )

private

Definition at line 619 of file PrismGeom.cpp.

References ASSERTL0, Nektar::SpatialDomains::Geometry::GetVertex(), Nektar::SpatialDomains::Geometry::GetVid(), Nektar::NekConstants::kNekZeroTol, kNfaces, kNqfaces, kNtfaces, Nektar::SpatialDomains::QuadGeom::kNverts, Nektar::SpatialDomains::Geometry::m_coordim, Nektar::SpatialDomains::Geometry3D::m_faces, Nektar::SpatialDomains::Geometry3D::m_forient, and Nektar::SpatialDomains::Geometry3D::m_verts.

Referenced by PrismGeom().

                                             {
             int f,i;
 
             // These arrays represent the vector of the A and B
             // coordinate of the local elemental coordinate system
             // where A corresponds with the coordinate direction xi_i
             // with the lowest index i (for that particular face)
             // Coordinate 'B' then corresponds to the other local
             // coordinate (i.e. with the highest index)
             Array<OneD,NekDouble> elementAaxis(m_coordim);
             Array<OneD,NekDouble> elementBaxis(m_coordim);
 
             // These arrays correspond to the local coordinate
             // system of the face itself (i.e. the Geometry2D)
             // faceAaxis correspond to the xi_0 axis
             // faceBaxis correspond to the xi_1 axis
             Array<OneD,NekDouble> faceAaxis(m_coordim);
             Array<OneD,NekDouble> faceBaxis(m_coordim);
 
             // This is the base vertex of the face (i.e. the Geometry2D)
             // This corresponds to thevertex with local ID 0 of the
             // Geometry2D
             unsigned int baseVertex;
 
             // The lenght of the vectors above
             NekDouble elementAaxis_length;
             NekDouble elementBaxis_length;
             NekDouble faceAaxis_length;
             NekDouble faceBaxis_length;
 
             // This 2D array holds the local id's of all the vertices
             // for every face. For every face, they are ordered in such
             // a way that the implementation below allows a unified approach
             // for all faces.
             const unsigned int faceVerts[kNfaces][QuadGeom::kNverts] =
                 { {0,1,2,3} ,
                   {0,1,4,0},  // This is triangle requires only three vertices
                   {1,2,5,4} ,
                   {3,2,5,0},  // This is triangle requires only three vertices
                   {0,3,5,4} ,};
 
             NekDouble dotproduct1 = 0.0;
             NekDouble dotproduct2 = 0.0;
 
             unsigned int orientation;
 
             // Loop over all the faces to set up the orientation
             for(f = 0; f < kNqfaces + kNtfaces; f++)
             {
                 // initialisation
                 elementAaxis_length = 0.0;
                 elementBaxis_length = 0.0;
                 faceAaxis_length = 0.0;
                 faceBaxis_length = 0.0;
 
                 dotproduct1 = 0.0;
                 dotproduct2 = 0.0;
 
                 baseVertex = m_faces[f]->GetVid(0);
 
                 // We are going to construct the vectors representing the A and B axis
                 // of every face. These vectors will be constructed as a vector-representation
                 // of the edges of the face. However, for both coordinate directions, we can
                 // represent the vectors by two different edges. That's why we need to make sure that
                 // we pick the edge to which the baseVertex of the Geometry2D-representation of the face
                 // belongs...
 
                 // Compute the length of edges on a base-face
                 if( f==1  ||  f==3 ) {   // Face is a Triangle
                     for(i = 0; i < m_coordim; i++)
                     {
                         elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
                         elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
                     }
                 }
                 else { // Face is a Quad
                     if( baseVertex == m_verts[ faceVerts[f][0] ]->GetVid() )
                     {
                         for(i = 0; i < m_coordim; i++)
                         {
                             elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
                             elementBaxis[i] = (*m_verts[ faceVerts[f][3] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
                         }
                     }
                     else if( baseVertex == m_verts[ faceVerts[f][1] ]->GetVid() )
                     {
                         for(i = 0; i < m_coordim; i++)
                         {
                             elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
                             elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
                         }
                     }
                     else if( baseVertex == m_verts[ faceVerts[f][2] ]->GetVid() )
                     {
                         for(i = 0; i < m_coordim; i++)
                         {
                             elementAaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][3] ])[i];
                             elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
                         }
                     }
                     else if( baseVertex == m_verts[ faceVerts[f][3] ]->GetVid() )
                     {
                         for(i = 0; i < m_coordim; i++)
                         {
                             elementAaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][3] ])[i];
                             elementBaxis[i] = (*m_verts[ faceVerts[f][3] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
                         }
                     }
                     else
                     {
                         ASSERTL0(false, "Could not find matching vertex for the face");
                     }
                 }
                 // Now, construct the edge-vectors of the local coordinates of
                 // the Geometry2D-representation of the face
                 for(i = 0; i < m_coordim; i++)
                 {
                     int v = m_faces[f]->GetNumVerts()-1;
                     faceAaxis[i] = (*m_faces[f]->GetVertex(1))[i] - (*m_faces[f]->GetVertex(0))[i];
                     faceBaxis[i] = (*m_faces[f]->GetVertex(v))[i] - (*m_faces[f]->GetVertex(0))[i];
 
                     elementAaxis_length += pow(elementAaxis[i],2);
                     elementBaxis_length += pow(elementBaxis[i],2);
                     faceAaxis_length += pow(faceAaxis[i],2);
                     faceBaxis_length += pow(faceBaxis[i],2);
                 }
 
                 elementAaxis_length = sqrt(elementAaxis_length);
                 elementBaxis_length = sqrt(elementBaxis_length);
                 faceAaxis_length = sqrt(faceAaxis_length);
                 faceBaxis_length = sqrt(faceBaxis_length);
 
                 // Calculate the inner product of both the A-axis
                 // (i.e. Elemental A axis and face A axis)
                 for(i = 0 ; i < m_coordim; i++)
                 {
                     dotproduct1 += elementAaxis[i]*faceAaxis[i];
                 }
 
                 orientation = 0;
                 // if the innerproduct is equal to the (absolute value of the ) products of the lengths
                 // of both vectors, then, the coordinate systems will NOT be transposed
                 if( fabs(elementAaxis_length*faceAaxis_length - fabs(dotproduct1)) < NekConstants::kNekZeroTol )
                 {
                     // if the inner product is negative, both A-axis point
                     // in reverse direction
                     if(dotproduct1 < 0.0)
                     {
                         orientation += 2;
                     }
 
                     // calculate the inner product of both B-axis
                     for(i = 0 ; i < m_coordim; i++)
                     {
                         dotproduct2 += elementBaxis[i]*faceBaxis[i];
                     }
 
 //                     // check that both these axis are indeed parallel
 //                     ASSERTL1(fabs(elementBaxis_length*faceBaxis_length - fabs(dotproduct2)) <
 //                              StdRegions::NekConstants::kEvaluateTol,
 //                              "These vectors should be parallel");
 
                     // if the inner product is negative, both B-axis point
                     // in reverse direction
                     if( dotproduct2 < 0.0 )
                     {
                         orientation++;
                     }
                 }
                 // The coordinate systems are transposed
                 else
                 {
                     orientation = 4;
 
                     // Calculate the inner product between the elemental A-axis
                     // and the B-axis of the face (which are now the corresponding axis)
                     dotproduct1 = 0.0;
                     for(i = 0 ; i < m_coordim; i++)
                     {
                         dotproduct1 += elementAaxis[i]*faceBaxis[i];
                     }
 
                     // check that both these axis are indeed parallel
                     if (fabs(elementAaxis_length*faceBaxis_length
                             - fabs(dotproduct1)) > NekConstants::kNekZeroTol)
                     {
                         cout << "Warning: Prism axes not parallel" << endl;
                     }
 
                     // if the result is negative, both axis point in reverse
                     // directions
                     if(dotproduct1 < 0.0)
                     {
                         orientation += 2;
                     }
 
                     // Do the same for the other two corresponding axis
                     dotproduct2 = 0.0;
                     for(i = 0 ; i < m_coordim; i++)
                     {
                         dotproduct2 += elementBaxis[i]*faceAaxis[i];
                     }
 
                     // check that both these axis are indeed parallel
                     if (fabs(elementBaxis_length*faceAaxis_length
                             - fabs(dotproduct2)) > NekConstants::kNekZeroTol)
                     {
                         cout << "Warning: Prism axes not parallel" << endl;
                     }
 
                     if( dotproduct2 < 0.0 )
                     {
                         orientation++;
                     }
                 }
 
                 orientation = orientation + 5;
                 // Fill the m_forient array
                 m_forient[f] = (StdRegions::Orientation) orientation;
             }
         }

void Nektar::SpatialDomains::PrismGeom::SetUpLocalEdges ( )

private

Definition at line 381 of file PrismGeom.cpp.

References ASSERTL0, Nektar::SpatialDomains::Geometry3D::GetEdge(), Nektar::SpatialDomains::Geometry3D::GetEid(), Nektar::SpatialDomains::Geometry::GetFid(), Nektar::SpatialDomains::Geometry3D::m_edges, and Nektar::SpatialDomains::Geometry3D::m_faces.

Referenced by PrismGeom().

                                        {
             // find edge 0
             int i,j;
             unsigned int check;
 
             SegGeomSharedPtr edge;
 
             // First set up the 4 bottom edges
             int f;  //  Connected face index
             for (f = 1; f < 5 ; f++)
             {
                 int nEdges = m_faces[f]->GetNumEdges();
                 check = 0;
                 for (i = 0; i < 4; i++)
                 {
                     for (j = 0; j < nEdges; j++)
                     {
                         if (m_faces[0]->GetEid(i) == m_faces[f]->GetEid(j))
                         {
                             edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[0])->GetEdge(i));
                             m_edges.push_back(edge);
                             check++;
                         }
                     }
                 }
 
                 if (check < 1)
                 {
                     std::ostringstream errstrm;
                     errstrm << "Connected faces do not share an edge. Faces ";
                     errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[f])->GetFid();
                     ASSERTL0(false, errstrm.str());
                 }
                 else if (check > 1)
                 {
                     std::ostringstream errstrm;
                     errstrm << "Connected faces share more than one edge. Faces ";
                     errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[f])->GetFid();
                     ASSERTL0(false, errstrm.str());
                 }
             }
 
             // Then, set up the 4 vertical edges
             check = 0;
             for(i = 0; i < 3; i++) //Set up the vertical edge :face(1) and face(4)
             {
                 for(j = 0; j < 4; j++)
                 {
                     if( (m_faces[1])->GetEid(i) == (m_faces[4])->GetEid(j) )
                     {
                         edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[1])->GetEdge(i));
                         m_edges.push_back(edge);
                         check++;
                     }
                 }
             }
             if( check < 1 )
             {
                 std::ostringstream errstrm;
                 errstrm << "Connected faces do not share an edge. Faces ";
                 errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[4])->GetFid();
                 ASSERTL0(false, errstrm.str());
             }
             else if( check > 1)
             {
                 std::ostringstream errstrm;
                 errstrm << "Connected faces share more than one edge. Faces ";
                 errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[4])->GetFid();
                 ASSERTL0(false, errstrm.str());
             }
             // Set up vertical edges: face(1) through face(4)
             for(f = 1; f < 4 ; f++)
             {
                 check = 0;
                 for(i = 0; i < m_faces[f]->GetNumEdges(); i++)
                 {
                     for(j = 0; j < m_faces[f+1]->GetNumEdges(); j++)
                     {
                         if( (m_faces[f])->GetEid(i) == (m_faces[f+1])->GetEid(j) )
                         {
                             edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[f])->GetEdge(i));
                             m_edges.push_back(edge);
                             check++;
                         }
                     }
                 }
 
                 if( check < 1 )
                 {
                     std::ostringstream errstrm;
                     errstrm << "Connected faces do not share an edge. Faces ";
                     errstrm << (m_faces[f])->GetFid() << ", " << (m_faces[f+1])->GetFid();
                     ASSERTL0(false, errstrm.str());
                 }
                 else if( check > 1)
                 {
                     std::ostringstream errstrm;
                     errstrm << "Connected faces share more than one edge. Faces ";
                     errstrm << (m_faces[f])->GetFid() << ", " << (m_faces[f+1])->GetFid();
                     ASSERTL0(false, errstrm.str());
                 }
             }
 
             // Finally, set up the 1 top edge
             check = 0;
             for(i = 0; i < 4; i++)
             {
                 for(j = 0; j < 4; j++)
                 {
                     if( (m_faces[2])->GetEid(i) == (m_faces[4])->GetEid(j) )
                     {
                         edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[2])->GetEdge(i));
                         m_edges.push_back(edge);
                         check++;
                     }
                 }
             }
 
             if( check < 1 )
             {
                 std::ostringstream errstrm;
                 errstrm << "Connected faces do not share an edge. Faces ";
                 errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[3])->GetFid();
                 ASSERTL0(false, errstrm.str());
             }
             else if( check > 1)
             {
                 std::ostringstream errstrm;
                 errstrm << "Connected faces share more than one edge. Faces ";
                 errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[3])->GetFid();
                 ASSERTL0(false, errstrm.str());
             }
         };

void Nektar::SpatialDomains::PrismGeom::SetUpLocalVertices ( )

private

Definition at line 516 of file PrismGeom.cpp.

References ASSERTL0, Nektar::SpatialDomains::Geometry::GetVertex(), Nektar::SpatialDomains::Geometry::GetVid(), Nektar::SpatialDomains::Geometry3D::m_edges, and Nektar::SpatialDomains::Geometry3D::m_verts.

Referenced by PrismGeom().

                                           {
 
             // Set up the first 2 vertices (i.e. vertex 0,1)
             if( ( m_edges[0]->GetVid(0) == m_edges[1]->GetVid(0) ) ||
                 ( m_edges[0]->GetVid(0) == m_edges[1]->GetVid(1) ) )
             {
                 m_verts.push_back(m_edges[0]->GetVertex(1));
                 m_verts.push_back(m_edges[0]->GetVertex(0));
             }
             else if( ( m_edges[0]->GetVid(1) == m_edges[1]->GetVid(0) ) ||
                      ( m_edges[0]->GetVid(1) == m_edges[1]->GetVid(1) ) )
             {
                 m_verts.push_back(m_edges[0]->GetVertex(0));
                 m_verts.push_back(m_edges[0]->GetVertex(1));
             }
             else
             {
                 std::ostringstream errstrm;
                 errstrm << "Connected edges do not share a vertex. Edges ";
                 errstrm << m_edges[0]->GetEid() << ", " << m_edges[1]->GetEid();
                 ASSERTL0(false, errstrm.str());
             }
 
             // set up the other bottom vertices (i.e. vertex 2,3)
             for(int i = 1; i < 3; i++)
             {
                 if( m_edges[i]->GetVid(0) == m_verts[i]->GetVid() )
                 {
                     m_verts.push_back(m_edges[i]->GetVertex(1));
                 }
                 else if( m_edges[i]->GetVid(1) == m_verts[i]->GetVid() )
                 {
                     m_verts.push_back(m_edges[i]->GetVertex(0));
                 }
                 else
                 {
                     std::ostringstream errstrm;
                     errstrm << "Connected edges do not share a vertex. Edges ";
                     errstrm << m_edges[i]->GetEid() << ", " << m_edges[i-1]->GetEid();
                     ASSERTL0(false, errstrm.str());
                 }
             }
 
             // set up top vertices
             // First, set up vertices 4,5
             if( (m_edges[8]->GetVid(0) == m_edges[4]->GetVid(0)) ||
                 (m_edges[8]->GetVid(0) == m_edges[4]->GetVid(1))  )
             {
                 m_verts.push_back(m_edges[8]->GetVertex(0));
                 m_verts.push_back(m_edges[8]->GetVertex(1));
             }
             else if( (m_edges[8]->GetVid(1) == m_edges[4]->GetVid(0)) ||
                      (m_edges[8]->GetVid(1) == m_edges[4]->GetVid(1))  )
             {
                 m_verts.push_back(m_edges[8]->GetVertex(1));
                 m_verts.push_back(m_edges[8]->GetVertex(0));
             }
             else
             {
                 std::ostringstream errstrm;
                 errstrm << "Connected edges do not share a vertex. Edges ";
                 errstrm << m_edges[8]->GetEid();
                 ASSERTL0(false, errstrm.str());
             }
         };

void Nektar::SpatialDomains::PrismGeom::SetUpXmap ( )

private

Set up the m_xmap object by determining the order of each direction from derived faces.

Definition at line 860 of file PrismGeom.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauMAlpha1Beta0, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::SpatialDomains::Geometry::GetXmap(), Nektar::SpatialDomains::Geometry3D::m_faces, Nektar::SpatialDomains::Geometry3D::m_forient, and Nektar::SpatialDomains::Geometry::m_xmap.

Referenced by PrismGeom(), and v_Reset().

         {
             vector<int> tmp;
 
             int order0, order1;
 
             if (m_forient[0] < 9)
             {
                 tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
                 tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
                 order0 = *max_element(tmp.begin(), tmp.end());
             }
             else
             {
                 tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
                 tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(3));
                 order0 = *max_element(tmp.begin(), tmp.end());
             }
 
             if (m_forient[0] < 9)
             {
                 tmp.clear();
                 tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
                 tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(3));
                 tmp.push_back(m_faces[2]->GetXmap()->GetEdgeNcoeffs(2));
                 order1 = *max_element(tmp.begin(), tmp.end());
             }
             else
             {
                 tmp.clear();
                 tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
                 tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
                 tmp.push_back(m_faces[2]->GetXmap()->GetEdgeNcoeffs(2));
                 order1 = *max_element(tmp.begin(), tmp.end());
             }
 
             tmp.clear();
             tmp.push_back(order0);
             tmp.push_back(order1);
             tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(1));
             tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(2));
             tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(1));
             tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(2));
             int order2 = *max_element(tmp.begin(), tmp.end());
 
             const LibUtilities::BasisKey A(
                 LibUtilities::eModified_A, order0,
                 LibUtilities::PointsKey(
                     order0+1, LibUtilities::eGaussLobattoLegendre));
             const LibUtilities::BasisKey B(
                 LibUtilities::eModified_A, order1,
                 LibUtilities::PointsKey(
                     order1+1, LibUtilities::eGaussLobattoLegendre));
             const LibUtilities::BasisKey C(
                 LibUtilities::eModified_B, order2,
                 LibUtilities::PointsKey(
                     order2, LibUtilities::eGaussRadauMAlpha1Beta0));
 
             m_xmap = MemoryManager<StdRegions::StdPrismExp>::AllocateSharedPtr(
                 A, B, C);
         }

bool Nektar::SpatialDomains::PrismGeom::v_ContainsPoint	(	const Array< OneD, const NekDouble > &	gloCoord,
		NekDouble	tol
	)

protectedvirtual

Determines if a point specified in global coordinates is located within this tetrahedral geometry.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 122 of file PrismGeom.cpp.

References Nektar::SpatialDomains::Geometry::GetCoordim().

Referenced by v_ContainsPoint().

         {
             Array<OneD,NekDouble> locCoord(GetCoordim(),0.0);
             return v_ContainsPoint(gloCoord,locCoord,tol);
         }

bool Nektar::SpatialDomains::PrismGeom::v_ContainsPoint	(	const Array< OneD, const NekDouble > &	gloCoord,
		Array< OneD, NekDouble > &	locCoord,
		NekDouble	tol
	)

protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 129 of file PrismGeom.cpp.

References v_ContainsPoint().

         {
             NekDouble resid;
             return v_ContainsPoint(gloCoord,locCoord,tol,resid);
         }

bool Nektar::SpatialDomains::PrismGeom::v_ContainsPoint	(	const Array< OneD, const NekDouble > &	gloCoord,
		Array< OneD, NekDouble > &	locCoord,
		NekDouble	tol,
		NekDouble &	resid
	)

protectedvirtual

Determines if a point specified in global coordinates is located within this tetrahedral geometry.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 142 of file PrismGeom.cpp.

References ASSERTL1, Nektar::SpatialDomains::eRegular, Nektar::SpatialDomains::Geometry::GetMetricInfo(), Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_xmap, npts, Nektar::SpatialDomains::Geometry3D::v_FillGeom(), v_GetLocCoords(), Vmath::Vmax(), and Vmath::Vmin().

         {
             // Validation checks
             ASSERTL1(gloCoord.num_elements() == 3,
                      "Three dimensional geometry expects three coordinates.");
 
             // find min, max point and check if within twice this
             // distance other false this is advisable since
             // GetLocCoord is expensive for non regular elements.
             if(GetMetricInfo()->GetGtype() !=  eRegular)
             {
                 int i;
                 Array<OneD, NekDouble> mincoord(3), maxcoord(3);
                 NekDouble diff = 0.0;
 
                 v_FillGeom();
 
                 const int npts = m_xmap->GetTotPoints();
                 Array<OneD, NekDouble> pts(npts);
 
                 for(i = 0; i < 3; ++i)
                 {
                     m_xmap->BwdTrans(m_coeffs[i], pts);
 
                     mincoord[i] = Vmath::Vmin(pts.num_elements(),pts,1);
                     maxcoord[i] = Vmath::Vmax(pts.num_elements(),pts,1);
 
                     diff = max(maxcoord[i] - mincoord[i],diff);
                 }
 
                 for(i = 0; i < 3; ++i)
                 {
                     if((gloCoord[i] < mincoord[i] - 0.2*diff)||
                        (gloCoord[i] > maxcoord[i] + 0.2*diff))
                     {
                         return false;
                     }
                 }
             }
 
             // Convert to the local Cartesian coordinates.
             resid = v_GetLocCoords(gloCoord, locCoord);
 
             // Check local coordinate is within std region bounds.
             if (locCoord[0] >= -(1+tol) && locCoord[1] >= -(1+tol) &&
                 locCoord[2] >= -(1+tol) && locCoord[1] <=  (1+tol) &&
                 locCoord[0] + locCoord[2] <= tol)
             {
                 return true;
             }
 
             // If out of range clamp locCoord to be within [-1,1]^3
             // since any larger value will be very oscillatory if
             // called by 'returnNearestElmt' option in
             // ExpList::GetExpIndex
             for(int i = 0; i < 3; ++i)
             {
                 if(locCoord[i] <-(1+tol))
                 {
                     locCoord[i] = -(1+tol);
                 }
 
                 if(locCoord[i] > (1+tol))
                 {
                     locCoord[i] = 1+tol;
                 }
             }
 
             return false;
         }

void Nektar::SpatialDomains::PrismGeom::v_GenGeomFactors ( )

protectedvirtual

Generate the geometry factors for this element.

Reimplemented from Nektar::SpatialDomains::Geometry3D.

Definition at line 217 of file PrismGeom.cpp.

         {
             if (m_geomFactorsState != ePtsFilled)
             {
                 int i,f;
                 GeomType Gtype = eRegular;
 
                 v_FillGeom();
 
                 // check to see if expansions are linear
                 for(i = 0; i < m_coordim; ++i)
                 {
                     if (m_xmap->GetBasisNumModes(0) != 2 ||
                         m_xmap->GetBasisNumModes(1) != 2 ||
                         m_xmap->GetBasisNumModes(2) != 2 )
                     {
                         Gtype = eDeformed;
                     }
                 }
 
                 // check to see if all quadrilateral faces are parallelograms
                 if(Gtype == eRegular)
                 {
                     // Vertex ids of quad faces
                     const unsigned int faceVerts[3][4] =
                         { {0,1,2,3} ,
                           {1,2,5,4} ,
                           {0,3,5,4} };
 
                     for(f = 0; f < 3; f++)
                     {
                         // Ensure each face is a parallelogram? Check this.
                         for (i = 0; i < m_coordim; i++)
                         {
                             if( fabs( (*m_verts[ faceVerts[f][0] ])(i) -
                                       (*m_verts[ faceVerts[f][1] ])(i) +
                                       (*m_verts[ faceVerts[f][2] ])(i) -
                                       (*m_verts[ faceVerts[f][3] ])(i) )
                                     > NekConstants::kNekZeroTol )
                             {
                                 Gtype = eDeformed;
                                 break;
                             }
                         }
 
                         if (Gtype == eDeformed)
                         {
                             break;
                         }
                     }
                 }
 
                 m_geomFactors = MemoryManager<GeomFactors>::AllocateSharedPtr(
                     Gtype, m_coordim, m_xmap, m_coeffs);
                 m_geomFactorsState = ePtsFilled;
             }
         }

int Nektar::SpatialDomains::PrismGeom::v_GetDir	(	const int	faceidx,
		const int	facedir
	)		const

protectedvirtual

Implements Nektar::SpatialDomains::Geometry3D.

Definition at line 102 of file PrismGeom.cpp.

         {
             if (faceidx == 0)
             {
                 return facedir;
             }
             else if (faceidx == 1 || faceidx == 3)
             {
                 return 2 * facedir;
             }
             else
             {
                 return 1 + facedir;
             }
         }

int Nektar::SpatialDomains::PrismGeom::v_GetEdgeFaceMap	(	const int	i,
		const int	j
	)		const

protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 375 of file PrismGeom.cpp.

References EdgeFaceConnectivity.

     {
         return EdgeFaceConnectivity[i][j];
     }

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

protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 276 of file PrismGeom.cpp.

Referenced by v_ContainsPoint().

         {
             NekDouble ptdist = 1e6;
 
             // calculate local coordinate for coord
             if(GetMetricInfo()->GetGtype() == eRegular)
             {
                 // Point inside tetrahedron
                 PointGeom r(m_coordim, 0, coords[0], coords[1], coords[2]);
 
                 // Edges
                 PointGeom er0, e10, e30, e40;
                 er0.Sub(r,*m_verts[0]);
                 e10.Sub(*m_verts[1],*m_verts[0]);
                 e30.Sub(*m_verts[3],*m_verts[0]);
                 e40.Sub(*m_verts[4],*m_verts[0]);
 
                 // Cross products (Normal times area)
                 PointGeom cp1030, cp3040, cp4010;
                 cp1030.Mult(e10,e30);
                 cp3040.Mult(e30,e40);
                 cp4010.Mult(e40,e10);
 
                 // Barycentric coordinates (relative volume)
                 NekDouble V = e40.dot(cp1030); // Prism Volume = {(e40)dot(e10)x(e30)}/2
                 NekDouble beta  = er0.dot(cp3040) / (2.0*V); // volume1 = {(er0)dot(e30)x(e40)}/4
                 NekDouble gamma = er0.dot(cp4010) / (3.0*V); // volume2 = {(er0)dot(e40)x(e10)}/6
                 NekDouble delta = er0.dot(cp1030) / (2.0*V); // volume3 = {(er0)dot(e10)x(e30)}/4
 
                 // Make Prism bigger
                 Lcoords[0] = 2.0*beta  - 1.0;
                 Lcoords[1] = 2.0*gamma - 1.0;
                 Lcoords[2] = 2.0*delta - 1.0;
 
                 // Set ptdist to distance to nearest vertex
                 for(int i = 0; i < 6; ++i)
                 {
                     ptdist = min(ptdist,r.dist(*m_verts[i]));
                 }
             }
             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);
 
                 // distance from coordinate to nearest point for return value.
                 ptdist = sqrt(tmp1[min_i]);
 
                 // Get collapsed coordinate
                 int qa = za.num_elements(), qb = zb.num_elements();
                 Lcoords[2] = zc[min_i/(qa*qb)];
                 min_i = min_i%(qa*qb);
                 Lcoords[1] = zb[min_i/qa];
                 Lcoords[0] = za[min_i%qa];
 
                 // recover cartesian coordinate from collapsed coordinate.
                 Lcoords[0] = (1.0+Lcoords[0])*(1.0-Lcoords[2])/2 - 1.0;
 
                 // Perform newton iteration to find local coordinates
                 NekDouble resid = 0.0;
                 NewtonIterationForLocCoord(coords, ptsx, ptsy, ptsz, Lcoords,resid);
             }
             return ptdist;
         }

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

protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 92 of file PrismGeom.cpp.

         {
             return 9;
         }

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

protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 97 of file PrismGeom.cpp.

         {
             return 5;
         }

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

protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 87 of file PrismGeom.cpp.

         {
             return 6;
         }

int Nektar::SpatialDomains::PrismGeom::v_GetVertexEdgeMap	(	const int	i,
		const int	j
	)		const

protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 365 of file PrismGeom.cpp.

References VertexEdgeConnectivity.

     {
         return VertexEdgeConnectivity[i][j];
     }

int Nektar::SpatialDomains::PrismGeom::v_GetVertexFaceMap	(	const int	i,
		const int	j
	)		const

protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 370 of file PrismGeom.cpp.

References VertexFaceConnectivity.

     {
         return VertexFaceConnectivity[i][j];
     }

void Nektar::SpatialDomains::PrismGeom::v_Reset	(	CurveMap &	curvedEdges,
		CurveMap &	curvedFaces
	)

protectedvirtual

Reset this geometry object: unset the current state and remove allocated GeomFactors.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 841 of file PrismGeom.cpp.

References Nektar::SpatialDomains::Geometry3D::m_faces, Nektar::SpatialDomains::Geometry::m_xmap, Nektar::SpatialDomains::Geometry::SetUpCoeffs(), SetUpXmap(), and Nektar::SpatialDomains::Geometry::v_Reset().

         {
             Geometry::v_Reset(curvedEdges, curvedFaces);
 
             for (int i = 0; i < 5; ++i)
             {
                 m_faces[i]->Reset(curvedEdges, curvedFaces);
             }
 
             SetUpXmap();
             SetUpCoeffs(m_xmap->GetNcoeffs());
         }