38 #include <boost/shared_ptr.hpp>
46 namespace SpatialDomains
48 Geometry2D::Geometry2D()
52 Geometry2D::Geometry2D(
const int coordim):
56 "Coordinate dimension should be at least 2 for a 2D geometry");
111 const int MaxIterations = 51;
125 NekDouble derx_1, derx_2, dery_1, dery_2,jac;
128 NekDouble init0 = Lcoords[0], init1 = Lcoords[1];
136 m_xmap->PhysDeriv(ptsx,DxD1,DxD2);
137 m_xmap->PhysDeriv(ptsy,DyD1,DyD2);
145 while(cnt++ < MaxIterations)
148 m_xmap->LocCoordToLocCollapsed(Lcoords,eta);
149 I[0] =
m_xmap->GetBasis(0)->GetI(eta);
150 I[1] =
m_xmap->GetBasis(1)->GetI(eta+1);
153 xmap =
m_xmap->PhysEvaluate(I, ptsx);
154 ymap =
m_xmap->PhysEvaluate(I, ptsy);
156 F1 = coords[0] - xmap;
157 F2 = coords[1] - ymap;
159 if(F1*F1 + F2*F2 < ScaledTol)
161 resid = sqrt(F1*F1 + F2*F2);
166 derx_1 =
m_xmap->PhysEvaluate(I, DxD1);
167 derx_2 =
m_xmap->PhysEvaluate(I, DxD2);
168 dery_1 =
m_xmap->PhysEvaluate(I, DyD1);
169 dery_2 =
m_xmap->PhysEvaluate(I, DyD2);
171 jac = dery_2*derx_1 - dery_1*derx_2;
175 Lcoords[0] = Lcoords[0] + (dery_2*(coords[0]-xmap) -
176 derx_2*(coords[1]-ymap))/jac;
178 Lcoords[1] = Lcoords[1] + ( - dery_1*(coords[0]-xmap)
179 + derx_1*(coords[1]-ymap))/jac;
181 if(fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv)
187 resid = sqrt(F1*F1 + F2*F2);
189 if(cnt >= MaxIterations)
192 m_xmap->LocCoordToLocCollapsed(Lcoords,collCoords);
195 if((collCoords[0] >= -1.0 && collCoords[0] <= 1.0)&&
196 (collCoords[1] >= -1.0 && collCoords[1] <= 1.0))
198 std::ostringstream ss;
200 ss <<
"Reached MaxIterations (" << MaxIterations
201 <<
") in Newton iteration ";
202 ss <<
"Init value ("<< setprecision(4) << init0 <<
","
203 << init1<<
"," <<
") ";
204 ss <<
"Fin value ("<<Lcoords[0] <<
"," << Lcoords[1]
206 ss <<
"Resid = " << resid <<
" Tolerance = "
217 "This function is only valid for shape type geometries");
224 "This function is only valid for shape type geometries");
233 "This function is only valid for shape type geometries");
240 "This function is only valid for shape type geometries");
248 "This function is only valid for shape type geometries");
256 "This function is only valid for shape type geometries");
263 "This function is only valid for shape type geometries");
271 "This function is only valid for shape type geometries");
278 "This function is only valid for shape type geometries");
285 "This function is only valid for shape type geometries");
292 "This function is only valid for shape type geometries");
306 "This function has not been defined for this geometry");
StdRegions::StdExpansionSharedPtr m_xmap
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
int WhichEdge(SegGeomSharedPtr edge)
Base class for shape geometry information.
GeomFactorsSharedPtr m_geomFactors
virtual int v_WhichFace(Geometry2DSharedPtr face)
const LibUtilities::BasisSharedPtr GetEdgeBasis(const int i)
StdRegions::Orientation GetCartesianEorient(const int i) const
const Geometry2DSharedPtr GetFace(int i) const
virtual int v_WhichEdge(SegGeomSharedPtr edge)
virtual PointGeomSharedPtr v_GetVertex(int i) const
virtual int v_GetShapeDim() const
virtual StdRegions::Orientation v_GetCartesianEorient(const int i) const
virtual const LibUtilities::BasisSharedPtr v_GetEdgeBasis(const int i)
boost::shared_ptr< SegGeom > SegGeomSharedPtr
virtual StdRegions::Orientation v_GetEorient(const int i) const
void NewtonIterationForLocCoord(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &ptsx, const Array< OneD, const NekDouble > &ptsy, Array< OneD, NekDouble > &Lcoords, NekDouble &resid)
virtual bool v_ContainsPoint(const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0)
virtual int v_GetFid() const
StdRegions::Orientation GetFaceOrient(const int i) const
int WhichFace(Geometry2DSharedPtr face)
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
virtual StdRegions::Orientation v_GetFaceOrient(const int i) const
virtual const Geometry2DSharedPtr v_GetFace(int i) const
#define WARNINGL1(condition, msg)
const Geometry1DSharedPtr GetEdge(int i) const
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
virtual const Geometry1DSharedPtr v_GetEdge(int i) const
boost::shared_ptr< Basis > BasisSharedPtr
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
int m_coordim
coordinate dimension
virtual int v_GetEid(int i) const
boost::shared_ptr< PointGeom > PointGeomSharedPtr