Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Attributes | Protected Member Functions | Private Member Functions | Static Private Attributes | List of all members
Nektar::SpatialDomains::PrismGeom Class Reference

#include <PrismGeom.h>

Inheritance diagram for Nektar::SpatialDomains::PrismGeom:
Inheritance graph
[legend]
Collaboration diagram for Nektar::SpatialDomains::PrismGeom:
Collaboration graph
[legend]

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< CompToElmtm_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 57 of file PrismGeom.cpp.

References Nektar::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 62 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().

62  :
63  Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
64  {
66 
67  /// Copy the face shared pointers.
68  m_faces.insert(m_faces.begin(), faces, faces+PrismGeom::kNfaces);
69 
70  /// Set up orientation vectors with correct amount of elements.
71  m_eorient.resize(kNedges);
72  m_forient.resize(kNfaces);
73 
74  /// Set up local objects.
79 
80  /// Determine necessary order for standard region.
81  SetUpXmap();
82  SetUpCoeffs(m_xmap->GetNcoeffs());
83  }
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry3D.h:92
void SetUpXmap()
Set up the m_xmap object by determining the order of each direction from derived faces.
Definition: PrismGeom.cpp:862
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:93
PointGeomSharedPtr GetVertex(int i) const
Definition: Geometry.h:348
const Geometry1DSharedPtr GetEdge(int i) const
LibUtilities::ShapeType m_shapeType
Definition: Geometry.h:177
void SetUpCoeffs(const int nCoeffs)
Initialise the m_coeffs array.
Definition: Geometry.h:484
Nektar::SpatialDomains::PrismGeom::~PrismGeom ( )

Definition at line 85 of file PrismGeom.cpp.

86  {
87  }

Member Function Documentation

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

Definition at line 584 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().

584  {
585 
586  // This 2D array holds the local id's of all the vertices
587  // for every edge. For every edge, they are ordered to what we
588  // define as being Forwards
589  const unsigned int edgeVerts[kNedges][2] =
590  { {0,1} ,
591  {1,2} ,
592  {3,2} ,
593  {0,3} ,
594  {0,4} ,
595  {1,4} ,
596  {2,5} ,
597  {3,5} ,
598  {4,5} };
599 
600 
601  int i;
602  for(i = 0; i < kNedges; i++)
603  {
604  if( m_edges[i]->GetVid(0) == m_verts[ edgeVerts[i][0] ]->GetVid() )
605  {
607  }
608  else if( m_edges[i]->GetVid(0) == m_verts[ edgeVerts[i][1] ]->GetVid() )
609  {
611  }
612  else
613  {
614  ASSERTL0(false,"Could not find matching vertex for the edge");
615  }
616  }
617 
618 
619  };
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry3D.h:92
int GetVid(int i) const
Definition: Geometry.h:319
void Nektar::SpatialDomains::PrismGeom::SetUpFaceOrientation ( )
private

Definition at line 621 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().

