46 {4, 5, 6, 7}, {1, 3, 6, 7}, {0, 2, 4, 7}, {1, 3, 4, 5}, {0, 2, 5, 6}};
54 :
Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
87 if (
m_xmap->GetBasisNumModes(0) != 2 ||
88 m_xmap->GetBasisNumModes(1) != 2 ||
89 m_xmap->GetBasisNumModes(2) != 2)
99 for (
int i = 0; i < 3; ++i)
142 else if (faceidx == 1 || faceidx == 3)
167 for (f = 1; f < 5; f++)
169 int nEdges =
m_faces[f]->GetNumEdges();
171 for (i = 0; i < 4; i++)
173 for (j = 0; j < nEdges; j++)
177 edge = std::dynamic_pointer_cast<SegGeom>(
187 std::ostringstream errstrm;
188 errstrm <<
"Connected faces do not share an edge. Faces ";
195 std::ostringstream errstrm;
196 errstrm <<
"Connected faces share more than one edge. Faces ";
205 for (i = 0; i < 3; i++)
207 for (j = 0; j < 3; j++)
211 edge = std::dynamic_pointer_cast<SegGeom>(
220 std::ostringstream errstrm;
221 errstrm <<
"Connected faces do not share an edge. Faces ";
228 std::ostringstream errstrm;
229 errstrm <<
"Connected faces share more than one edge. Faces ";
236 for (f = 1; f < 4; f++)
239 for (i = 0; i <
m_faces[f]->GetNumEdges(); i++)
241 for (j = 0; j <
m_faces[f + 1]->GetNumEdges(); j++)
245 edge = std::dynamic_pointer_cast<SegGeom>(
255 std::ostringstream errstrm;
256 errstrm <<
"Connected faces do not share an edge. Faces ";
263 std::ostringstream errstrm;
264 errstrm <<
"Connected faces share more than one edge. Faces ";
289 std::ostringstream errstrm;
290 errstrm <<
"Connected edges do not share a vertex. Edges ";
291 errstrm <<
m_edges[0]->GetGlobalID() <<
", "
297 for (
int i = 1; i < 3; i++)
309 std::ostringstream errstrm;
310 errstrm <<
"Connected edges do not share a vertex. Edges ";
311 errstrm <<
m_edges[i]->GetGlobalID() <<
", "
312 <<
m_edges[i - 1]->GetGlobalID();
328 for (
int i = 5; i < 8; ++i)
340 std::ostringstream errstrm;
341 errstrm <<
"Connected edges do not share a vertex. Edges ";
342 errstrm <<
m_edges[3]->GetGlobalID() <<
", "
353 const unsigned int edgeVerts[
kNedges][2] = {{0, 1}, {1, 2}, {3, 2}, {0, 3},
354 {0, 4}, {1, 4}, {2, 4}, {3, 4}};
370 ASSERTL0(
false,
"Could not find matching vertex for the edge");
396 unsigned int baseVertex;
407 const unsigned int faceVerts[
kNfaces][4] = {
417 unsigned int orientation;
423 elementAaxis_length = 0.0;
424 elementBaxis_length = 0.0;
425 faceAaxis_length = 0.0;
426 faceBaxis_length = 0.0;
431 baseVertex =
m_faces[f]->GetVid(0);
448 elementAaxis[i] = (*
m_verts[faceVerts[f][1]])[i] -
449 (*
m_verts[faceVerts[f][0]])[i];
450 elementBaxis[i] = (*
m_verts[faceVerts[f][2]])[i] -
451 (*
m_verts[faceVerts[f][0]])[i];
458 elementAaxis[i] = (*
m_verts[faceVerts[f][1]])[i] -
459 (*
m_verts[faceVerts[f][0]])[i];
460 elementBaxis[i] = (*
m_verts[faceVerts[f][2]])[i] -
461 (*
m_verts[faceVerts[f][1]])[i];
468 elementAaxis[i] = (*
m_verts[faceVerts[f][1]])[i] -
469 (*
m_verts[faceVerts[f][2]])[i];
470 elementBaxis[i] = (*
m_verts[faceVerts[f][2]])[i] -
471 (*
m_verts[faceVerts[f][0]])[i];
476 ASSERTL0(
false,
"Could not find matching vertex for the face");
485 elementAaxis[i] = (*
m_verts[faceVerts[f][1]])[i] -
486 (*
m_verts[faceVerts[f][0]])[i];
487 elementBaxis[i] = (*
m_verts[faceVerts[f][3]])[i] -
488 (*
m_verts[faceVerts[f][0]])[i];
495 elementAaxis[i] = (*
m_verts[faceVerts[f][1]])[i] -
496 (*
m_verts[faceVerts[f][0]])[i];
497 elementBaxis[i] = (*
m_verts[faceVerts[f][2]])[i] -
498 (*
m_verts[faceVerts[f][1]])[i];
505 elementAaxis[i] = (*
m_verts[faceVerts[f][2]])[i] -
506 (*
m_verts[faceVerts[f][3]])[i];
507 elementBaxis[i] = (*
m_verts[faceVerts[f][2]])[i] -
508 (*
m_verts[faceVerts[f][1]])[i];
515 elementAaxis[i] = (*
m_verts[faceVerts[f][2]])[i] -
516 (*
m_verts[faceVerts[f][3]])[i];
517 elementBaxis[i] = (*
m_verts[faceVerts[f][3]])[i] -
518 (*
m_verts[faceVerts[f][0]])[i];
523 ASSERTL0(
false,
"Could not find matching vertex for the face");
531 int v =
m_faces[f]->GetNumVerts() - 1;
537 elementAaxis_length += pow(elementAaxis[i], 2);
538 elementBaxis_length += pow(elementBaxis[i], 2);
539 faceAaxis_length += pow(faceAaxis[i], 2);
540 faceBaxis_length += pow(faceBaxis[i], 2);
543 elementAaxis_length =
sqrt(elementAaxis_length);
544 elementBaxis_length =
sqrt(elementBaxis_length);
545 faceAaxis_length =
sqrt(faceAaxis_length);
546 faceBaxis_length =
sqrt(faceBaxis_length);
552 dotproduct1 += elementAaxis[i] * faceAaxis[i];
556 fabs(dotproduct1) / elementAaxis_length / faceAaxis_length;
566 if (dotproduct1 < 0.0)
574 dotproduct2 += elementBaxis[i] * faceBaxis[i];
577 norm = fabs(dotproduct2) / elementBaxis_length / faceBaxis_length;
581 "These vectors should be parallel");
585 if (dotproduct2 < 0.0)
600 dotproduct1 += elementAaxis[i] * faceBaxis[i];
603 norm = fabs(dotproduct1) / elementAaxis_length / faceBaxis_length;
605 "These vectors should be parallel");
609 if (dotproduct1 < 0.0)
618 dotproduct2 += elementBaxis[i] * faceAaxis[i];
621 norm = fabs(dotproduct2) / elementBaxis_length / faceAaxis_length;
625 "These vectors should be parallel");
627 if (dotproduct2 < 0.0)
633 orientation = orientation + 5;
639 "Orientation of triangular face (id = " +
641 ") is inconsistent with face " + std::to_string(f) +
642 " of pyramid element (id = " + std::to_string(
m_globalID) +
643 ") since Dir2 is aligned with Dir1. Mesh setup "
644 "needs investigation");
656 for (
int i = 0; i < 5; ++i)
658 m_faces[i]->Reset(curvedEdges, curvedFaces);
669 for (
int i = 0; i < 5; ++i)
685 std::vector<int> tmp;
692 order0 = *max_element(tmp.begin(), tmp.end());
698 order0 = *max_element(tmp.begin(), tmp.end());
707 order1 = *max_element(tmp.begin(), tmp.end());
715 order1 = *max_element(tmp.begin(), tmp.end());
719 tmp.push_back(order0);
720 tmp.push_back(order1);
725 int order2 = *max_element(tmp.begin(), tmp.end());
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Describes the specification for a Basis.
Defines a specification for a set of points.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::vector< StdRegions::Orientation > m_forient
void v_CalculateInverseIsoParam() override
void v_FillGeom() override
Put all quadrature information into face/edge structure and backward transform.
std::vector< StdRegions::Orientation > m_eorient
bool m_setupState
Wether or not the setup routines have been run.
PointGeomSharedPtr GetVertex(int i) const
Returns vertex i of this object.
void SetUpCoeffs(const int nCoeffs)
Initialise the Geometry::m_coeffs array.
Geometry1DSharedPtr GetEdge(int i) const
Returns edge i of this object.
int GetVid(int i) const
Get the ID of vertex i of this object.
Array< OneD, Array< OneD, NekDouble > > m_isoParameter
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated Geo...
int GetGlobalID(void) const
Get the ID of this object.
LibUtilities::ShapeType m_shapeType
Type of shape.
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Array containing expansion coefficients of m_xmap.
GeomState m_geomFactorsState
State of the geometric factors.
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
StdRegions::StdExpansionSharedPtr GetXmap() const
Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standar...
GeomFactorsSharedPtr m_geomFactors
Geometric factors.
int m_coordim
Coordinate dimension of this geometry object.
int GetEid(int i) const
Get the ID of edge i of this object.
void v_GenGeomFactors() override
void SetUpFaceOrientation()
int v_GetDir(const int faceidx, const int facedir) const override
Returns the element coordinate direction corresponding to a given face coordinate direction.
void SetUpEdgeOrientation()
void SetUpLocalVertices()
void SetUpXmap()
Set up the m_xmap object by determining the order of each direction from derived faces.
void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces) override
Reset this geometry object: unset the current state, zero Geometry::m_coeffs and remove allocated Geo...
static const int kNtfaces
static const int kNqfaces
static const unsigned int EdgeNormalToFaceVert[5][4]
int v_GetEdgeNormalToFaceVert(const int i, const int j) const override
Returns the standard lement edge IDs that are normal to a given face vertex.
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
@ eModifiedPyr_C
Principle Modified Functions.
@ eModified_A
Principle Modified Functions .
static const NekDouble kNekZeroTol
std::unordered_map< int, CurveSharedPtr > CurveMap
GeomType
Indicates the type of element geometry.
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< SegGeom > SegGeomSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
@ ePtsFilled
Geometric information has been generated.
@ eDir1FwdDir2_Dir2FwdDir1
scalarT< T > sqrt(scalarT< T > in)