46 namespace SpatialDomains
48 const unsigned int TetGeom::VertexEdgeConnectivity[4][3] = {
49 {0,2,3},{0,1,4},{1,2,5},{3,4,5}};
50 const unsigned int TetGeom::VertexFaceConnectivity[4][3] = {
51 {0,1,3},{0,1,2},{0,2,3},{1,2,3}};
52 const unsigned int TetGeom::EdgeFaceConnectivity [6][2] = {
53 {0,1},{0,2},{0,3},{1,3},{1,2},{2,3}};
61 Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
114 ASSERTL1(gloCoord.num_elements() == 3,
115 "Three dimensional geometry expects three coordinates.");
131 for(i = 0; i < 3; ++i)
135 mincoord[i] =
Vmath::Vmin(pts.num_elements(),pts,1);
136 maxcoord[i] =
Vmath::Vmax(pts.num_elements(),pts,1);
138 diff = max(maxcoord[i] - mincoord[i],diff);
141 for(i = 0; i < 3; ++i)
143 if((gloCoord[i] < mincoord[i] - 0.2*diff)||
144 (gloCoord[i] > maxcoord[i] + 0.2*diff))
155 if (locCoord[0] >= -(1+tol) && locCoord[1] >= -(1+tol) &&
156 locCoord[2] >= -(1+tol) &&
157 locCoord[0] + locCoord[1] + locCoord[2] <= -1+tol)
166 for(
int i = 0; i < 3; ++i)
168 if(locCoord[i] <-(1+tol))
170 locCoord[i] = -(1+tol);
173 if(locCoord[i] > (1+tol))
206 cp1020.
Mult(e10,e20);
207 cp2030.
Mult(e20,e30);
208 cp3010.
Mult(e30,e10);
218 Lcoords[0] = 2.0*beta - 1.0;
219 Lcoords[1] = 2.0*gamma - 1.0;
220 Lcoords[2] = 2.0*delta - 1.0;
223 for(
int i = 0; i < 4; ++i)
256 ptdist = sqrt(tmp1[min_i]);
259 int qa = za.num_elements(), qb = zb.num_elements();
260 Lcoords[2] = zc[min_i/(qa*qb)];
261 min_i = min_i%(qa*qb);
262 Lcoords[1] = zb[min_i/qa];
263 Lcoords[0] = za[min_i%qa];
266 Lcoords[1] = (1.0+Lcoords[0])*(1.0-Lcoords[2])/2 -1.0;
267 Lcoords[0] = (1.0+Lcoords[0])*(-Lcoords[1]-Lcoords[2])/2 -1.0;
297 else if (faceidx == 1)
335 std::ostringstream errstrm;
336 errstrm <<
"Local edge 0 (eid=" <<
m_faces[0]->GetEid(0);
337 errstrm <<
") on face " <<
m_faces[0]->GetFid();
338 errstrm <<
" must be the same as local edge 0 (eid="<<
m_faces[1]->GetEid(0);
339 errstrm <<
") on face " <<
m_faces[1]->GetFid();
344 for(faceConnected = 1; faceConnected < 4 ; faceConnected++)
347 for(i = 0; i < 3; i++)
359 std::ostringstream errstrm;
360 errstrm <<
"Face 0 does not share an edge with first edge of adjacent face. Faces ";
366 std::ostringstream errstrm;
367 errstrm <<
"Connected faces share more than one edge. Faces ";
376 for(i = 0; i < 3; i++)
378 for(j = 0; j < 3; j++)
390 std::ostringstream errstrm;
391 errstrm <<
"Connected faces do not share an edge. Faces ";
397 std::ostringstream errstrm;
398 errstrm <<
"Connected faces share more than one edge. Faces ";
403 for(faceConnected = 1; faceConnected < 3 ; faceConnected++)
406 for(i = 0; i < 3; i++)
408 for(j = 0; j < 3; j++)
421 std::ostringstream errstrm;
422 errstrm <<
"Connected faces do not share an edge. Faces ";
428 std::ostringstream errstrm;
429 errstrm <<
"Connected faces share more than one edge. Faces ";
453 std::ostringstream errstrm;
454 errstrm <<
"Connected edges do not share a vertex. Edges ";
460 for(
int i = 1; i < 2; i++)
472 std::ostringstream errstrm;
473 errstrm <<
"Connected edges do not share a vertex. Edges ";
474 errstrm <<
m_edges[i]->GetEid() <<
", " <<
m_edges[i-1]->GetEid();
490 for (
int i = 4; i < 6; ++i)
501 std::ostringstream errstrm;
502 errstrm <<
"Connected edges do not share a vertex. Edges ";
513 const unsigned int edgeVerts[
kNedges][2] =
534 ASSERTL0(
false,
"Could not find matching vertex for the edge");
563 unsigned int baseVertex;
584 unsigned int orientation;
590 elementAaxis_length = 0.0;
591 elementBaxis_length = 0.0;
592 faceAaxis_length = 0.0;
593 faceBaxis_length = 0.0;
598 baseVertex =
m_faces[f]->GetVid(0);
612 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
613 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
616 else if( baseVertex ==
m_verts[ faceVerts[f][1] ]->
GetVid() )
620 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
621 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][1] ])[i];
624 else if( baseVertex ==
m_verts[ faceVerts[f][2] ]->
GetVid() )
628 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][2] ])[i];
629 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
634 ASSERTL0(
false,
"Could not find matching vertex for the face");
644 elementAaxis_length += pow(elementAaxis[i],2);
645 elementBaxis_length += pow(elementBaxis[i],2);
646 faceAaxis_length += pow(faceAaxis[i],2);
647 faceBaxis_length += pow(faceBaxis[i],2);
650 elementAaxis_length = sqrt(elementAaxis_length);
651 elementBaxis_length = sqrt(elementBaxis_length);
652 faceAaxis_length = sqrt(faceAaxis_length);
653 faceBaxis_length = sqrt(faceBaxis_length);
659 dotproduct1 += elementAaxis[i]*faceAaxis[i];
669 if(dotproduct1 < 0.0)
677 dotproduct2 += elementBaxis[i]*faceBaxis[i];
681 ASSERTL1(fabs(elementBaxis_length*faceBaxis_length - fabs(dotproduct2)) <
683 "These vectors should be parallel");
687 if( dotproduct2 < 0.0 )
702 dotproduct1 += elementAaxis[i]*faceBaxis[i];
706 ASSERTL1(fabs(elementAaxis_length*faceBaxis_length - fabs(dotproduct1)) <
708 "These vectors should be parallel");
712 if(dotproduct1 < 0.0)
721 dotproduct2 += elementBaxis[i]*faceAaxis[i];
725 ASSERTL1(fabs(elementBaxis_length*faceAaxis_length - fabs(dotproduct2)) <
727 "These vectors should be parallel");
729 if( dotproduct2 < 0.0 )
735 orientation = orientation + 5;
748 for (
int i = 0; i < 4; ++i)
750 m_faces[i]->Reset(curvedEdges, curvedFaces);
765 int order0 = *max_element(tmp.begin(), tmp.end());
768 tmp.push_back(order0);
771 int order1 = *max_element(tmp.begin(), tmp.end());
774 tmp.push_back(order0);
775 tmp.push_back(order1);
779 int order2 = *max_element(tmp.begin(), tmp.end());
virtual int v_GetEdgeFaceMap(const int i, const int j) const
StdRegions::StdExpansionSharedPtr m_xmap
#define ASSERTL0(condition, msg)
void SetUpFaceOrientation()
std::vector< StdRegions::Orientation > m_eorient
Principle Modified Functions .
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
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...
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
NekDouble dot(PointGeom &a)
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_GetDir(const int faceidx, const int facedir) const
void Sub(PointGeom &a, PointGeom &b)
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
static const int kNtfaces
StdRegions::StdExpansionSharedPtr GetXmap() const
std::vector< StdRegions::Orientation > m_forient
Gauss Radau pinned at x=-1, .
void SetUpEdgeOrientation()
static const unsigned int EdgeFaceConnectivity[6][2]
virtual int v_GetNumFaces() const
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
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Principle Modified Functions .
virtual int v_GetNumEdges() const
static const unsigned int VertexFaceConnectivity[4][3]
void SetUpXmap()
Set up the m_xmap object by determining the order of each direction from derived faces.
Defines a specification for a set of points.
virtual int v_GetVertexEdgeMap(const int i, const int j) const
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
static const int kNqfaces
PointGeomSharedPtr GetVertex(int i) const
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
virtual NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Get Local cartesian points.
void Mult(PointGeom &a, PointGeom &b)
const Geometry1DSharedPtr GetEdge(int i) const
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Geometry is straight-sided with constant geometric factors.
virtual int v_GetVertexFaceMap(const int i, const int j) const
Gauss Radau pinned at x=-1, .
virtual int v_GetNumVerts() const
boost::shared_ptr< TriGeom > TriGeomSharedPtr
LibUtilities::ShapeType m_shapeType
NekDouble dist(PointGeom &a)
boost::unordered_map< int, CurveSharedPtr > CurveMap
void SetUpCoeffs(const int nCoeffs)
Initialise the m_coeffs array.
static const unsigned int VertexEdgeConnectivity[4][3]
GeomFactorsSharedPtr GetMetricInfo()
void SetUpLocalVertices()
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
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.