621  {
622  int f,i;
623 
624  // These arrays represent the vector of the A and B
625  // coordinate of the local elemental coordinate system
626  // where A corresponds with the coordinate direction xi_i
627  // with the lowest index i (for that particular face)
628  // Coordinate 'B' then corresponds to the other local
629  // coordinate (i.e. with the highest index)
630  Array<OneD,NekDouble> elementAaxis(m_coordim);
631  Array<OneD,NekDouble> elementBaxis(m_coordim);
632 
633  // These arrays correspond to the local coordinate
634  // system of the face itself (i.e. the Geometry2D)
635  // faceAaxis correspond to the xi_0 axis
636  // faceBaxis correspond to the xi_1 axis
637  Array<OneD,NekDouble> faceAaxis(m_coordim);
638  Array<OneD,NekDouble> faceBaxis(m_coordim);
639 
640  // This is the base vertex of the face (i.e. the Geometry2D)
641  // This corresponds to thevertex with local ID 0 of the
642  // Geometry2D
643  unsigned int baseVertex;
644 
645  // The lenght of the vectors above
646  NekDouble elementAaxis_length;
647  NekDouble elementBaxis_length;
648  NekDouble faceAaxis_length;
649  NekDouble faceBaxis_length;
650 
651  // This 2D array holds the local id's of all the vertices
652  // for every face. For every face, they are ordered in such
653  // a way that the implementation below allows a unified approach
654  // for all faces.
655  const unsigned int faceVerts[kNfaces][QuadGeom::kNverts] =
656  { {0,1,2,3} ,
657  {0,1,4,0}, // This is triangle requires only three vertices
658  {1,2,5,4} ,
659  {3,2,5,0}, // This is triangle requires only three vertices
660  {0,3,5,4} ,};
661 
662  NekDouble dotproduct1 = 0.0;
663  NekDouble dotproduct2 = 0.0;
664 
665  unsigned int orientation;
666 
667  // Loop over all the faces to set up the orientation
668  for(f = 0; f < kNqfaces + kNtfaces; f++)
669  {
670  // initialisation
671  elementAaxis_length = 0.0;
672  elementBaxis_length = 0.0;
673  faceAaxis_length = 0.0;
674  faceBaxis_length = 0.0;
675 
676  dotproduct1 = 0.0;
677  dotproduct2 = 0.0;
678 
679  baseVertex = m_faces[f]->GetVid(0);
680 
681  // We are going to construct the vectors representing the A and B axis
682  // of every face. These vectors will be constructed as a vector-representation
683  // of the edges of the face. However, for both coordinate directions, we can
684  // represent the vectors by two different edges. That's why we need to make sure that
685  // we pick the edge to which the baseVertex of the Geometry2D-representation of the face
686  // belongs...
687 
688  // Compute the length of edges on a base-face
689  if( f==1 || f==3 ) { // Face is a Triangle
690  for(i = 0; i < m_coordim; i++)
691  {
692  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
693  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
694  }
695  }
696  else { // Face is a Quad
697  if( baseVertex == m_verts[ faceVerts[f][0] ]->GetVid() )
698  {
699  for(i = 0; i < m_coordim; i++)
700  {
701  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
702  elementBaxis[i] = (*m_verts[ faceVerts[f][3] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
703  }
704  }
705  else if( baseVertex == m_verts[ faceVerts[f][1] ]->GetVid() )
706  {
707  for(i = 0; i < m_coordim; i++)
708  {
709  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
710  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
711  }
712  }
713  else if( baseVertex == m_verts[ faceVerts[f][2] ]->GetVid() )
714  {
715  for(i = 0; i < m_coordim; i++)
716  {
717  elementAaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][3] ])[i];
718  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
719  }
720  }
721  else if( baseVertex == m_verts[ faceVerts[f][3] ]->GetVid() )
722  {
723  for(i = 0; i < m_coordim; i++)
724  {
725  elementAaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][3] ])[i];
726  elementBaxis[i] = (*m_verts[ faceVerts[f][3] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
727  }
728  }
729  else
730  {
731  ASSERTL0(false, "Could not find matching vertex for the face");
732  }
733  }
734  // Now, construct the edge-vectors of the local coordinates of
735  // the Geometry2D-representation of the face
736  for(i = 0; i < m_coordim; i++)
737  {
738  int v = m_faces[f]->GetNumVerts()-1;
739  faceAaxis[i] = (*m_faces[f]->GetVertex(1))[i] - (*m_faces[f]->GetVertex(0))[i];
740  faceBaxis[i] = (*m_faces[f]->GetVertex(v))[i] - (*m_faces[f]->GetVertex(0))[i];
741 
742  elementAaxis_length += pow(elementAaxis[i],2);
743  elementBaxis_length += pow(elementBaxis[i],2);
744  faceAaxis_length += pow(faceAaxis[i],2);
745  faceBaxis_length += pow(faceBaxis[i],2);
746  }
747 
748  elementAaxis_length = sqrt(elementAaxis_length);
749  elementBaxis_length = sqrt(elementBaxis_length);
750  faceAaxis_length = sqrt(faceAaxis_length);
751  faceBaxis_length = sqrt(faceBaxis_length);
752 
753  // Calculate the inner product of both the A-axis
754  // (i.e. Elemental A axis and face A axis)
755  for(i = 0 ; i < m_coordim; i++)
756  {
757  dotproduct1 += elementAaxis[i]*faceAaxis[i];
758  }
759 
760  orientation = 0;
761  // if the innerproduct is equal to the (absolute value of the ) products of the lengths
762  // of both vectors, then, the coordinate systems will NOT be transposed
763  if( fabs(elementAaxis_length*faceAaxis_length - fabs(dotproduct1)) < NekConstants::kNekZeroTol )
764  {
765  // if the inner product is negative, both A-axis point
766  // in reverse direction
767  if(dotproduct1 < 0.0)
768  {
769  orientation += 2;
770  }
771 
772  // calculate the inner product of both B-axis
773  for(i = 0 ; i < m_coordim; i++)
774  {
775  dotproduct2 += elementBaxis[i]*faceBaxis[i];
776  }
777 
778 // // check that both these axis are indeed parallel
779 // ASSERTL1(fabs(elementBaxis_length*faceBaxis_length - fabs(dotproduct2)) <
780 // StdRegions::NekConstants::kEvaluateTol,
781 // "These vectors should be parallel");
782 
783  // if the inner product is negative, both B-axis point
784  // in reverse direction
785  if( dotproduct2 < 0.0 )
786  {
787  orientation++;
788  }
789  }
790  // The coordinate systems are transposed
791  else
792  {
793  orientation = 4;
794 
795  // Calculate the inner product between the elemental A-axis
796  // and the B-axis of the face (which are now the corresponding axis)
797  dotproduct1 = 0.0;
798  for(i = 0 ; i < m_coordim; i++)
799  {
800  dotproduct1 += elementAaxis[i]*faceBaxis[i];
801  }
802 
803  // check that both these axis are indeed parallel
804  if (fabs(elementAaxis_length*faceBaxis_length
805  - fabs(dotproduct1)) > NekConstants::kNekZeroTol)
806  {
807  cout << "Warning: Prism axes not parallel" << endl;
808  }
809 
810  // if the result is negative, both axis point in reverse
811  // directions
812  if(dotproduct1 < 0.0)
813  {
814  orientation += 2;
815  }
816 
817  // Do the same for the other two corresponding axis
818  dotproduct2 = 0.0;
819  for(i = 0 ; i < m_coordim; i++)
820  {
821  dotproduct2 += elementBaxis[i]*faceAaxis[i];
822  }
823 
824  // check that both these axis are indeed parallel
825  if (fabs(elementBaxis_length*faceAaxis_length
826  - fabs(dotproduct2)) > NekConstants::kNekZeroTol)
827  {
828  cout << "Warning: Prism axes not parallel" << endl;
829  }
830 
831  if( dotproduct2 < 0.0 )
832  {
833  orientation++;
834  }
835  }
836 
837  orientation = orientation + 5;
838  // Fill the m_forient array
839  m_forient[f] = (StdRegions::Orientation) orientation;
840  }
841  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:93
static const NekDouble kNekZeroTol
double NekDouble
PointGeomSharedPtr GetVertex(int i) const
Definition: Geometry.h:348
int GetVid(int i) const
Definition: Geometry.h:319
int m_coordim
coordinate dimension
Definition: Geometry.h:169
void Nektar::SpatialDomains::PrismGeom::SetUpLocalEdges ( )
private

