48 namespace SpatialDomains
50 const unsigned int PrismGeom::VertexEdgeConnectivity[6][3] = {
51 {0,3,4},{0,1,5},{1,2,6},{2,3,7},{4,5,8},{6,7,8}};
52 const unsigned int PrismGeom::VertexFaceConnectivity[6][3] = {
53 {0,1,4},{0,1,2},{0,2,3},{0,3,4},{1,2,4},{2,3,4}};
54 const unsigned int PrismGeom::EdgeFaceConnectivity [9][2] = {
55 {0,1},{0,2},{0,3},{0,4},{1,4},{1,2},{2,3},{3,4},{2,4}};
57 PrismGeom::PrismGeom()
63 Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
110 else if (faceidx == 1 || faceidx == 3)
151 ASSERTL1(gloCoord.num_elements() == 3,
152 "Three dimensional geometry expects three coordinates.");
168 for(i = 0; i < 3; ++i)
172 mincoord[i] =
Vmath::Vmin(pts.num_elements(),pts,1);
173 maxcoord[i] =
Vmath::Vmax(pts.num_elements(),pts,1);
175 diff = max(maxcoord[i] - mincoord[i],diff);
178 for(i = 0; i < 3; ++i)
180 if((gloCoord[i] < mincoord[i] - 0.2*diff)||
181 (gloCoord[i] > maxcoord[i] + 0.2*diff))
192 if (locCoord[0] >= -(1+tol) && locCoord[1] >= -(1+tol) &&
193 locCoord[2] >= -(1+tol) && locCoord[1] <= (1+tol) &&
194 locCoord[0] + locCoord[2] <= tol)
203 for(
int i = 0; i < 3; ++i)
205 if(locCoord[i] <-(1+tol))
207 locCoord[i] = -(1+tol);
210 if(locCoord[i] > (1+tol))
231 if (
m_xmap->GetBasisNumModes(0) != 2 ||
232 m_xmap->GetBasisNumModes(1) != 2 ||
233 m_xmap->GetBasisNumModes(2) != 2 )
243 const unsigned int faceVerts[3][4] =
248 for(f = 0; f < 3; f++)
253 if( fabs( (*
m_verts[ faceVerts[f][0] ])(i) -
254 (*
m_verts[ faceVerts[f][1] ])(i) +
255 (*
m_verts[ faceVerts[f][2] ])(i) -
256 (*
m_verts[ faceVerts[f][3] ])(i) )
299 cp1030.
Mult(e10,e30);
300 cp3040.
Mult(e30,e40);
301 cp4010.
Mult(e40,e10);
310 Lcoords[0] = 2.0*beta - 1.0;
311 Lcoords[1] = 2.0*gamma - 1.0;
312 Lcoords[2] = 2.0*delta - 1.0;
315 for(
int i = 0; i < 6; ++i)
348 ptdist = sqrt(tmp1[min_i]);
351 int qa = za.num_elements(), qb = zb.num_elements();
352 Lcoords[2] = zc[min_i/(qa*qb)];
353 min_i = min_i%(qa*qb);
354 Lcoords[1] = zb[min_i/qa];
355 Lcoords[0] = za[min_i%qa];
358 Lcoords[0] = (1.0+Lcoords[0])*(1.0-Lcoords[2])/2 - 1.0;
392 for (f = 1; f < 5 ; f++)
394 int nEdges =
m_faces[f]->GetNumEdges();
396 for (i = 0; i < 4; i++)
398 for (j = 0; j < nEdges; j++)
411 std::ostringstream errstrm;
412 errstrm <<
"Connected faces do not share an edge. Faces ";
418 std::ostringstream errstrm;
419 errstrm <<
"Connected faces share more than one edge. Faces ";
427 for(i = 0; i < 3; i++)
429 for(j = 0; j < 4; j++)
441 std::ostringstream errstrm;
442 errstrm <<
"Connected faces do not share an edge. Faces ";
448 std::ostringstream errstrm;
449 errstrm <<
"Connected faces share more than one edge. Faces ";
454 for(f = 1; f < 4 ; f++)
457 for(i = 0; i <
m_faces[f]->GetNumEdges(); i++)
459 for(j = 0; j <
m_faces[f+1]->GetNumEdges(); j++)
472 std::ostringstream errstrm;
473 errstrm <<
"Connected faces do not share an edge. Faces ";
479 std::ostringstream errstrm;
480 errstrm <<
"Connected faces share more than one edge. Faces ";
488 for(i = 0; i < 4; i++)
490 for(j = 0; j < 4; j++)
503 std::ostringstream errstrm;
504 errstrm <<
"Connected faces do not share an edge. Faces ";
510 std::ostringstream errstrm;
511 errstrm <<
"Connected faces share more than one edge. Faces ";
535 std::ostringstream errstrm;
536 errstrm <<
"Connected edges do not share a vertex. Edges ";
542 for(
int i = 1; i < 3; i++)
554 std::ostringstream errstrm;
555 errstrm <<
"Connected edges do not share a vertex. Edges ";
556 errstrm <<
m_edges[i]->GetEid() <<
", " <<
m_edges[i-1]->GetEid();
577 std::ostringstream errstrm;
578 errstrm <<
"Connected edges do not share a vertex. Edges ";
579 errstrm <<
m_edges[8]->GetEid();
589 const unsigned int edgeVerts[
kNedges][2] =
614 ASSERTL0(
false,
"Could not find matching vertex for the edge");
643 unsigned int baseVertex;
665 unsigned int orientation;
671 elementAaxis_length = 0.0;
672 elementBaxis_length = 0.0;
673 faceAaxis_length = 0.0;
674 faceBaxis_length = 0.0;
679 baseVertex =
m_faces[f]->GetVid(0);
692 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
693 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
701 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
702 elementBaxis[i] = (*
m_verts[ faceVerts[f][3] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
705 else if( baseVertex ==
m_verts[ faceVerts[f][1] ]->
GetVid() )
709 elementAaxis[i] = (*
m_verts[ faceVerts[f][1] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
710 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][1] ])[i];
713 else if( baseVertex ==
m_verts[ faceVerts[f][2] ]->
GetVid() )
717 elementAaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][3] ])[i];
718 elementBaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][1] ])[i];
721 else if( baseVertex ==
m_verts[ faceVerts[f][3] ]->
GetVid() )
725 elementAaxis[i] = (*
m_verts[ faceVerts[f][2] ])[i] - (*
m_verts[ faceVerts[f][3] ])[i];
726 elementBaxis[i] = (*
m_verts[ faceVerts[f][3] ])[i] - (*
m_verts[ faceVerts[f][0] ])[i];
731 ASSERTL0(
false,
"Could not find matching vertex for the face");
738 int v =
m_faces[f]->GetNumVerts()-1;
742 elementAaxis_length += pow(elementAaxis[i],2);
743 elementBaxis_length += pow(elementBaxis[i],2);
744 faceAaxis_length += pow(faceAaxis[i],2);
745 faceBaxis_length += pow(faceBaxis[i],2);
748 elementAaxis_length = sqrt(elementAaxis_length);
749 elementBaxis_length = sqrt(elementBaxis_length);
750 faceAaxis_length = sqrt(faceAaxis_length);
751 faceBaxis_length = sqrt(faceBaxis_length);
757 dotproduct1 += elementAaxis[i]*faceAaxis[i];
767 if(dotproduct1 < 0.0)
775 dotproduct2 += elementBaxis[i]*faceBaxis[i];
785 if( dotproduct2 < 0.0 )
800 dotproduct1 += elementAaxis[i]*faceBaxis[i];
804 if (fabs(elementAaxis_length*faceBaxis_length
807 cout <<
"Warning: Prism axes not parallel" << endl;
812 if(dotproduct1 < 0.0)
821 dotproduct2 += elementBaxis[i]*faceAaxis[i];
825 if (fabs(elementBaxis_length*faceAaxis_length
828 cout <<
"Warning: Prism axes not parallel" << endl;
831 if( dotproduct2 < 0.0 )
837 orientation = orientation + 5;
849 for (
int i = 0; i < 5; ++i)
851 m_faces[i]->Reset(curvedEdges, curvedFaces);
872 order0 = *max_element(tmp.begin(), tmp.end());
878 order0 = *max_element(tmp.begin(), tmp.end());
887 order1 = *max_element(tmp.begin(), tmp.end());
895 order1 = *max_element(tmp.begin(), tmp.end());
899 tmp.push_back(order0);
900 tmp.push_back(order1);
905 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.