Nektar++
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:
[legend]

Public Member Functions

 PrismGeom ()
 
 PrismGeom (int id, const Geometry2DSharedPtr faces[])
 
 ~PrismGeom ()
 
- Public Member Functions inherited from Nektar::SpatialDomains::Geometry3D
 Geometry3D ()
 
 Geometry3D (const int coordim)
 
virtual ~Geometry3D ()
 
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 ()
 Default constructor. More...
 
 Geometry (int coordim)
 Constructor when supplied a coordinate dimension. More...
 
virtual ~Geometry ()
 Default destructor. More...
 
int GetCoordim () const
 Return the coordinate dimension of this object (i.e. the dimension of the space in which this object is embedded). More...
 
void SetCoordim (int coordim)
 Sets the coordinate dimension of this object (i.e. the dimension of the space in which this object is embedded). More...
 
GeomFactorsSharedPtr GetGeomFactors ()
 Get the geometric factors for this object, generating them if required. More...
 
GeomFactorsSharedPtr GetRefGeomFactors (const Array< OneD, const LibUtilities::BasisSharedPtr > &tbasis)
 
GeomFactorsSharedPtr GetMetricInfo ()
 Get the geometric factors for this object. More...
 
LibUtilities::ShapeType GetShapeType (void)
 Get the geometric shape type of this object. More...
 
int GetGlobalID (void) const
 Get the ID of this object. More...
 
void SetGlobalID (int globalid)
 Set the ID of this object. More...
 
int GetVid (int i) const
 Get the ID of vertex i of this object. More...
 
int GetEid (int i) const
 Get the ID of edge i of this object. More...
 
int GetFid (int i) const
 Get the ID of face i of this object. More...
 
int GetTid (int i) const
 Get the ID of trace i of this object. More...
 
PointGeomSharedPtr GetVertex (int i) const
 Returns vertex i of this object. More...
 
Geometry1DSharedPtr GetEdge (int i) const
 Returns edge i of this object. More...
 
Geometry2DSharedPtr GetFace (int i) const
 Returns face i of this object. More...
 
StdRegions::Orientation GetEorient (const int i) const
 Returns the orientation of edge i with respect to the ordering of edges in the standard element. More...
 
StdRegions::Orientation GetForient (const int i) const
 Returns the orientation of face i with respect to the ordering of faces in the standard element. More...
 
int GetNumVerts () const
 Get the number of vertices of this object. More...
 
int GetNumEdges () const
 Get the number of edges of this object. More...
 
int GetNumFaces () const
 Get the number of faces of this object. More...
 
int GetShapeDim () const
 Get the object's shape dimension. More...
 
StdRegions::StdExpansionSharedPtr GetXmap () const
 Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standard element to physical element. More...
 
const Array< OneD, const NekDouble > & GetCoeffs (const int i) const
 Return the coefficients of the transformation Geometry::m_xmap in coordinate direction i. More...
 
void FillGeom ()
 Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements. More...
 
std::array< NekDouble, 6 > GetBoundingBox ()
 Generates the bounding box for the element. More...
 
bool ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0)
 Determine whether an element contains a particular Cartesian coordinate \((x,y,z)\). More...
 
bool ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol)
 Determine whether an element contains a particular Cartesian coordinate \((x,y,z)\). More...
 
bool ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol, NekDouble &resid)
 Determine whether an element contains a particular Cartesian coordinate \(\vec{x} = (x,y,z)\). More...
 