Definition at line 383 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().

383  {
384  // find edge 0
385  int i,j;
386  unsigned int check;
387 
388  SegGeomSharedPtr edge;
389 
390  // First set up the 4 bottom edges
391  int f; // Connected face index
392  for (f = 1; f < 5 ; f++)
393  {
394  int nEdges = m_faces[f]->GetNumEdges();
395  check = 0;
396  for (i = 0; i < 4; i++)
397  {
398  for (j = 0; j < nEdges; j++)
399  {
400  if (m_faces[0]->GetEid(i) == m_faces[f]->GetEid(j))
401  {
402  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[0])->GetEdge(i));
403  m_edges.push_back(edge);
404  check++;
405  }
406  }
407  }
408 
409  if (check < 1)
410  {
411  std::ostringstream errstrm;
412  errstrm << "Connected faces do not share an edge. Faces ";
413  errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[f])->GetFid();
414  ASSERTL0(false, errstrm.str());
415  }
416  else if (check > 1)
417  {
418  std::ostringstream errstrm;
419  errstrm << "Connected faces share more than one edge. Faces ";
420  errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[f])->GetFid();
421  ASSERTL0(false, errstrm.str());
422  }
423  }
424 
425  // Then, set up the 4 vertical edges
426  check = 0;
427  for(i = 0; i < 3; i++) //Set up the vertical edge :face(1) and face(4)
428  {
429  for(j = 0; j < 4; j++)
430  {
431  if( (m_faces[1])->GetEid(i) == (m_faces[4])->GetEid(j) )
432  {
433  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[1])->GetEdge(i));
434  m_edges.push_back(edge);
435  check++;
436  }
437  }
438  }
439  if( check < 1 )
440  {
441  std::ostringstream errstrm;
442  errstrm << "Connected faces do not share an edge. Faces ";
443  errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[4])->GetFid();
444  ASSERTL0(false, errstrm.str());
445  }
446  else if( check > 1)
447  {
448  std::ostringstream errstrm;
449  errstrm << "Connected faces share more than one edge. Faces ";
450  errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[4])->GetFid();
451  ASSERTL0(false, errstrm.str());
452  }
453  // Set up vertical edges: face(1) through face(4)
454  for(f = 1; f < 4 ; f++)
455  {
456  check = 0;
457  for(i = 0; i < m_faces[f]->GetNumEdges(); i++)
458  {
459  for(j = 0; j < m_faces[f+1]->GetNumEdges(); j++)
460  {
461  if( (m_faces[f])->GetEid(i) == (m_faces[f+1])->GetEid(j) )
462  {
463  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[f])->GetEdge(i));
464  m_edges.push_back(edge);
465  check++;
466  }
467  }
468  }
469 
470  if( check < 1 )
471  {
472  std::ostringstream errstrm;
473  errstrm << "Connected faces do not share an edge. Faces ";
474  errstrm << (m_faces[f])->GetFid() << ", " << (m_faces[f+1])->GetFid();
475  ASSERTL0(false, errstrm.str());
476  }
477  else if( check > 1)
478  {
479  std::ostringstream errstrm;
480  errstrm << "Connected faces share more than one edge. Faces ";
481  errstrm << (m_faces[f])->GetFid() << ", " << (m_faces[f+1])->GetFid();
482  ASSERTL0(false, errstrm.str());
483  }
484  }
485 
486  // Finally, set up the 1 top edge
487  check = 0;
488  for(i = 0; i < 4; i++)
489  {
490  for(j = 0; j < 4; j++)
491  {
492  if( (m_faces[2])->GetEid(i) == (m_faces[4])->GetEid(j) )
493  {
494  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[2])->GetEdge(i));
495  m_edges.push_back(edge);
496  check++;
497  }
498  }
499  }
500 
501  if( check < 1 )
502  {
503  std::ostringstream errstrm;
504  errstrm << "Connected faces do not share an edge. Faces ";
505  errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[3])->GetFid();
506  ASSERTL0(false, errstrm.str());
507  }
508  else if( check > 1)
509  {
510  std::ostringstream errstrm;
511  errstrm << "Connected faces share more than one edge. Faces ";
512  errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[3])->GetFid();
513  ASSERTL0(false, errstrm.str());
514  }
515  };
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
int GetEid(int i) const
Return the ID of edge i in this element.
Definition: Geometry3D.cpp:74
int GetFid(int i) const
Definition: Geometry.h:329
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
const Geometry1DSharedPtr GetEdge(int i) const
void Nektar::SpatialDomains::PrismGeom::SetUpLocalVertices ( )
private

