46 namespace SpatialDomains
49 {0,3,4},{0,1,5},{1,2,6},{2,3,7},{4,5,8},{6,7,8}};
51 {0,1,4},{0,1,2},{0,2,3},{0,3,4},{1,2,4},{2,3,4}};
53 {0,1},{0,2},{0,3},{0,4},{1,4},{1,2},{2,3},{3,4},{2,4}};
61 Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
108 else if (faceidx == 1 || faceidx == 3)
149 ASSERTL1(gloCoord.num_elements() == 3,
150 "Three dimensional geometry expects three coordinates.");
166 for(i = 0; i < 3; ++i)
170 mincoord[i] =
Vmath::Vmin(pts.num_elements(),pts,1);
171 maxcoord[i] =
Vmath::Vmax(pts.num_elements(),pts,1);
173 diff = max(maxcoord[i] - mincoord[i],diff);
176 for(i = 0; i < 3; ++i)
178 if((gloCoord[i] < mincoord[i] - 0.2*diff)||
179 (gloCoord[i] > maxcoord[i] + 0.2*diff))
190 if (locCoord[0] >= -(1+tol) && locCoord[1] >= -(1+tol) &&
191 locCoord[2] >= -(1+tol) && locCoord[1] <= (1+tol) &&
192 locCoord[0] + locCoord[2] <= tol)
201 for(
int i = 0; i < 3; ++i)
203 if(locCoord[i] <-(1+tol))
205 locCoord[i] = -(1+tol);
208 if(locCoord[i] > (1+tol))
229 if (
m_xmap->GetBasisNumModes(0) != 2 ||
230 m_xmap->GetBasisNumModes(1) != 2 ||
231 m_xmap->GetBasisNumModes(2) != 2 )
241 const unsigned int faceVerts[3][4] =
246 for(f = 0; f < 3; f++)
251 if( fabs( (*
m_verts[ faceVerts[f][0] ])(i) -
252 (*
m_verts[ faceVerts[f][1] ])(i) +
253 (*
m_verts[ faceVerts[f][2] ])(i) -
254 (*
m_verts[ faceVerts[f][3] ])(i) )
297 cp1030.
Mult(e10,e30);
298 cp3040.
Mult(e30,e40);
299 cp4010.
Mult(e40,e10);
308 Lcoords[0] = 2.0*beta - 1.0;
309 Lcoords[1] = 2.0*gamma - 1.0;
310 Lcoords[2] = 2.0*delta - 1.0;
313 for(
int i = 0; i < 6; ++i)
346 ptdist = sqrt(tmp1[min_i]);
349 int qa = za.num_elements(), qb = zb.num_elements();
350 Lcoords[2] = zc[min_i/(qa*qb)];
351 min_i = min_i%(qa*qb);
352 Lcoords[1] = zb[min_i/qa];
353 Lcoords[0] = za[min_i%qa];
356 Lcoords[0] = (1.0+Lcoords[0])*(1.0-Lcoords[2])/2 - 1.0;
390 for (f = 1; f < 5 ; f++)
392 int nEdges =
m_faces[f]->GetNumEdges();
394 for (i = 0; i < 4; i++)
396 for (j = 0; j < nEdges; j++)
409 std::ostringstream errstrm;
410 errstrm <<
"Connected faces do not share an edge. Faces ";
416 std::ostringstream errstrm;
417 errstrm <<
"Connected faces share more than one edge. Faces ";
425 for(i = 0; i < 3; i++)
427 for(j = 0; j < 4; j++)
439 std::ostringstream errstrm;
440 errstrm <<
"Connected faces do not share an edge. Faces ";
446 std::ostringstream errstrm;
447 errstrm <<
"Connected faces share more than one edge. Faces ";
452 for(f = 1; f < 4 ; f++)
455 for(i = 0; i <
m_faces[f]->GetNumEdges(); i++)
457 for(j = 0; j <
m_faces[f+1]->GetNumEdges(); j++)
470 std::ostringstream errstrm;
471 errstrm <<
"Connected faces do not share an edge. Faces ";
477 std::ostringstream errstrm;
478 errstrm <<
"Connected faces share more than one edge. Faces ";
486 for(i = 0; i < 4; i++)
488 for(j = 0; j < 4; j++)
501 std::ostringstream errstrm;
502 errstrm <<
"Connected faces do not share an edge. Faces ";
508 std::ostringstream errstrm;
509 errstrm <<
"Connected faces share more than one edge. Faces ";
533 std::ostringstream errstrm;
534 errstrm <<
"Connected edges do not share a vertex. Edges ";
540 for(
int i = 1; i < 3; i++)
552 std::ostringstream errstrm;
553 errstrm <<
"Connected edges do not share a vertex. Edges ";
554 errstrm <<
m_edges[i]->GetEid() <<
", " <<
m_edges[i-1]->GetEid();
575 std::ostringstream errstrm;
576 errstrm <<
"Connected edges do not share a vertex. Edges ";
577 errstrm <<
m_edges[8]->GetEid();
587 const unsigned int edgeVerts[
kNedges][2] =
612 ASSERTL0(
false,
"Could not find matching vertex for the edge");
641 unsigned int baseVertex;
663 unsigned int orientation;
669 elementAaxis_length = 0.0;
670 elementBaxis_length = 0.0;
671 faceAaxis_length = 0.0;
672 faceBaxis_length = 0.0;
677 baseVertex =
m_faces[f]->GetVid(0);
690 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
691 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
699 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
700 elementBaxis[i] = (*
m_verts[ faceVerts[f][3] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
703 else if( baseVertex ==
m_verts[ faceVerts[f][1] ]->
GetVid() )
707 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
708 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][1] ])[i];
711 else if( baseVertex ==
m_verts[ faceVerts[f][2] ]->
GetVid() )
715 elementAaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][3] ])[i];
716 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][1] ])[i];
719 else if( baseVertex ==
m_verts[ faceVerts[f][3] ]->
GetVid() )
723 elementAaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][3] ])[i];
724 elementBaxis[i] = (*
m_verts[ faceVerts[f][3] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
729 ASSERTL0(
false,
"Could not find matching vertex for the face");
736 int v =
m_faces[f]->GetNumVerts()-1;
740 elementAaxis_length += pow(elementAaxis[i],2);
741 elementBaxis_length += pow(elementBaxis[i],2);
742 faceAaxis_length += pow(faceAaxis[i],2);
743 faceBaxis_length += pow(faceBaxis[i],2);
746 elementAaxis_length = sqrt(elementAaxis_length);
747 elementBaxis_length = sqrt(elementBaxis_length);
748 faceAaxis_length = sqrt(faceAaxis_length);
749 faceBaxis_length = sqrt(faceBaxis_length);
755 dotproduct1 += elementAaxis[i]*faceAaxis[i];
765 if(dotproduct1 < 0.0)
773 dotproduct2 += elementBaxis[i]*faceBaxis[i];
783 if( dotproduct2 < 0.0 )
798 dotproduct1 += elementAaxis[i]*faceBaxis[i];
802 if (fabs(elementAaxis_length*faceBaxis_length
805 cout <<
"Warning: Prism axes not parallel" << endl;
810 if(dotproduct1 < 0.0)
819 dotproduct2 += elementBaxis[i]*faceAaxis[i];
823 if (fabs(elementBaxis_length*faceAaxis_length
826 cout <<
"Warning: Prism axes not parallel" << endl;
829 if( dotproduct2 < 0.0 )
835 orientation = orientation + 5;
847 for (
int i = 0; i < 5; ++i)
849 m_faces[i]->Reset(curvedEdges, curvedFaces);
870 order0 = *max_element(tmp.begin(), tmp.end());
876 order0 = *max_element(tmp.begin(), tmp.end());
885 order1 = *max_element(tmp.begin(), tmp.end());
893 order1 = *max_element(tmp.begin(), tmp.end());
897 tmp.push_back(order0);
898 tmp.push_back(order1);
903 int order2 = *max_element(tmp.begin(), tmp.end());
StdRegions::StdExpansionSharedPtr m_xmap
#define ASSERTL0(condition, msg)
std::vector< StdRegions::Orientation > m_eorient
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
GeomFactorsSharedPtr m_geomFactors
void SetUpFaceOrientation()
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
virtual int v_GetVertexFaceMap(const int i, const int j) const
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)
void SetUpLocalVertices()
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 .
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.
void SetUpXmap()
Set up the m_xmap object by determining the order of each direction from derived faces.
virtual NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
StdRegions::StdExpansionSharedPtr GetXmap() const
GeomState m_geomFactorsState
virtual int v_GetDir(const int faceidx, const int facedir) const
static const unsigned int VertexFaceConnectivity[6][3]
std::vector< StdRegions::Orientation > m_forient
Gauss Radau pinned at x=-1, .
virtual int v_GetNumEdges() 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
virtual int v_GetEdgeFaceMap(const int i, const int j) const
virtual int v_GetNumVerts() const
static const int kNqfaces
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Principle Modified Functions .
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
void SetUpEdgeOrientation()
Defines a specification for a set of points.
static const unsigned int EdgeFaceConnectivity[9][2]
virtual void v_GenGeomFactors()
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
PointGeomSharedPtr GetVertex(int i) const
virtual int v_GetVertexEdgeMap(const int i, const int j) const
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
void Mult(PointGeom &a, PointGeom &b)
const Geometry1DSharedPtr GetEdge(int i) const
Array< OneD, Array< OneD, NekDouble > > m_coeffs
virtual int v_GetNumFaces() const
Geometry is straight-sided with constant geometric factors.
static const unsigned int VertexEdgeConnectivity[6][3]
LibUtilities::ShapeType m_shapeType
NekDouble dist(PointGeom &a)
boost::unordered_map< int, CurveSharedPtr > CurveMap
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...
Geometric information has been generated.
void SetUpCoeffs(const int nCoeffs)
Initialise the m_coeffs array.
GeomFactorsSharedPtr GetMetricInfo()
GeomType
Indicates the type of element geometry.
static const int kNtfaces
#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.