46 namespace SpatialDomains
48 const unsigned int HexGeom::VertexEdgeConnectivity[8][3] = {
49 {0,3,4},{0,1,5},{1,2,6},{2,3,7},
50 {4,8,11},{5,8,9},{6,9,10},{7,10,11}};
51 const unsigned int HexGeom::VertexFaceConnectivity[8][3] = {
52 {0,1,4},{0,1,2},{0,2,3},{0,3,4},
53 {1,4,5},{1,2,5},{2,3,5},{3,4,5}};
54 const unsigned int HexGeom::EdgeFaceConnectivity[12][2] = {
55 {0,1},{0,2},{0,3},{0,4},{1,4},{1,2},{2,3},{3,4},
56 {1,5},{2,5},{3,5},{4,5}};
64 Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
100 if (
m_xmap->GetBasisNumModes(0) != 2 ||
101 m_xmap->GetBasisNumModes(1) != 2 ||
102 m_xmap->GetBasisNumModes(2) != 2 )
124 if( fabs( (*
m_verts[ faceVerts[f][0] ])(i) -
125 (*
m_verts[ faceVerts[f][1] ])(i) +
126 (*
m_verts[ faceVerts[f][2] ])(i) -
127 (*
m_verts[ faceVerts[f][3] ])(i) )
173 nq0 =
m_xmap->GetNumPoints(0);
174 nq1 =
m_xmap->GetNumPoints(1);
175 nq2 =
m_xmap->GetNumPoints(2);
180 len0 += (pts[nq0-1]-pts[0])*(pts[nq0-1]-pts[0]);
181 xi0 += (coords[i] -pts[0])*(pts[nq0-1]-pts[0]);
184 len1 += (pts[nq0*(nq1-1)]-pts[0])*(pts[nq0*(nq1-1)]-pts[0]);
185 xi1 += (coords[i] -pts[0])*(pts[nq0*(nq1-1)]-pts[0]);
188 len2 += (pts[nq0*nq1*(nq2-1)]-pts[0])*(pts[nq0*nq1*(nq2-1)]-pts[0]);
189 xi2 += (coords[i] -pts[0])*(pts[nq0*nq1*(nq2-1)]-pts[0]);
192 Lcoords[0] = 2*xi0/len0-1.0;
193 Lcoords[1] = 2*xi1/len1-1.0;
194 Lcoords[2] = 2*xi2/len2-1.0;
198 PointGeom r(m_coordim, 0, coords[0], coords[1], coords[2]);
199 for(
int i = 0; i < 8; ++i)
230 ptdist = sqrt(tmp1[min_i]);
233 int qa = za.num_elements(), qb = zb.num_elements();
234 Lcoords[2] = zc[min_i/(qa*qb)];
235 min_i = min_i%(qa*qb);
236 Lcoords[1] = zb[min_i/qa];
237 Lcoords[0] = za[min_i%qa];
274 ASSERTL1(gloCoord.num_elements() == 3,
275 "Three dimensional geometry expects three coordinates.");
291 for(i = 0; i < 3; ++i)
295 mincoord[i] =
Vmath::Vmin(pts.num_elements(),pts,1);
296 maxcoord[i] =
Vmath::Vmax(pts.num_elements(),pts,1);
298 diff = max(maxcoord[i] - mincoord[i],diff);
301 for(i = 0; i < 3; ++i)
303 if((gloCoord[i] < mincoord[i] - 0.2*diff)||
304 (gloCoord[i] > maxcoord[i] + 0.2*diff))
313 if (locCoord[0] >= -(1+tol) && locCoord[0] <= 1+tol
314 && locCoord[1] >= -(1+tol) && locCoord[1] <= 1+tol
315 && locCoord[2] >= -(1+tol) && locCoord[2] <= 1+tol)
324 for(
int i = 0; i < 3; ++i)
326 if(locCoord[i] <-(1+tol))
328 locCoord[i] = -(1+tol);
331 if(locCoord[i] > (1+tol))
372 if (faceidx == 0 || faceidx == 5)
376 else if (faceidx == 1 || faceidx == 3)
396 for(f = 1; f < 5 ; f++)
399 for(i = 0; i < 4; i++)
401 for(j = 0; j < 4; j++)
414 std::ostringstream errstrm;
415 errstrm <<
"Connected faces do not share an edge. Faces ";
421 std::ostringstream errstrm;
422 errstrm <<
"Connected faces share more than one edge. Faces ";
430 for(i = 0; i < 4; i++)
432 for(j = 0; j < 4; j++)
444 std::ostringstream errstrm;
445 errstrm <<
"Connected faces do not share an edge. Faces ";
451 std::ostringstream errstrm;
452 errstrm <<
"Connected faces share more than one edge. Faces ";
456 for(f = 1; f < 4 ; f++)
459 for(i = 0; i < 4; i++)
461 for(j = 0; j < 4; j++)
474 std::ostringstream errstrm;
475 errstrm <<
"Connected faces do not share an edge. Faces ";
481 std::ostringstream errstrm;
482 errstrm <<
"Connected faces share more than one edge. Faces ";
489 for(f = 1; f < 5 ; f++)
492 for(i = 0; i < 4; i++)
494 for(j = 0; j < 4; j++)
507 std::ostringstream errstrm;
508 errstrm <<
"Connected faces do not share an edge. Faces ";
514 std::ostringstream errstrm;
515 errstrm <<
"Connected faces share more than one edge. Faces ";
539 std::ostringstream errstrm;
540 errstrm <<
"Connected edges do not share a vertex. Edges ";
547 for(i = 1; i < 3; i++)
559 std::ostringstream errstrm;
560 errstrm <<
"Connected edges do not share a vertex. Edges ";
561 errstrm <<
m_edges[i]->GetEid() <<
", " <<
m_edges[i-1]->GetEid();
582 std::ostringstream errstrm;
583 errstrm <<
"Connected edges do not share a vertex. Edges ";
589 for(i = 9; i < 11; i++)
601 std::ostringstream errstrm;
602 errstrm <<
"Connected edges do not share a vertex. Edges ";
603 errstrm <<
m_edges[i]->GetEid() <<
", " <<
m_edges[i-1]->GetEid();
632 unsigned int baseVertex;
655 unsigned int orientation;
661 elementAaxis_length = 0.0;
662 elementBaxis_length = 0.0;
663 faceAaxis_length = 0.0;
664 faceBaxis_length = 0.0;
669 baseVertex =
m_faces[f]->GetVid(0);
681 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
682 elementBaxis[i] = (*
m_verts[ faceVerts[f][3] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
685 else if( baseVertex ==
m_verts[ faceVerts[f][1] ]->
GetVid() )
689 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
690 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][1] ])[i];
693 else if( baseVertex ==
m_verts[ faceVerts[f][2] ]->
GetVid() )
697 elementAaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][3] ])[i];
698 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][1] ])[i];
701 else if( baseVertex ==
m_verts[ faceVerts[f][3] ]->
GetVid() )
705 elementAaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][3] ])[i];
706 elementBaxis[i] = (*
m_verts[ faceVerts[f][3] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
711 ASSERTL0(
false,
"Could not find matching vertex for the face");
721 elementAaxis_length += pow(elementAaxis[i],2);
722 elementBaxis_length += pow(elementBaxis[i],2);
723 faceAaxis_length += pow(faceAaxis[i],2);
724 faceBaxis_length += pow(faceBaxis[i],2);
727 elementAaxis_length = sqrt(elementAaxis_length);
728 elementBaxis_length = sqrt(elementBaxis_length);
729 faceAaxis_length = sqrt(faceAaxis_length);
730 faceBaxis_length = sqrt(faceBaxis_length);
736 dotproduct1 += elementAaxis[i]*faceAaxis[i];
746 if(dotproduct1 < 0.0)
754 dotproduct2 += elementBaxis[i]*faceBaxis[i];
758 ASSERTL1(fabs(elementBaxis_length*faceBaxis_length - fabs(dotproduct2)) <
760 "These vectors should be parallel");
764 if( dotproduct2 < 0.0 )
779 dotproduct1 += elementAaxis[i]*faceBaxis[i];
783 ASSERTL1(fabs(elementAaxis_length*faceBaxis_length - fabs(dotproduct1)) <
785 "These vectors should be parallel");
789 if(dotproduct1 < 0.0)
798 dotproduct2 += elementBaxis[i]*faceAaxis[i];
802 ASSERTL1(fabs(elementBaxis_length*faceAaxis_length - fabs(dotproduct2)) <
804 "These vectors should be parallel");
806 if( dotproduct2 < 0.0 )
812 orientation = orientation + 5;
824 const unsigned int edgeVerts[
kNedges][2] =
851 ASSERTL0(
false,
"Could not find matching vertex for the edge");
862 for (
int i = 0; i < 6; ++i)
864 m_faces[i]->Reset(curvedEdges, curvedFaces);
879 vector<int> tmp1, tmp2;
911 int order0 = *max_element(tmp1.begin(), tmp1.end());
946 int order1 = *max_element(tmp1.begin(), tmp1.end());
981 int order2 = *max_element(tmp1.begin(), tmp1.end());
static const unsigned int VertexEdgeConnectivity[8][3]
StdRegions::StdExpansionSharedPtr m_xmap
void SetUpXmap()
Set up the m_xmap object by determining the order of each direction from derived faces.
#define ASSERTL0(condition, msg)
std::vector< StdRegions::Orientation > m_eorient
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual void v_GenGeomFactors()
GeomFactorsSharedPtr m_geomFactors
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
static const int kNtfaces
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
int GetEid(int i) const
Return the ID of edge i in this element.
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
Principle Modified Functions .
virtual int v_GetNumFaces() const
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
virtual int v_GetNumEdges() const
StdRegions::StdExpansionSharedPtr GetXmap() const
GeomState m_geomFactorsState
std::vector< StdRegions::Orientation > m_forient
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)
static const NekDouble kNekZeroTol
virtual bool v_ContainsPoint(const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0)
Determines if a point specified in global coordinates is located within this tetrahedral geometry...
void SetUpLocalVertices()
boost::shared_ptr< SegGeom > SegGeomSharedPtr
virtual NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
virtual int v_GetVertexEdgeMap(const int i, const int j) const
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
Defines a specification for a set of points.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
void SetUpEdgeOrientation()
PointGeomSharedPtr GetVertex(int i) const
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
virtual int v_GetEdgeFaceMap(const int i, const int j) const
const Geometry1DSharedPtr GetEdge(int i) const
virtual int v_GetVertexFaceMap(const int i, const int j) const
Array< OneD, Array< OneD, NekDouble > > m_coeffs
static const unsigned int EdgeFaceConnectivity[12][2]
virtual int v_GetNumVerts() const
Geometry is straight-sided with constant geometric factors.
void SetUpFaceOrientation()
virtual int v_GetDir(const int faceidx, const int facedir) const
LibUtilities::ShapeType m_shapeType
static const unsigned int VertexFaceConnectivity[8][3]
NekDouble dist(PointGeom &a)
boost::unordered_map< int, CurveSharedPtr > CurveMap
Geometric information has been generated.
void SetUpCoeffs(const int nCoeffs)
Initialise the m_coeffs array.
static const int kNqfaces
GeomFactorsSharedPtr GetMetricInfo()
GeomType
Indicates the type of element geometry.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Geometry is curved or has non-constant factors.
int m_coordim
coordinate dimension
Describes the specification for a Basis.
1D Gauss-Lobatto-Legendre quadrature points
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.