Definition at line 518 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().

518  {
519 
520  // Set up the first 2 vertices (i.e. vertex 0,1)
521  if( ( m_edges[0]->GetVid(0) == m_edges[1]->GetVid(0) ) ||
522  ( m_edges[0]->GetVid(0) == m_edges[1]->GetVid(1) ) )
523  {
524  m_verts.push_back(m_edges[0]->GetVertex(1));
525  m_verts.push_back(m_edges[0]->GetVertex(0));
526  }
527  else if( ( m_edges[0]->GetVid(1) == m_edges[1]->GetVid(0) ) ||
528  ( m_edges[0]->GetVid(1) == m_edges[1]->GetVid(1) ) )
529  {
530  m_verts.push_back(m_edges[0]->GetVertex(0));
531  m_verts.push_back(m_edges[0]->GetVertex(1));
532  }
533  else
534  {
535  std::ostringstream errstrm;
536  errstrm << "Connected edges do not share a vertex. Edges ";
537  errstrm << m_edges[0]->GetEid() << ", " << m_edges[1]->GetEid();
538  ASSERTL0(false, errstrm.str());
539  }
540 
541  // set up the other bottom vertices (i.e. vertex 2,3)
542  for(int i = 1; i < 3; i++)
543  {
544  if( m_edges[i]->GetVid(0) == m_verts[i]->GetVid() )
545  {
546  m_verts.push_back(m_edges[i]->GetVertex(1));
547  }
548  else if( m_edges[i]->GetVid(1) == m_verts[i]->GetVid() )
549  {
550  m_verts.push_back(m_edges[i]->GetVertex(0));
551  }
552  else
553  {
554  std::ostringstream errstrm;
555  errstrm << "Connected edges do not share a vertex. Edges ";
556  errstrm << m_edges[i]->GetEid() << ", " << m_edges[i-1]->GetEid();
557  ASSERTL0(false, errstrm.str());
558  }
559  }
560 
561  // set up top vertices
562  // First, set up vertices 4,5
563  if( (m_edges[8]->GetVid(0) == m_edges[4]->GetVid(0)) ||
564  (m_edges[8]->GetVid(0) == m_edges[4]->GetVid(1)) )
565  {
566  m_verts.push_back(m_edges[8]->GetVertex(0));
567  m_verts.push_back(m_edges[8]->GetVertex(1));
568  }
569  else if( (m_edges[8]->GetVid(1) == m_edges[4]->GetVid(0)) ||
570  (m_edges[8]->GetVid(1) == m_edges[4]->GetVid(1)) )
571  {
572  m_verts.push_back(m_edges[8]->GetVertex(1));
573  m_verts.push_back(m_edges[8]->GetVertex(0));
574  }
575  else
576  {
577  std::ostringstream errstrm;
578  errstrm << "Connected edges do not share a vertex. Edges ";
579  errstrm << m_edges[8]->GetEid();
580  ASSERTL0(false, errstrm.str());
581  }
582  };
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
PointGeomSharedPtr GetVertex(int i) const
Definition: Geometry.h:348
int GetVid(int i) const
Definition: Geometry.h:319
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 862 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().

