49 namespace SpatialDomains
51 Geometry3D::Geometry3D()
55 Geometry3D::Geometry3D(
const int coordim):
59 "Coordinate dimension should be at least 3 for a 3D geometry.");
108 const int MaxIterations = 51;
123 NekDouble derx_1, derx_2, derx_3, dery_1, dery_2, dery_3,
124 derz_1, derz_2, derz_3, jac;
127 NekDouble init0 = Lcoords[0], init1 = Lcoords[1], init2 = Lcoords[2];
140 m_xmap->PhysDeriv(ptsx,DxD1,DxD2,DxD3);
141 m_xmap->PhysDeriv(ptsy,DyD1,DyD2,DyD3);
142 m_xmap->PhysDeriv(ptsz,DzD1,DzD2,DzD3);
150 while(cnt++ < MaxIterations)
153 m_xmap->LocCoordToLocCollapsed(Lcoords,eta);
154 I[0] =
m_xmap->GetBasis(0)->GetI(eta );
155 I[1] =
m_xmap->GetBasis(1)->GetI(eta+1);
156 I[2] =
m_xmap->GetBasis(2)->GetI(eta+2);
159 xmap =
m_xmap->PhysEvaluate(I, ptsx);
160 ymap =
m_xmap->PhysEvaluate(I, ptsy);
161 zmap =
m_xmap->PhysEvaluate(I, ptsz);
163 F1 = coords[0] - xmap;
164 F2 = coords[1] - ymap;
165 F3 = coords[2] - zmap;
167 if(F1*F1 + F2*F2 + F3*F3 < ScaledTol)
169 resid = sqrt(F1*F1 + F2*F2 + F3*F3);
174 derx_1 =
m_xmap->PhysEvaluate(I, DxD1);
175 derx_2 =
m_xmap->PhysEvaluate(I, DxD2);
176 derx_3 =
m_xmap->PhysEvaluate(I, DxD3);
177 dery_1 =
m_xmap->PhysEvaluate(I, DyD1);
178 dery_2 =
m_xmap->PhysEvaluate(I, DyD2);
179 dery_3 =
m_xmap->PhysEvaluate(I, DyD3);
180 derz_1 =
m_xmap->PhysEvaluate(I, DzD1);
181 derz_2 =
m_xmap->PhysEvaluate(I, DzD2);
182 derz_3 =
m_xmap->PhysEvaluate(I, DzD3);
184 jac = derx_1*(dery_2*derz_3 - dery_3*derz_2)
185 - derx_2*(dery_1*derz_3 - dery_3*derz_1)
186 + derx_3*(dery_1*derz_2 - dery_2*derz_1);
190 Lcoords[0] = Lcoords[0]
191 +((dery_2*derz_3 - dery_3*derz_2)*(coords[0]-xmap)
192 - (derx_2*derz_3 - derx_3*derz_2)*(coords[1]-ymap)
193 + (derx_2*dery_3 - derx_3*dery_2)*(coords[2]-zmap)
196 Lcoords[1] = Lcoords[1]
197 -((dery_1*derz_3 - dery_3*derz_1)*(coords[0]-xmap)
198 - (derx_1*derz_3 - derx_3*derz_1)*(coords[1]-ymap)
199 + (derx_1*dery_3 - derx_3*dery_1)*(coords[2]-zmap)
202 Lcoords[2] = Lcoords[2]
203 +((dery_1*derz_2 - dery_2*derz_1)*(coords[0]-xmap)
204 - (derx_1*derz_2 - derx_2*derz_1)*(coords[1]-ymap)
205 + (derx_1*dery_2 - derx_2*dery_1)*(coords[2]-zmap)
208 if (fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv ||
209 fabs(Lcoords[0]) > LcoordDiv)
215 resid = sqrt(F1*F1 + F2*F2 + F3*F3);
217 if(cnt >= MaxIterations)
220 m_xmap->LocCoordToLocCollapsed(Lcoords,collCoords);
223 if((collCoords[0] >= -1.0 && collCoords[0] <= 1.0)&&
224 (collCoords[1] >= -1.0 && collCoords[1] <= 1.0)&&
225 (collCoords[2] >= -1.0 && collCoords[2] <= 1.0))
227 std::ostringstream ss;
229 ss <<
"Reached MaxIterations (" << MaxIterations <<
") in Newton iteration ";
230 ss <<
"Init value ("<< setprecision(4) << init0 <<
"," << init1<<
"," << init2 <<
") ";
231 ss <<
"Fin value ("<<Lcoords[0] <<
"," << Lcoords[1]<<
"," << Lcoords[2] <<
") ";
232 ss <<
"Resid = " << resid <<
" Tolerance = " << sqrt(ScaledTol) ;
258 int nFaceCoeffs =
m_faces[i]->GetXmap()->GetNcoeffs();
265 m_xmap->GetFaceToElementMap(
272 m_xmap->GetFaceToElementMap(
283 for (k = 0; k < nFaceCoeffs; k++)
308 if (
m_xmap->GetBasisNumModes(0) != 2 ||
309 m_xmap->GetBasisNumModes(1) != 2 ||
310 m_xmap->GetBasisNumModes(2) != 2)
330 "Geometry is not in physical space");
335 return m_xmap->PhysEvaluate(Lcoord, tmp);
357 "Vertex ID must be between 0 and "+
358 boost::lexical_cast<
string>(
m_verts.size() - 1));
376 "Edge ID must be between 0 and "+
377 boost::lexical_cast<
string>(
m_edges.size() - 1));
388 "Edge ID must be between 0 and "+
389 boost::lexical_cast<
string>(
m_edges.size() - 1));
399 "Edge ID must be between 0 and "+
400 boost::lexical_cast<
string>(
m_edges.size() - 1));
409 ASSERTL2((i >=0) && (i <= 5),
"Edge id must be between 0 and 4");
419 "Face ID must be between 0 and "+
420 boost::lexical_cast<
string>(
m_faces.size() - 1));
430 "Face ID must be between 0 and "+
431 boost::lexical_cast<
string>(
m_faces.size() - 1));
449 return m_xmap->GetBasis(i);
465 for (i=0,edgeIter =
m_edges.begin(); edgeIter !=
m_edges.end(); ++edgeIter,++i)
467 if (*edgeIter == edge)
515 std::list<CompToElmt>::const_iterator def;
StdRegions::StdExpansionSharedPtr m_xmap
#define ASSERTL0(condition, msg)
std::vector< StdRegions::Orientation > m_eorient
Base class for shape geometry information.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
int GetDir(const int faceidx, const int facedir) const
Returns the element coordinate direction corresponding to a given face coordinate direction...
GeomFactorsSharedPtr m_geomFactors
virtual PointGeomSharedPtr v_GetVertex(int i) const
Return vertex i in this element.
std::list< CompToElmt > m_elmtmap
virtual int v_GetShapeDim() const
Return the dimension of this element.
Structure holding graphvertexobject id and local element facet id.
virtual void v_SetOwnData()
int GetEid(int i) const
Return the ID of edge i in this element.
Geometry2DSharedPtr GetFace(int i)
Return face i in this element.
StdRegions::StdExpansionSharedPtr GetXmap() const
GeomState m_geomFactorsState
std::vector< StdRegions::Orientation > m_forient
virtual int v_WhichEdge(SegGeomSharedPtr edge)
Return the local ID of a given edge.
virtual int v_GetFid(int i) const
Return the ID of face i in this element.
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)
virtual int v_WhichFace(Geometry2DSharedPtr face)
Return the local ID of a given face.
virtual void v_AddElmtConnected(int gvo_id, int locid)
boost::shared_ptr< SegGeom > SegGeomSharedPtr
virtual bool v_IsElmtConnected(int gvo_id, int locid) const
virtual int v_GetDir(const int faceidx, const int facedir) const =0
virtual void v_GenGeomFactors()
virtual int v_GetVid(int i) const
Return the vertex ID of vertex i.
virtual StdRegions::Orientation v_GetEorient(const int i) const
Return the orientation of edge i in this element.
virtual const LibUtilities::BasisSharedPtr v_GetBasis(const int i)
Return the j-th basis of the i-th co-ordinate dimension.
virtual int v_GetEid() const
Return the ID of this element.
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
virtual StdRegions::Orientation v_GetForient(const int i) const
Return the orientation of face i in this element.
virtual const Geometry2DSharedPtr v_GetFace(int i) const
Return face i in this element.
Array< OneD, Array< OneD, NekDouble > > m_coeffs
#define WARNINGL1(condition, msg)
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Geometry is straight-sided with constant geometric factors.
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
virtual NekDouble v_GetCoord(const int i, const Array< OneD, const NekDouble > &Lcoord)
Given local collapsed coordinate Lcoord return the value of physical coordinate in direction i...
virtual const SegGeomSharedPtr v_GetEdge(int i) const
Return edge i of this element.
Geometric information has been generated.
GeomType
Indicates the type of element geometry.
boost::shared_ptr< Basis > BasisSharedPtr
GeomState m_state
enum identifier to determine if quad points are filled
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
#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
boost::shared_ptr< PointGeom > PointGeomSharedPtr
virtual int v_NumElmtConnected() const