NekDouble GetLocCoords (const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
 Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this geometry object. More...
 
NekDouble GetCoord (const int i, const Array< OneD, const NekDouble > &Lcoord)
 Given local collapsed coordinate Lcoord, return the value of physical coordinate in direction i. More...
 
bool MinMaxCheck (const Array< OneD, const NekDouble > &gloCoord)
 Check if given global coord is within twice the min/max distance of the element. More...
 
void ClampLocCoords (Array< OneD, NekDouble > &locCoord, NekDouble tol)
 Clamp local coords to be within standard regions [-1, 1]^dim. More...
 
int GetVertexEdgeMap (int i, int j) const
 Returns the standard element edge IDs that are connected to a given vertex. More...
 
int GetVertexFaceMap (int i, int j) const
 Returns the standard element face IDs that are connected to a given vertex. More...
 
int GetEdgeFaceMap (int i, int j) const
 Returns the standard element edge IDs that are connected to a given face. More...
 
void Reset (CurveMap &curvedEdges, CurveMap &curvedFaces)
 Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated GeomFactors. More...
 
void Setup ()
 

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 = kNqfaces + kNtfaces
 
static const std::string XMLElementType
 
- Static Public Attributes inherited from Nektar::SpatialDomains::Geometry3D
static const int kDim = 3
 

Protected Member Functions

virtual void v_GenGeomFactors ()
 
virtual NekDouble v_GetLocCoords (const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
 Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this geometry object. More...
 
virtual bool v_ContainsPoint (const Array< OneD, const NekDouble > &gloCoord, Array< OneD, NekDouble > &locCoord, NekDouble tol, NekDouble &resid)
 
virtual int v_GetVertexEdgeMap (const int i, const int j) const
 Returns the standard element edge IDs that are connected to a given vertex. More...
 
virtual int v_GetVertexFaceMap (const int i, const int j) const
 Returns the standard element face IDs that are connected to a given vertex. More...
 
virtual int v_GetEdgeFaceMap (const int i, const int j) const
 Returns the standard element edge IDs that are connected to a given face. More...
 
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, zero Geometry::m_coeffs and remove allocated GeomFactors. More...
 
virtual void v_Setup ()
 
- 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
 Get the object's shape dimension. More...
 
virtual int v_GetNumVerts () const
 Get the number of vertices of this object. More...
 
virtual int v_GetNumEdges () const
 Get the number of edges of this object. More...
 
virtual int v_GetNumFaces () const
 Get the number of faces of this object. More...
 
virtual PointGeomSharedPtr v_GetVertex (int i) const
 
virtual Geometry1DSharedPtr v_GetEdge (int i) const
 Returns edge i of this object. More...
 
virtual Geometry2DSharedPtr v_GetFace (int i) const
 Returns face i of this object. More...
 
virtual StdRegions::Orientation v_GetEorient (const int i) const
 Returns the orientation of edge i with respect to the ordering of edges in the standard element. More...
 
virtual StdRegions::Orientation v_GetForient (const int i) const
 Returns the orientation of face i with respect to the ordering of faces in the standard element. More...
 
- Protected Member Functions inherited from Nektar::SpatialDomains::Geometry
void GenGeomFactors ()
 Handles generation of geometry factors. More...
 
virtual StdRegions::StdExpansionSharedPtr v_GetXmap () const
 Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standard element to physical element. More...
 
void SetUpCoeffs (const int nCoeffs)
 Initialise the Geometry::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)
 Check to see if a geometric factor has already been created that contains the same regular information. More...
 
- Protected Attributes inherited from Nektar::SpatialDomains::Geometry3D
PointGeomVector m_verts
 
SegGeomVector m_edges
 
Geometry2DVector m_faces
 
std::vector< StdRegions::Orientationm_eorient
 
std::vector< StdRegions::Orientationm_forient
 
int m_eid
 
bool m_ownverts
 
- Protected Attributes inherited from Nektar::SpatialDomains::Geometry
int m_coordim
 Coordinate dimension of this geometry object. More...
 
GeomFactorsSharedPtr m_geomFactors
 Geometric factors. More...
 
GeomState m_geomFactorsState
 State of the geometric factors. More...
 
StdRegions::StdExpansionSharedPtr m_xmap
 \(\chi\) mapping containing isoparametric transformation. More...
 
GeomState m_state
 Enumeration to dictate whether coefficients are filled. More...
 
bool m_setupState
 Wether or not the setup routines have been run. More...
 
GeomType m_geomType
 Type of geometry. More...
 
LibUtilities::ShapeType m_shapeType
 Type of shape. More...
 
int m_globalID
 Global ID. More...
 
Array< OneD, Array< OneD, NekDouble > > m_coeffs
 Array containing expansion coefficients of m_xmap. More...
 
- Static Protected Attributes inherited from Nektar::SpatialDomains::Geometry
static GeomFactorsVector m_regGeomFactorsManager
 

Detailed Description

Definition at line 47 of file PrismGeom.h.

Constructor & Destructor Documentation

◆ PrismGeom() [1/2]

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

Definition at line 56 of file PrismGeom.cpp.

References Nektar::LibUtilities::ePrism.

57 {
59 }
LibUtilities::ShapeType m_shapeType
Type of shape.
Definition: Geometry.h:197

◆ PrismGeom() [2/2]

Nektar::SpatialDomains::PrismGeom::PrismGeom ( int  id,
const Geometry2DSharedPtr  faces[] 
)

Copy the face shared pointers.

Set up orientation vectors with correct amount of elements.

Set up local objects.

Definition at line 61 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_globalID, Nektar::SpatialDomains::Geometry::m_shapeType, SetUpEdgeOrientation(), SetUpFaceOrientation(), SetUpLocalEdges(), and SetUpLocalVertices().

62  : Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
63 {
65  m_globalID = id;
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 }
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry3D.h:86
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:87
int GetCoordim() const
Return the coordinate dimension of this object (i.e. the dimension of the space in which this object ...
Definition: Geometry.h:271
Geometry1DSharedPtr GetEdge(int i) const
Returns edge i of this object.
Definition: Geometry.h:351
LibUtilities::ShapeType m_shapeType
Type of shape.
Definition: Geometry.h:197
PointGeomSharedPtr GetVertex(int i) const
Returns vertex i of this object.
Definition: Geometry.h:343

◆ ~PrismGeom()

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

Definition at line 81 of file PrismGeom.cpp.

82 {
83 }

Member Function Documentation

◆ SetUpEdgeOrientation()

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

Definition at line 515 of file PrismGeom.cpp.

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

Referenced by PrismGeom().

516 {
517 
518  // This 2D array holds the local id's of all the vertices
519  // for every edge. For every edge, they are ordered to what we
520  // define as being Forwards
521  const unsigned int edgeVerts[kNedges][2] = {
522  {0, 1}, {1, 2}, {3, 2}, {0, 3}, {0, 4}, {1, 4}, {2, 5}, {3, 5}, {4, 5}};
523 
524  int i;
525  for (i = 0; i < kNedges; i++)
526  {
527  if (m_edges[i]->GetVid(0) == m_verts[edgeVerts[i][0]]->GetGlobalID())
528  {
530  }
531  else if (m_edges[i]->GetVid(0) == m_verts[edgeVerts[i][1]]->GetGlobalID())
532  {
534  }
535  else
536  {
537  ASSERTL0(false, "Could not find matching vertex for the edge");
538  }
539  }
540 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry3D.h:86
int GetGlobalID(void) const
Get the ID of this object.
Definition: Geometry.h:314
int GetVid(int i) const
Get the ID of vertex i of this object.
Definition: Geometry.cpp:137

◆ SetUpFaceOrientation()

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

Definition at line 542 of file PrismGeom.cpp.

References ASSERTL0, ASSERTL1, Nektar::SpatialDomains::Geometry::GetGlobalID(), 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().

543 {
544  int f, i;
545 
546  // These arrays represent the vector of the A and B
547  // coordinate of the local elemental coordinate system
548  // where A corresponds with the coordinate direction xi_i
549  // with the lowest index i (for that particular face)
550  // Coordinate 'B' then corresponds to the other local
551  // coordinate (i.e. with the highest index)
552  Array<OneD, NekDouble> elementAaxis(m_coordim);
553  Array<OneD, NekDouble> elementBaxis(m_coordim);
554 
555  // These arrays correspond to the local coordinate
556  // system of the face itself (i.e. the Geometry2D)
557  // faceAaxis correspond to the xi_0 axis
558  // faceBaxis correspond to the xi_1 axis
559  Array<OneD, NekDouble> faceAaxis(m_coordim);
560  Array<OneD, NekDouble> faceBaxis(m_coordim);
561 
562  // This is the base vertex of the face (i.e. the Geometry2D)
563  // This corresponds to thevertex with local ID 0 of the
564  // Geometry2D
565  unsigned int baseVertex;
566 
567  // The lenght of the vectors above
568  NekDouble elementAaxis_length;
569  NekDouble elementBaxis_length;
570  NekDouble faceAaxis_length;
571  NekDouble faceBaxis_length;
572 
573  // This 2D array holds the local id's of all the vertices
574  // for every face. For every face, they are ordered in such
575  // a way that the implementation below allows a unified approach
576  // for all faces.
577  const unsigned int faceVerts[kNfaces][QuadGeom::kNverts] = {
578  {0, 1, 2, 3},
579  {0, 1, 4, 0}, // This is triangle requires only three vertices
580  {1, 2, 5, 4},
581  {3, 2, 5, 0}, // This is triangle requires only three vertices
582  {0, 3, 5, 4},
583  };
584 
585  NekDouble dotproduct1 = 0.0;
586  NekDouble dotproduct2 = 0.0;
587 
588  unsigned int orientation;
589 
590  // Loop over all the faces to set up the orientation
591  for (f = 0; f < kNqfaces + kNtfaces; f++)
592  {
593  // initialisation
594  elementAaxis_length = 0.0;
595  elementBaxis_length = 0.0;
596  faceAaxis_length = 0.0;
597  faceBaxis_length = 0.0;
598 
599  dotproduct1 = 0.0;
600  dotproduct2 = 0.0;
601 
602  baseVertex = m_faces[f]->GetVid(0);
603 
604  // We are going to construct the vectors representing the A and B axis
605  // of every face. These vectors will be constructed as a
606  // vector-representation
607  // of the edges of the face. However, for both coordinate directions, we
608  // can
609  // represent the vectors by two different edges. That's why we need to
610  // make sure that
611  // we pick the edge to which the baseVertex of the
612  // Geometry2D-representation of the face
613  // belongs...
614 
615  // Compute the length of edges on a base-face
616  if (f == 1 || f == 3)
617  { // Face is a Triangle
618  if (baseVertex == m_verts[faceVerts[f][0]]->GetVid())
619  {
620  for (i = 0; i < m_coordim; i++)
621  {
622  elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
623  (*m_verts[faceVerts[f][0]])[i];
624  elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
625  (*m_verts[faceVerts[f][0]])[i];
626  }
627  }
628  else if (baseVertex == m_verts[faceVerts[f][1]]->GetVid())
629  {
630  for (i = 0; i < m_coordim; i++)
631  {
632  elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
633  (*m_verts[faceVerts[f][0]])[i];
634  elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
635  (*m_verts[faceVerts[f][1]])[i];
636  }
637  }
638  else if (baseVertex == m_verts[faceVerts[f][2]]->GetVid())
639  {
640  for (i = 0; i < m_coordim; i++)
641  {
642  elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
643  (*m_verts[faceVerts[f][2]])[i];
644  elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
645  (*m_verts[faceVerts[f][0]])[i];
646  }
647  }
648  else
649  {
650  ASSERTL0(false, "Could not find matching vertex for the face");
651  }
652  }
653  else
654  { // Face is a Quad
655  if (baseVertex == m_verts[faceVerts[f][0]]->GetGlobalID())
656  {
657  for (i = 0; i < m_coordim; i++)
658  {
659  elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
660  (*m_verts[faceVerts[f][0]])[i];
661  elementBaxis[i] = (*m_verts[faceVerts[f][3]])[i] -
662  (*m_verts[faceVerts[f][0]])[i];
663  }
664  }
665  else if (baseVertex == m_verts[faceVerts[f][1]]->GetGlobalID())
666  {
667  for (i = 0; i < m_coordim; i++)
668  {
669  elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
670  (*m_verts[faceVerts[f][0]])[i];
671  elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
672  (*m_verts[faceVerts[f][1]])[i];
673  }
674  }
675  else if (baseVertex == m_verts[faceVerts[f][2]]->GetGlobalID())
676  {
677  for (i = 0; i < m_coordim; i++)
678  {
679  elementAaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
680  (*m_verts[faceVerts[f][3]])[i];
681  elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
682  (*m_verts[faceVerts[f][1]])[i];
683  }
684  }
685  else if (baseVertex == m_verts[faceVerts[f][3]]->GetGlobalID())
686  {
687  for (i = 0; i < m_coordim; i++)
688  {
689  elementAaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
690  (*m_verts[faceVerts[f][3]])[i];
691  elementBaxis[i] = (*m_verts[faceVerts[f][3]])[i] -
692  (*m_verts[faceVerts[f][0]])[i];
693  }
694  }
695  else
696  {
697  ASSERTL0(false, "Could not find matching vertex for the face");
698  }
699  }
700  // Now, construct the edge-vectors of the local coordinates of
701  // the Geometry2D-representation of the face
702  for (i = 0; i < m_coordim; i++)
703  {
704  int v = m_faces[f]->GetNumVerts() - 1;
705  faceAaxis[i] =
706  (*m_faces[f]->GetVertex(1))[i] - (*m_faces[f]->GetVertex(0))[i];
707  faceBaxis[i] =
708  (*m_faces[f]->GetVertex(v))[i] - (*m_faces[f]->GetVertex(0))[i];
709 
710  elementAaxis_length += pow(elementAaxis[i], 2);
711  elementBaxis_length += pow(elementBaxis[i], 2);
712  faceAaxis_length += pow(faceAaxis[i], 2);
713  faceBaxis_length += pow(faceBaxis[i], 2);
714  }
715 
716  elementAaxis_length = sqrt(elementAaxis_length);
717  elementBaxis_length = sqrt(elementBaxis_length);
718  faceAaxis_length = sqrt(faceAaxis_length);
719  faceBaxis_length = sqrt(faceBaxis_length);
720 
721  // Calculate the inner product of both the A-axis
722  // (i.e. Elemental A axis and face A axis)
723  for (i = 0; i < m_coordim; i++)
724  {
725  dotproduct1 += elementAaxis[i] * faceAaxis[i];
726  }
727 
728  orientation = 0;
729 
730  // if the innerproduct is equal to the (absolute value of the ) products
731  // of the lengths of both vectors, then, the coordinate systems will NOT
732  // be transposed
733  if (fabs(elementAaxis_length * faceAaxis_length - fabs(dotproduct1)) <
735  {
736  // if the inner product is negative, both A-axis point
737  // in reverse direction
738  if (dotproduct1 < 0.0)
739  {
740  orientation += 2;
741  }
742 
743  // calculate the inner product of both B-axis
744  for (i = 0; i < m_coordim; i++)
745  {
746  dotproduct2 += elementBaxis[i] * faceBaxis[i];
747  }
748 
749  ASSERTL1(
750  fabs(fabs(dotproduct2 / elementBaxis_length / faceBaxis_length)
751  - 1.0) < NekConstants::kNekZeroTol,
752  "These vectors should be parallel");
753 
754  // if the inner product is negative, both B-axis point
755  // in reverse direction
756  if (dotproduct2 < 0.0)
757  {
758  orientation++;
759  }
760  }
761  // The coordinate systems are transposed
762  else
763  {
764  orientation = 4;
765 
766  // Calculate the inner product between the elemental A-axis
767  // and the B-axis of the face (which are now the corresponding axis)
768  dotproduct1 = 0.0;
769  for (i = 0; i < m_coordim; i++)
770  {
771  dotproduct1 += elementAaxis[i] * faceBaxis[i];
772  }
773 
774  // check that both these axis are indeed parallel
775  ASSERTL1(
776  fabs(fabs(dotproduct1) / elementAaxis_length / faceBaxis_length
777  - 1.0) < NekConstants::kNekZeroTol,
778  "These vectors should be parallel");
779 
780  // if the result is negative, both axis point in reverse
781  // directions
782  if (dotproduct1 < 0.0)
783  {
784  orientation += 2;
785  }
786 
787  // Do the same for the other two corresponding axis
788  dotproduct2 = 0.0;
789  for (i = 0; i < m_coordim; i++)
790  {
791  dotproduct2 += elementBaxis[i] * faceAaxis[i];
792  }
793 
794  ASSERTL1(
795  fabs(fabs(dotproduct2) / elementBaxis_length / faceAaxis_length
796  - 1.0) < NekConstants::kNekZeroTol,
797  "These vectors should be parallel");
798 
799  if (dotproduct2 < 0.0)
800  {
801  orientation++;
802  }
803  }
804 
805  orientation = orientation + 5;
806  // Fill the m_forient array
807  m_forient[f] = (StdRegions::Orientation)orientation;
808  }
809 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
static const int kNverts
Definition: QuadGeom.h:77
int GetGlobalID(void) const
Get the ID of this object.
Definition: Geometry.h:314
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:87
static const NekDouble kNekZeroTol
double NekDouble
int GetVid(int i) const
Get the ID of vertex i of this object.
Definition: Geometry.cpp:137
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
PointGeomSharedPtr GetVertex(int i) const
Returns vertex i of this object.
Definition: Geometry.h:343
int m_coordim
Coordinate dimension of this geometry object.
Definition: Geometry.h:183

◆ SetUpLocalEdges()

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

Definition at line 299 of file PrismGeom.cpp.

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

Referenced by PrismGeom().

300 {
301  // find edge 0
302  int i, j;
303  unsigned int check;
304 
305  SegGeomSharedPtr edge;
306 
307  // First set up the 4 bottom edges
308  int f; // Connected face index
309  for (f = 1; f < 5; f++)
310  {
311  int nEdges = m_faces[f]->GetNumEdges();
312  check = 0;
313  for (i = 0; i < 4; i++)
314  {
315  for (j = 0; j < nEdges; j++)
316  {
317  if (m_faces[0]->GetEid(i) == m_faces[f]->GetEid(j))
318  {
319  edge = dynamic_pointer_cast<SegGeom>(
320  (m_faces[0])->GetEdge(i));
321  m_edges.push_back(edge);
322  check++;
323  }
324  }
325  }
326 
327  if (check < 1)
328  {
329  std::ostringstream errstrm;
330  errstrm << "Connected faces do not share an edge. Faces ";
331  errstrm << (m_faces[0])->GetGlobalID() << ", "
332  << (m_faces[f])->GetGlobalID();
333  ASSERTL0(false, errstrm.str());
334  }
335  else if (check > 1)
336  {
337  std::ostringstream errstrm;
338  errstrm << "Connected faces share more than one edge. Faces ";
339  errstrm << (m_faces[0])->GetGlobalID() << ", "
340  << (m_faces[f])->GetGlobalID();
341  ASSERTL0(false, errstrm.str());
342  }
343  }
344 
345  // Then, set up the 4 vertical edges
346  check = 0;
347  for (i = 0; i < 3; i++) // Set up the vertical edge :face(1) and face(4)
348  {
349  for (j = 0; j < 4; j++)
350  {
351  if ((m_faces[1])->GetEid(i) == (m_faces[4])->GetEid(j))
352  {
353  edge = dynamic_pointer_cast<SegGeom>(
354  (m_faces[1])->GetEdge(i));
355  m_edges.push_back(edge);
356  check++;
357  }
358  }
359  }
360  if (check < 1)
361  {
362  std::ostringstream errstrm;
363  errstrm << "Connected faces do not share an edge. Faces ";
364  errstrm << (m_faces[1])->GetGlobalID() << ", "
365  << (m_faces[4])->GetGlobalID();
366  ASSERTL0(false, errstrm.str());
367  }
368  else if (check > 1)
369  {
370  std::ostringstream errstrm;
371  errstrm << "Connected faces share more than one edge. Faces ";
372  errstrm << (m_faces[1])->GetGlobalID() << ", "
373  << (m_faces[4])->GetGlobalID();
374  ASSERTL0(false, errstrm.str());
375  }
376  // Set up vertical edges: face(1) through face(4)
377  for (f = 1; f < 4; f++)
378  {
379  check = 0;
380  for (i = 0; i < m_faces[f]->GetNumEdges(); i++)
381  {
382  for (j = 0; j < m_faces[f + 1]->GetNumEdges(); j++)
383  {
384  if ((m_faces[f])->GetEid(i) == (m_faces[f + 1])->GetEid(j))
385  {
386  edge = dynamic_pointer_cast<SegGeom>(
387  (m_faces[f])->GetEdge(i));
388  m_edges.push_back(edge);
389  check++;
390  }
391  }
392  }
393 
394  if (check < 1)
395  {
396  std::ostringstream errstrm;
397  errstrm << "Connected faces do not share an edge. Faces ";
398  errstrm << (m_faces[f])->GetGlobalID() << ", "
399  << (m_faces[f + 1])->GetGlobalID();
400  ASSERTL0(false, errstrm.str());
401  }
402  else if (check > 1)
403  {
404  std::ostringstream errstrm;
405  errstrm << "Connected faces share more than one edge. Faces ";
406  errstrm << (m_faces[f])->GetGlobalID() << ", "
407  << (m_faces[f + 1])->GetGlobalID();
408  ASSERTL0(false, errstrm.str());
409  }
410  }
411 
412  // Finally, set up the 1 top edge
413  check = 0;
414  for (i = 0; i < 4; i++)
415  {
416  for (j = 0; j < 4; j++)
417  {
418  if ((m_faces[2])->GetEid(i) == (m_faces[4])->GetEid(j))
419  {
420  edge = dynamic_pointer_cast<SegGeom>(
421  (m_faces[2])->GetEdge(i));
422  m_edges.push_back(edge);
423  check++;
424  }
425  }
426  }
427 
428  if (check < 1)
429  {
430  std::ostringstream errstrm;
431  errstrm << "Connected faces do not share an edge. Faces ";
432  errstrm << (m_faces[1])->GetGlobalID() << ", "
433  << (m_faces[3])->GetGlobalID();
434  ASSERTL0(false, errstrm.str());
435  }
436  else if (check > 1)
437  {
438  std::ostringstream errstrm;
439  errstrm << "Connected faces share more than one edge. Faces ";
440  errstrm << (m_faces[1])->GetGlobalID() << ", "
441  << (m_faces[3])->GetGlobalID();
442  ASSERTL0(false, errstrm.str());
443  }
444 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
int GetEid(int i) const
Get the ID of edge i of this object.
Definition: Geometry.cpp:145
int GetGlobalID(void) const
Get the ID of this object.
Definition: Geometry.h:314
Geometry1DSharedPtr GetEdge(int i) const
Returns edge i of this object.
Definition: Geometry.h:351
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:62

◆ SetUpLocalVertices()

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

Definition at line 446 of file PrismGeom.cpp.

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

Referenced by PrismGeom().

447 {
448 
449  // Set up the first 2 vertices (i.e. vertex 0,1)
450  if ((m_edges[0]->GetVid(0) == m_edges[1]->GetVid(0)) ||
451  (m_edges[0]->GetVid(0) == m_edges[1]->GetVid(1)))
452  {
453  m_verts.push_back(m_edges[0]->GetVertex(1));
454  m_verts.push_back(m_edges[0]->GetVertex(0));
455  }
456  else if ((m_edges[0]->GetVid(1) == m_edges[1]->GetVid(0)) ||
457  (m_edges[0]->GetVid(1) == m_edges[1]->GetVid(1)))
458  {
459  m_verts.push_back(m_edges[0]->GetVertex(0));
460  m_verts.push_back(m_edges[0]->GetVertex(1));
461  }
462  else
463  {
464  std::ostringstream errstrm;
465  errstrm << "Connected edges do not share a vertex. Edges ";
466  errstrm << m_edges[0]->GetGlobalID() << ", "
467  << m_edges[1]->GetGlobalID();
468  ASSERTL0(false, errstrm.str());
469  }
470 
471  // set up the other bottom vertices (i.e. vertex 2,3)
472  for (int i = 1; i < 3; i++)
473  {
474  if (m_edges[i]->GetVid(0) == m_verts[i]->GetGlobalID())
475  {
476  m_verts.push_back(m_edges[i]->GetVertex(1));
477  }
478  else if (m_edges[i]->GetVid(1) == m_verts[i]->GetGlobalID())
479  {
480  m_verts.push_back(m_edges[i]->GetVertex(0));
481  }
482  else
483  {
484  std::ostringstream errstrm;
485  errstrm << "Connected edges do not share a vertex. Edges ";
486  errstrm << m_edges[i]->GetGlobalID() << ", "
487  << m_edges[i - 1]->GetGlobalID();
488  ASSERTL0(false, errstrm.str());
489  }
490  }
491 
492  // set up top vertices
493  // First, set up vertices 4,5
494  if ((m_edges[8]->GetVid(0) == m_edges[4]->GetVid(0)) ||
495  (m_edges[8]->GetVid(0) == m_edges[4]->GetVid(1)))
496  {
497  m_verts.push_back(m_edges[8]->GetVertex(0));
498  m_verts.push_back(m_edges[8]->GetVertex(1));
499  }
500  else if ((m_edges[8]->GetVid(1) == m_edges[4]->GetVid(0)) ||
501  (m_edges[8]->GetVid(1) == m_edges[4]->GetVid(1)))
502  {
503  m_verts.push_back(m_edges[8]->GetVertex(1));
504  m_verts.push_back(m_edges[8]->GetVertex(0));
505  }
506  else
507  {
508  std::ostringstream errstrm;
509  errstrm << "Connected edges do not share a vertex. Edges ";
510  errstrm << m_edges[8]->GetGlobalID();
511  ASSERTL0(false, errstrm.str());
512  }
513 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
int GetGlobalID(void) const
Get the ID of this object.
Definition: Geometry.h:314
int GetVid(int i) const
Get the ID of vertex i of this object.
Definition: Geometry.cpp:137
PointGeomSharedPtr GetVertex(int i) const
Returns vertex i of this object.
Definition: Geometry.h:343

◆ SetUpXmap()

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 839 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 v_Setup().

840 {
841  vector<int> tmp;
842  int order0, order1;
843 
844  if (m_forient[0] < 9)
845  {
846  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
847  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
848  order0 = *max_element(tmp.begin(), tmp.end());
849  }
850  else
851  {
852  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
853  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(3));
854  order0 = *max_element(tmp.begin(), tmp.end());
855  }
856 
857  if (m_forient[0] < 9)
858  {
859  tmp.clear();
860  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
861  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(3));
862  tmp.push_back(m_faces[2]->GetXmap()->GetEdgeNcoeffs(2));
863  order1 = *max_element(tmp.begin(), tmp.end());
864  }
865  else
866  {
867  tmp.clear();
868  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
869  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
870  tmp.push_back(m_faces[2]->GetXmap()->GetEdgeNcoeffs(2));
871  order1 = *max_element(tmp.begin(), tmp.end());
872  }
873 
874  tmp.clear();
875  tmp.push_back(order0);
876  tmp.push_back(order1);
877  tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(1));
878  tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(2));
879  tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(1));
880  tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(2));
881  int order2 = *max_element(tmp.begin(), tmp.end());
882 
883  const LibUtilities::BasisKey A(
885  order0,
886  LibUtilities::PointsKey(order0+1, LibUtilities::eGaussLobattoLegendre));
887  const LibUtilities::BasisKey B(
889  order1,
890  LibUtilities::PointsKey(order1+1, LibUtilities::eGaussLobattoLegendre));
891  const LibUtilities::BasisKey C(
893  order2,
894  LibUtilities::PointsKey(order2,
896 
898 }
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
Definition: Geometry.h:189
StdRegions::StdExpansionSharedPtr GetXmap() const
Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standar...
Definition: Geometry.h:421
Principle Modified Functions .
Definition: BasisType.h:48
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:87
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
Principle Modified Functions .
Definition: BasisType.h:49
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51

◆ v_ContainsPoint()

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

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

The parameter tol which is by default 0, can be used to expand the coordinate range of the standard element from \([-1,1]^d\) to \([-1-\epsilon,1+\epsilon\) to handle challenging edge cases. The function also returns the local coordinates corresponding to gloCoord that can be used to speed up subsequent searches.

Parameters
gloCoordThe coordinate \( (x,y,z) \).
locCoordOn exit, this is the local coordinate \(\vec{\xi}\) such that \(\chi(\vec{\xi}) = \vec{x}\).
tolThe tolerance used to dictate the bounding box of the standard coordinates \(\vec{\xi}\).
residOn exit, returns the minimum distance between gloCoord and the quadrature points inside the element.
Returns
true if the coordinate gloCoord is contained in the element; false otherwise.
See also
Geometry::GetLocCoords.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 101 of file PrismGeom.cpp.

References Nektar::SpatialDomains::Geometry::ClampLocCoords(), Nektar::SpatialDomains::eRegular, Nektar::SpatialDomains::Geometry::GetLocCoords(), Nektar::SpatialDomains::Geometry::GetMetricInfo(), and Nektar::SpatialDomains::Geometry::MinMaxCheck().

105 {
106  //Rough check if within twice min/max point
107  if (GetMetricInfo()->GetGtype() != eRegular)
108  {
109  if (!MinMaxCheck(gloCoord))
110  {
111  return false;
112  }
113  }
114 
115  // Convert to the local Cartesian coordinates.
116  resid = GetLocCoords(gloCoord, locCoord);
117 
118  // Check local coordinate is within std region bounds.
119  if (locCoord[0] >= -(1 + tol) && locCoord[1] >= -(1 + tol) &&
120  locCoord[2] >= -(1 + tol) && locCoord[1] <= (1 + tol) &&
121  locCoord[0] + locCoord[2] <= tol)
122  {
123  return true;
124  }
125 
126  //Clamp local coords
127  ClampLocCoords(locCoord, tol);
128 
129  return false;
130 }
bool MinMaxCheck(const Array< OneD, const NekDouble > &gloCoord)
Check if given global coord is within twice the min/max distance of the element.
Definition: Geometry.cpp:435
void ClampLocCoords(Array< OneD, NekDouble > &locCoord, NekDouble tol)
Clamp local coords to be within standard regions [-1, 1]^dim.
Definition: Geometry.cpp:478
Geometry is straight-sided with constant geometric factors.
NekDouble GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this ge...
Definition: Geometry.h:534
GeomFactorsSharedPtr GetMetricInfo()
Get the geometric factors for this object.
Definition: Geometry.h:298

◆ v_GenGeomFactors()

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

Implements Nektar::SpatialDomains::Geometry.

Definition at line 132 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::Geometry::m_setupState, Nektar::SpatialDomains::Geometry3D::m_verts, Nektar::SpatialDomains::Geometry::m_xmap, Nektar::SpatialDomains::Geometry3D::v_FillGeom(), and v_Setup().

133 {
134  if(!m_setupState)
135  {
137  }
138 
140  {
141  int i, f;
142  GeomType Gtype = eRegular;
143 
144  v_FillGeom();
145 
146  // check to see if expansions are linear
147  for (i = 0; i < m_coordim; ++i)
148  {
149  if (m_xmap->GetBasisNumModes(0) != 2 ||
150  m_xmap->GetBasisNumModes(1) != 2 ||
151  m_xmap->GetBasisNumModes(2) != 2)
152  {
153  Gtype = eDeformed;
154  }
155  }
156 
157  // check to see if all quadrilateral faces are parallelograms
158  if (Gtype == eRegular)
159  {
160  // Vertex ids of quad faces
161  const unsigned int faceVerts[3][4] = {
162  {0, 1, 2, 3}, {1, 2, 5, 4}, {0, 3, 5, 4}};
163 
164  for (f = 0; f < 3; f++)
165  {
166  // Ensure each face is a parallelogram? Check this.
167  for (i = 0; i < m_coordim; i++)
168  {
169  if (fabs((*m_verts[faceVerts[f][0]])(i) -
170  (*m_verts[faceVerts[f][1]])(i) +
171  (*m_verts[faceVerts[f][2]])(i) -
172  (*m_verts[faceVerts[f][3]])(i)) >
174  {
175  Gtype = eDeformed;
176  break;
177  }
178  }
179 
180  if (Gtype == eDeformed)
181  {
182  break;
183  }
184  }
185  }
186 
188  Gtype, m_coordim, m_xmap, m_coeffs);
190  }
191 }
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
Definition: Geometry.h:189
GeomFactorsSharedPtr m_geomFactors
Geometric factors.
Definition: Geometry.h:185
GeomState m_geomFactorsState
State of the geometric factors.
Definition: Geometry.h:187
static const NekDouble kNekZeroTol
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
bool m_setupState
Wether or not the setup routines have been run.
Definition: Geometry.h:193
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
Definition: Geometry3D.cpp:235
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Array containing expansion coefficients of m_xmap.
Definition: Geometry.h:201
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 of this geometry object.
Definition: Geometry.h:183

◆ v_GetDir()

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

Implements Nektar::SpatialDomains::Geometry3D.

Definition at line 85 of file PrismGeom.cpp.

86 {
87  if (faceidx == 0)
88  {
89  return facedir;
90  }
91  else if (faceidx == 1 || faceidx == 3)
92  {
93  return 2 * facedir;
94  }
95  else
96  {
97  return 1 + facedir;
98  }
99 }

◆ v_GetEdgeFaceMap()

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

Returns the standard element edge IDs that are connected to a given face.

For example, on a prism, edge 0 is connnected to faces 0 and 1; GetEdgeFaceMap(0,j) would therefore return the values 0 and 1 respectively. We assume that j runs between 0 and 1 inclusive, since every face is connected to precisely two faces for all 3D elements.

This function is used in the construction of the low-energy preconditioner.

Parameters
iThe edge to query connectivity for.
jThe local face index between 0 and 1 connected to this element.
See also
MultiRegions::PreconditionerLowEnergy

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 294 of file PrismGeom.cpp.

References EdgeFaceConnectivity.

295 {
296  return EdgeFaceConnectivity[i][j];
297 }
static const unsigned int EdgeFaceConnectivity[9][2]
Definition: PrismGeom.h:85

◆ v_GetLocCoords()

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

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

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

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

Parameters
coordsInput Cartesian global coordinates
LcoordsCorresponding local coordinates
Returns
Distance between obtained coordinates and provided ones.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 193 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(), Vmath::Sadd(), Nektar::SpatialDomains::PointGeom::Sub(), Nektar::SpatialDomains::Geometry3D::v_FillGeom(), Vmath::Vmul(), and Vmath::Vvtvp().

195 {
196  NekDouble ptdist = 1e6;
197 
198  // calculate local coordinate for coord
199  if (GetMetricInfo()->GetGtype() == eRegular)
200  {
201  // Point inside tetrahedron
202  PointGeom r(m_coordim, 0, coords[0], coords[1], coords[2]);
203 
204  // Edges
205  PointGeom er0, e10, e30, e40;
206  er0.Sub(r, *m_verts[0]);
207  e10.Sub(*m_verts[1], *m_verts[0]);
208  e30.Sub(*m_verts[3], *m_verts[0]);
209  e40.Sub(*m_verts[4], *m_verts[0]);
210 
211  // Cross products (Normal times area)
212  PointGeom cp1030, cp3040, cp4010;
213  cp1030.Mult(e10, e30);
214  cp3040.Mult(e30, e40);
215  cp4010.Mult(e40, e10);
216 
217  // Barycentric coordinates (relative volume)
218  NekDouble V = e40.dot(cp1030); // Prism Volume = {(e40)dot(e10)x(e30)}/2
219  NekDouble beta =
220  er0.dot(cp3040) / (2.0 * V); // volume1 = {(er0)dot(e30)x(e40)}/4
221  NekDouble gamma =
222  er0.dot(cp4010) / (3.0 * V); // volume2 = {(er0)dot(e40)x(e10)}/6
223  NekDouble delta =
224  er0.dot(cp1030) / (2.0 * V); // volume3 = {(er0)dot(e10)x(e30)}/4
225 
226  // Make Prism bigger
227  Lcoords[0] = 2.0 * beta - 1.0;
228  Lcoords[1] = 2.0 * gamma - 1.0;
229  Lcoords[2] = 2.0 * delta - 1.0;
230 
231  // Set ptdist to distance to nearest vertex
232  for (int i = 0; i < 6; ++i)
233  {
234  ptdist = min(ptdist, r.dist(*m_verts[i]));
235  }
236  }
237  else
238  {
239  v_FillGeom();
240 
241  // Determine nearest point of coords to values in m_xmap
242  int npts = m_xmap->GetTotPoints();
243  Array<OneD, NekDouble> ptsx(npts), ptsy(npts), ptsz(npts);
244  Array<OneD, NekDouble> tmp1(npts), tmp2(npts);
245 
246  m_xmap->BwdTrans(m_coeffs[0], ptsx);
247  m_xmap->BwdTrans(m_coeffs[1], ptsy);
248  m_xmap->BwdTrans(m_coeffs[2], ptsz);
249 
250  const Array<OneD, const NekDouble> za = m_xmap->GetPoints(0);
251  const Array<OneD, const NekDouble> zb = m_xmap->GetPoints(1);
252  const Array<OneD, const NekDouble> zc = m_xmap->GetPoints(2);
253 
254  // guess the first local coords based on nearest point
255  Vmath::Sadd(npts, -coords[0], ptsx, 1, tmp1, 1);
256  Vmath::Vmul(npts, tmp1, 1, tmp1, 1, tmp1, 1);
257  Vmath::Sadd(npts, -coords[1], ptsy, 1, tmp2, 1);
258  Vmath::Vvtvp(npts, tmp2, 1, tmp2, 1, tmp1, 1, tmp1, 1);
259  Vmath::Sadd(npts, -coords[2], ptsz, 1, tmp2, 1);
260  Vmath::Vvtvp(npts, tmp2, 1, tmp2, 1, tmp1, 1, tmp1, 1);
261 
262  int min_i = Vmath::Imin(npts, tmp1, 1);
263 
264  // distance from coordinate to nearest point for return value.
265  ptdist = sqrt(tmp1[min_i]);
266 
267  // Get collapsed coordinate
268  int qa = za.num_elements(), qb = zb.num_elements();
269  Lcoords[2] = zc[min_i / (qa * qb)];
270  min_i = min_i % (qa * qb);
271  Lcoords[1] = zb[min_i / qa];
272  Lcoords[0] = za[min_i % qa];
273 
274  // recover cartesian coordinate from collapsed coordinate.
275  Lcoords[0] = (1.0 + Lcoords[0]) * (1.0 - Lcoords[2]) / 2 - 1.0;
276 
277  // Perform newton iteration to find local coordinates
278  NekDouble resid = 0.0;
279  NewtonIterationForLocCoord(coords, ptsx, ptsy, ptsz, Lcoords, resid);
280  }
281  return ptdist;
282 }
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
Definition: Geometry.h:189
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:850
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:445
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:81
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:318
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
Definition: Geometry3D.cpp:235
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Array containing expansion coefficients of m_xmap.
Definition: Geometry.h:201
Geometry is straight-sided with constant geometric factors.
GeomFactorsSharedPtr GetMetricInfo()
Get the geometric factors for this object.
Definition: Geometry.h:298
int m_coordim
Coordinate dimension of this geometry object.
Definition: Geometry.h:183
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:186

◆ v_GetVertexEdgeMap()

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

Returns the standard element edge IDs that are connected to a given vertex.

For example, on a prism, vertex 0 is connnected to edges 0, 3, and 4; GetVertexEdgeMap(0,j) would therefore return the values 0, 1 and 4 respectively. We assume that j runs between 0 and 2 inclusive, which is true for every 3D element asides from the pyramid.

This function is used in the construction of the low-energy preconditioner.

Parameters
iThe vertex to query connectivity for.
jThe local edge index between 0 and 2 connected to this element.
Todo:
Expand to work with pyramid elements.
See also
MultiRegions::PreconditionerLowEnergy

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 284 of file PrismGeom.cpp.

References VertexEdgeConnectivity.

285 {
286  return VertexEdgeConnectivity[i][j];
287 }
static const unsigned int VertexEdgeConnectivity[6][3]
Definition: PrismGeom.h:83

◆ v_GetVertexFaceMap()

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

Returns the standard element face IDs that are connected to a given vertex.

For example, on a hexahedron, vertex 0 is connnected to faces 0, 1, and 4; GetVertexFaceMap(0,j) would therefore return the values 0, 1 and 4 respectively. We assume that j runs between 0 and 2 inclusive, which is true for every 3D element asides from the pyramid.

This is used in the construction of the low-energy preconditioner.

Parameters
iThe vertex to query connectivity for.
jThe local face index between 0 and 2 connected to this element.
Todo:
Expand to work with pyramid elements.
See also
MultiRegions::PreconditionerLowEnergy

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 289 of file PrismGeom.cpp.

References VertexFaceConnectivity.

290 {
291  return VertexFaceConnectivity[i][j];
292 }
static const unsigned int VertexFaceConnectivity[6][3]
Definition: PrismGeom.h:84

◆ v_Reset()

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

Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated GeomFactors.

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 811 of file PrismGeom.cpp.

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

812 {
813  Geometry::v_Reset(curvedEdges, curvedFaces);
814 
815  for (int i = 0; i < 5; ++i)
816  {
817  m_faces[i]->Reset(curvedEdges, curvedFaces);
818  }
819 }
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated Geo...
Definition: Geometry.cpp:335

◆ v_Setup()

void Nektar::SpatialDomains::PrismGeom::v_Setup ( )
protectedvirtual

Reimplemented from Nektar::SpatialDomains::Geometry.

Definition at line 821 of file PrismGeom.cpp.

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

Referenced by v_GenGeomFactors().

822 {
823  if(!m_setupState)
824  {
825  for (int i = 0; i < 5; ++i)
826  {
827  m_faces[i]->Setup();
828  }
829  SetUpXmap();
830  SetUpCoeffs(m_xmap->GetNcoeffs());
831  m_setupState = true;
832  }
833 }
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
Definition: Geometry.h:189
void SetUpXmap()
Set up the m_xmap object by determining the order of each direction from derived faces.
Definition: PrismGeom.cpp:839
bool m_setupState
Wether or not the setup routines have been run.
Definition: Geometry.h:193
void SetUpCoeffs(const int nCoeffs)
Initialise the Geometry::m_coeffs array.
Definition: Geometry.h:643

Member Data Documentation

◆ EdgeFaceConnectivity

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 85 of file PrismGeom.h.

Referenced by v_GetEdgeFaceMap().

◆ kNedges

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

Definition at line 55 of file PrismGeom.h.

Referenced by PrismGeom(), and SetUpEdgeOrientation().

◆ kNfaces

const int Nektar::SpatialDomains::PrismGeom::kNfaces = kNqfaces + kNtfaces
static

Definition at line 58 of file PrismGeom.h.

Referenced by PrismGeom(), and SetUpFaceOrientation().

◆ kNqfaces

const int Nektar::SpatialDomains::PrismGeom::kNqfaces = 3
static

Definition at line 56 of file PrismGeom.h.

Referenced by SetUpFaceOrientation().

◆ kNtfaces

const int Nektar::SpatialDomains::PrismGeom::kNtfaces = 2
static

Definition at line 57 of file PrismGeom.h.

Referenced by SetUpFaceOrientation().

◆ kNverts

const int Nektar::SpatialDomains::PrismGeom::kNverts = 6
static

Definition at line 54 of file PrismGeom.h.

◆ VertexEdgeConnectivity

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 83 of file PrismGeom.h.

Referenced by v_GetVertexEdgeMap().

◆ VertexFaceConnectivity

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 84 of file PrismGeom.h.

Referenced by v_GetVertexFaceMap().

◆ XMLElementType

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

Definition at line 59 of file PrismGeom.h.