863  {
864  vector<int> tmp;
865 
866  int order0, order1;
867 
868  if (m_forient[0] < 9)
869  {
870  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
871  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
872  order0 = *max_element(tmp.begin(), tmp.end());
873  }
874  else
875  {
876  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
877  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(3));
878  order0 = *max_element(tmp.begin(), tmp.end());
879  }
880 
881  if (m_forient[0] < 9)
882  {
883  tmp.clear();
884  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
885  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(3));
886  tmp.push_back(m_faces[2]->GetXmap()->GetEdgeNcoeffs(2));
887  order1 = *max_element(tmp.begin(), tmp.end());
888  }
889  else
890  {
891  tmp.clear();
892  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
893  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
894  tmp.push_back(m_faces[2]->GetXmap()->GetEdgeNcoeffs(2));
895  order1 = *max_element(tmp.begin(), tmp.end());
896  }
897 
898  tmp.clear();
899  tmp.push_back(order0);
900  tmp.push_back(order1);
901  tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(1));
902  tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(2));
903  tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(1));
904  tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(2));
905  int order2 = *max_element(tmp.begin(), tmp.end());
906 
907  const LibUtilities::BasisKey A(
909  LibUtilities::PointsKey(
911  const LibUtilities::BasisKey B(
913  LibUtilities::PointsKey(
915  const LibUtilities::BasisKey C(
917  LibUtilities::PointsKey(
919 
921  A, B, C);
922  }
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Principle Modified Functions .
Definition: BasisType.h:49
StdRegions::StdExpansionSharedPtr GetXmap() const
Definition: Geometry.h:383
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:93
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
Principle Modified Functions .
Definition: BasisType.h:50
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
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 124 of file PrismGeom.cpp.

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

Referenced by v_ContainsPoint().

126  {
127  Array<OneD,NekDouble> locCoord(GetCoordim(),0.0);
128  return v_ContainsPoint(gloCoord,locCoord,tol);
129  }
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...
Definition: PrismGeom.cpp:124
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 131 of file PrismGeom.cpp.

References v_ContainsPoint().

135  {
136  NekDouble resid;
137  return v_ContainsPoint(gloCoord,locCoord,tol,resid);
138  }
double NekDouble
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...
Definition: PrismGeom.cpp:124
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 144 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().

149  {
150  // Validation checks
151  ASSERTL1(gloCoord.num_elements() == 3,
152  "Three dimensional geometry expects three coordinates.");
153 
154  // find min, max point and check if within twice this
155  // distance other false this is advisable since
156  // GetLocCoord is expensive for non regular elements.
157  if(GetMetricInfo()->GetGtype() != eRegular)
158  {
159  int i;
160  Array<OneD, NekDouble> mincoord(3), maxcoord(3);
161  NekDouble diff = 0.0;
162 
163  v_FillGeom();
164 
165  const int npts = m_xmap->GetTotPoints();
166  Array<OneD, NekDouble> pts(npts);
167 
168  for(i = 0; i < 3; ++i)
169  {
170  m_xmap->BwdTrans(m_coeffs[i], pts);
171 
172  mincoord[i] = Vmath::Vmin(pts.num_elements(),pts,1);
173  maxcoord[i] = Vmath::Vmax(pts.num_elements(),pts,1);
174 
175  diff = max(maxcoord[i] - mincoord[i],diff);
176  }
177 
178  for(i = 0; i < 3; ++i)
179  {
180  if((gloCoord[i] < mincoord[i] - 0.2*diff)||
181  (gloCoord[i] > maxcoord[i] + 0.2*diff))
182  {
183  return false;
184  }
185  }
186  }
187 
188  // Convert to the local Cartesian coordinates.
189  resid = v_GetLocCoords(gloCoord, locCoord);
190 
191  // Check local coordinate is within std region bounds.
192  if (locCoord[0] >= -(1+tol) && locCoord[1] >= -(1+tol) &&
193  locCoord[2] >= -(1+tol) && locCoord[1] <= (1+tol) &&
194  locCoord[0] + locCoord[2] <= tol)
195  {
196  return true;
197  }
198 
199  // If out of range clamp locCoord to be within [-1,1]^3
200  // since any larger value will be very oscillatory if
201  // called by 'returnNearestElmt' option in
202  // ExpList::GetExpIndex
203  for(int i = 0; i < 3; ++i)
204  {
205  if(locCoord[i] <-(1+tol))
206  {
207  locCoord[i] = -(1+tol);
208  }
209 
210  if(locCoord[i] > (1+tol))
211  {
212  locCoord[i] = 1+tol;
213  }
214  }
215 
216  return false;
217  }
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:765
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:857
virtual NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Definition: PrismGeom.cpp:278
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
Definition: Geometry3D.cpp:247
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Definition: Geometry.h:180
Geometry is straight-sided with constant geometric factors.
GeomFactorsSharedPtr GetMetricInfo()
Definition: Geometry.h:299
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Nektar::SpatialDomains::PrismGeom::v_GenGeomFactors ( )
protectedvirtual

Generate the geometry factors for this element.

Reimplemented from Nektar::SpatialDomains::Geometry3D.

Definition at line 219 of file PrismGeom.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::SpatialDomains::eDeformed, Nektar::SpatialDomains::ePtsFilled, Nektar::SpatialDomains::eRegular, Nektar::NekConstants::kNekZeroTol, Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_coordim, Nektar::SpatialDomains::Geometry::m_geomFactors, Nektar::SpatialDomains::Geometry::m_geomFactorsState, Nektar::SpatialDomains::Geometry3D::m_verts, Nektar::SpatialDomains::Geometry::m_xmap, and Nektar::SpatialDomains::Geometry3D::v_FillGeom().

220  {
222  {
223  int i,f;
224  GeomType Gtype = eRegular;
225 
226  v_FillGeom();
227 
228  // check to see if expansions are linear
229  for(i = 0; i < m_coordim; ++i)
230  {
231  if (m_xmap->GetBasisNumModes(0) != 2 ||
232  m_xmap->GetBasisNumModes(1) != 2 ||
233  m_xmap->GetBasisNumModes(2) != 2 )
234  {
235  Gtype = eDeformed;
236  }
237  }
238 
239  // check to see if all quadrilateral faces are parallelograms
240  if(Gtype == eRegular)
241  {
242  // Vertex ids of quad faces
243  const unsigned int faceVerts[3][4] =
244  { {0,1,2,3} ,
245  {1,2,5,4} ,
246  {0,3,5,4} };
247 
248  for(f = 0; f < 3; f++)
249  {
250  // Ensure each face is a parallelogram? Check this.
251  for (i = 0; i < m_coordim; i++)
252  {
253  if( fabs( (*m_verts[ faceVerts[f][0] ])(i) -
254  (*m_verts[ faceVerts[f][1] ])(i) +
255  (*m_verts[ faceVerts[f][2] ])(i) -
256  (*m_verts[ faceVerts[f][3] ])(i) )
258  {
259  Gtype = eDeformed;
260  break;
261  }
262  }
263 
264  if (Gtype == eDeformed)
265  {
266  break;
267  }
268  }
269  }
270 
272  Gtype, m_coordim, m_xmap, m_coeffs);
274  }
275  }
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
GeomFactorsSharedPtr m_geomFactors
Definition: Geometry.h:170
static const NekDouble kNekZeroTol
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
Definition: Geometry3D.cpp:247
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Definition: Geometry.h:180
Geometry is straight-sided with constant geometric factors.
Geometric information has been generated.
GeomType
Indicates the type of element geometry.
Geometry is curved or has non-constant factors.
int m_coordim
coordinate dimension
Definition: Geometry.h:169
int Nektar::SpatialDomains::PrismGeom::v_GetDir ( const int  faceidx,
const int  facedir 
) const
protectedvirtual

Implements Nektar::SpatialDomains::Geometry3D.

Definition at line 104 of file PrismGeom.cpp.

105  {
106  if (faceidx == 0)
107  {
108  return facedir;
109  }
110  else if (faceidx == 1 || faceidx == 3)
111  {
112  return 2 * facedir;
113  }
114  else
115  {
116  return 1 + facedir;
117  }
118  }
int Nektar::SpatialDomains::PrismGeom::v_GetEdgeFaceMap ( const int  i,
const int  j 
) const
protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 377 of file PrismGeom.cpp.

References EdgeFaceConnectivity.

378  {
379  return EdgeFaceConnectivity[i][j];
380  }
static const unsigned int EdgeFaceConnectivity[9][2]
Definition: PrismGeom.h:106
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 278 of file PrismGeom.cpp.

References Nektar::SpatialDomains::PointGeom::dist(), Nektar::SpatialDomains::PointGeom::dot(), Nektar::SpatialDomains::eRegular, Nektar::SpatialDomains::Geometry::GetMetricInfo(), Vmath::Imin(), Nektar::SpatialDomains::Geometry::m_coeffs, Nektar::SpatialDomains::Geometry::m_coordim, Nektar::SpatialDomains::Geometry3D::m_verts, Nektar::SpatialDomains::Geometry::m_xmap, Nektar::SpatialDomains::PointGeom::Mult(), Nektar::SpatialDomains::Geometry3D::NewtonIterationForLocCoord(), npts, Vmath::Sadd(), Nektar::SpatialDomains::PointGeom::Sub(), Nektar::SpatialDomains::Geometry3D::v_FillGeom(), Vmath::Vmul(), and Vmath::Vvtvp().

Referenced by v_ContainsPoint().

281  {
282  NekDouble ptdist = 1e6;
283 
284  // calculate local coordinate for coord
285  if(GetMetricInfo()->GetGtype() == eRegular)
286  {
287  // Point inside tetrahedron
288  PointGeom r(m_coordim, 0, coords[0], coords[1], coords[2]);
289 
290  // Edges
291  PointGeom er0, e10, e30, e40;
292  er0.Sub(r,*m_verts[0]);
293  e10.Sub(*m_verts[1],*m_verts[0]);
294  e30.Sub(*m_verts[3],*m_verts[0]);
295  e40.Sub(*m_verts[4],*m_verts[0]);
296 
297  // Cross products (Normal times area)
298  PointGeom cp1030, cp3040, cp4010;
299  cp1030.Mult(e10,e30);
300  cp3040.Mult(e30,e40);
301  cp4010.Mult(e40,e10);
302 
303  // Barycentric coordinates (relative volume)
304  NekDouble V = e40.dot(cp1030); // Prism Volume = {(e40)dot(e10)x(e30)}/2
305  NekDouble beta = er0.dot(cp3040) / (2.0*V); // volume1 = {(er0)dot(e30)x(e40)}/4
306  NekDouble gamma = er0.dot(cp4010) / (3.0*V); // volume2 = {(er0)dot(e40)x(e10)}/6
307  NekDouble delta = er0.dot(cp1030) / (2.0*V); // volume3 = {(er0)dot(e10)x(e30)}/4
308 
309  // Make Prism bigger
310  Lcoords[0] = 2.0*beta - 1.0;
311  Lcoords[1] = 2.0*gamma - 1.0;
312  Lcoords[2] = 2.0*delta - 1.0;
313 
314  // Set ptdist to distance to nearest vertex
315  for(int i = 0; i < 6; ++i)
316  {
317  ptdist = min(ptdist,r.dist(*m_verts[i]));
318  }
319  }
320  else
321  {
322  v_FillGeom();
323 
324  // Determine nearest point of coords to values in m_xmap
325  int npts = m_xmap->GetTotPoints();
326  Array<OneD, NekDouble> ptsx(npts), ptsy(npts), ptsz(npts);
327  Array<OneD, NekDouble> tmp1(npts), tmp2(npts);
328 
329  m_xmap->BwdTrans(m_coeffs[0], ptsx);
330  m_xmap->BwdTrans(m_coeffs[1], ptsy);
331  m_xmap->BwdTrans(m_coeffs[2], ptsz);
332 
333  const Array<OneD, const NekDouble> za = m_xmap->GetPoints(0);
334  const Array<OneD, const NekDouble> zb = m_xmap->GetPoints(1);
335  const Array<OneD, const NekDouble> zc = m_xmap->GetPoints(2);
336 
337  //guess the first local coords based on nearest point
338  Vmath::Sadd(npts, -coords[0], ptsx,1,tmp1,1);
339  Vmath::Vmul (npts, tmp1,1,tmp1,1,tmp1,1);
340  Vmath::Sadd(npts, -coords[1], ptsy,1,tmp2,1);
341  Vmath::Vvtvp(npts, tmp2,1,tmp2,1,tmp1,1,tmp1,1);
342  Vmath::Sadd(npts, -coords[2], ptsz,1,tmp2,1);
343  Vmath::Vvtvp(npts, tmp2,1,tmp2,1,tmp1,1,tmp1,1);
344 
345  int min_i = Vmath::Imin(npts,tmp1,1);
346 
347  // distance from coordinate to nearest point for return value.
348  ptdist = sqrt(tmp1[min_i]);
349 
350  // Get collapsed coordinate
351  int qa = za.num_elements(), qb = zb.num_elements();
352  Lcoords[2] = zc[min_i/(qa*qb)];
353  min_i = min_i%(qa*qb);
354  Lcoords[1] = zb[min_i/qa];
355  Lcoords[0] = za[min_i%qa];
356 
357  // recover cartesian coordinate from collapsed coordinate.
358  Lcoords[0] = (1.0+Lcoords[0])*(1.0-Lcoords[2])/2 - 1.0;
359 
360  // Perform newton iteration to find local coordinates
361  NekDouble resid = 0.0;
362  NewtonIterationForLocCoord(coords, ptsx, ptsy, ptsz, Lcoords,resid);
363  }
364  return ptdist;
365  }
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:833
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
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)
Definition: Geometry3D.cpp:99
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
Definition: Geometry3D.cpp:247
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Definition: Geometry.h:180
Geometry is straight-sided with constant geometric factors.
GeomFactorsSharedPtr GetMetricInfo()
Definition: Geometry.h:299
int m_coordim
coordinate dimension
Definition: Geometry.h:169
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
int Nektar::SpatialDomains::PrismGeom::v_GetNumEdges ( ) const
protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 94 of file PrismGeom.cpp.

95  {
96  return 9;
97  }
int Nektar::SpatialDomains::PrismGeom::v_GetNumFaces ( ) const
protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 99 of file PrismGeom.cpp.

100  {
101  return 5;
102  }
int Nektar::SpatialDomains::PrismGeom::v_GetNumVerts ( ) const
protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 89 of file PrismGeom.cpp.

90  {
91  return 6;
92  }
int Nektar::SpatialDomains::PrismGeom::v_GetVertexEdgeMap ( const int  i,
const int  j 
) const
protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 367 of file PrismGeom.cpp.

References VertexEdgeConnectivity.

368  {
369  return VertexEdgeConnectivity[i][j];
370  }
static const unsigned int VertexEdgeConnectivity[6][3]
Definition: PrismGeom.h:104
int Nektar::SpatialDomains::PrismGeom::v_GetVertexFaceMap ( const int  i,
const int  j 
) const
protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 372 of file PrismGeom.cpp.

References VertexFaceConnectivity.

373  {
374  return VertexFaceConnectivity[i][j];
375  }
static const unsigned int VertexFaceConnectivity[6][3]
Definition: PrismGeom.h:105
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 843 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().

846  {
847  Geometry::v_Reset(curvedEdges, curvedFaces);
848 
849  for (int i = 0; i < 5; ++i)
850  {
851  m_faces[i]->Reset(curvedEdges, curvedFaces);
852  }
853 
854  SetUpXmap();
855  SetUpCoeffs(m_xmap->GetNcoeffs());
856  }
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
Definition: Geometry.cpp:307
void SetUpXmap()
Set up the m_xmap object by determining the order of each direction from derived faces.
Definition: PrismGeom.cpp:862
void SetUpCoeffs(const int nCoeffs)
Initialise the m_coeffs array.
Definition: Geometry.h:484

Member Data Documentation

const unsigned int Nektar::SpatialDomains::PrismGeom::EdgeFaceConnectivity
staticprivate
Initial value:
= {
{0,1},{0,2},{0,3},{0,4},{1,4},{1,2},{2,3},{3,4},{2,4}}

Definition at line 106 of file PrismGeom.h.

Referenced by v_GetEdgeFaceMap().

const int Nektar::SpatialDomains::PrismGeom::kNedges = 9
static

Definition at line 58 of file PrismGeom.h.

Referenced by PrismGeom(), and SetUpEdgeOrientation().

const int Nektar::SpatialDomains::PrismGeom::kNfaces
static
const int Nektar::SpatialDomains::PrismGeom::kNqfaces = 3
static
const int Nektar::SpatialDomains::PrismGeom::kNtfaces = 2
static
const int Nektar::SpatialDomains::PrismGeom::kNverts = 6
static

Definition at line 57 of file PrismGeom.h.

const unsigned int Nektar::SpatialDomains::PrismGeom::VertexEdgeConnectivity
staticprivate
Initial value:
= {
{0,3,4},{0,1,5},{1,2,6},{2,3,7},{4,5,8},{6,7,8}}

Definition at line 104 of file PrismGeom.h.

Referenced by v_GetVertexEdgeMap().

const unsigned int Nektar::SpatialDomains::PrismGeom::VertexFaceConnectivity
staticprivate
Initial value:
= {
{0,1,4},{0,1,2},{0,2,3},{0,3,4},{1,2,4},{2,3,4}}

Definition at line 105 of file PrismGeom.h.

Referenced by v_GetVertexFaceMap().

const std::string Nektar::SpatialDomains::PrismGeom::XMLElementType
static

Definition at line 63 of file PrismGeom.h.