Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Geometry2D.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Geometry2D.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: 2D geometry information
33 //
34 //
35 ////////////////////////////////////////////////////////////////////////////////
37 #include <SpatialDomains/SegGeom.h>
38 #include <boost/shared_ptr.hpp>
39 
40 #include <iomanip>
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46  namespace SpatialDomains
47  {
48  Geometry2D::Geometry2D()
49  {
50  }
51 
52  Geometry2D::Geometry2D(const int coordim):
53  Geometry(coordim)
54  {
55  ASSERTL0(m_coordim > 1,
56  "Coordinate dimension should be at least 2 for a 2D geometry");
57  }
58 
60  {
61  }
62 
63  int Geometry2D::GetFid() const
64  {
65  return v_GetFid();
66  }
67 
69  {
70  return v_GetEdgeBasis(i);
71  }
72 
74  {
75  return v_GetFace(i);
76  }
77 
79  {
80  return v_GetFaceOrient(i);
81  }
82 
84  {
85  return v_GetEdge(i);
86  }
87 
89  {
90  return v_GetCartesianEorient(i);
91  }
92 
94  {
95  return v_WhichEdge(edge);
96  }
97 
99  {
100  return v_WhichFace(face);
101  }
102 
104  const Array<OneD, const NekDouble> &coords,
105  const Array<OneD, const NekDouble> &ptsx,
106  const Array<OneD, const NekDouble> &ptsy,
107  Array<OneD, NekDouble> &Lcoords,
108  NekDouble &resid)
109  {
110  // Maximum iterations for convergence
111  const int MaxIterations = 51;
112  // |x-xp|^2 < EPSILON error tolerance
113  const NekDouble Tol = 1.e-8;
114  // |r,s| > LcoordDIV stop the search
115  const NekDouble LcoordDiv = 15.0;
116 
118  m_geomFactors->GetJac(m_xmap->GetPointsKeys());
119 
120  NekDouble ScaledTol = Vmath::Vsum(Jac.num_elements(),Jac,1)/
121  ((NekDouble)Jac.num_elements());
122  ScaledTol *= Tol;
123 
124  NekDouble xmap,ymap, F1,F2;
125  NekDouble derx_1, derx_2, dery_1, dery_2,jac;
126 
127  // save intiial guess for later reference if required.
128  NekDouble init0 = Lcoords[0], init1 = Lcoords[1];
129 
130  Array<OneD, NekDouble> DxD1(ptsx.num_elements());
131  Array<OneD, NekDouble> DxD2(ptsx.num_elements());
132  Array<OneD, NekDouble> DyD1(ptsx.num_elements());
133  Array<OneD, NekDouble> DyD2(ptsx.num_elements());
134 
135  // Ideally this will be stored in m_geomfactors
136  m_xmap->PhysDeriv(ptsx,DxD1,DxD2);
137  m_xmap->PhysDeriv(ptsy,DyD1,DyD2);
138 
139  int cnt=0;
141  Array<OneD, NekDouble> eta(2);
142 
143  F1 = F2 = 2000; // Starting value of Function
144 
145  while(cnt++ < MaxIterations)
146  {
147  // evaluate lagrange interpolant at Lcoords
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);
151 
152  //calculate the global point `corresponding to Lcoords
153  xmap = m_xmap->PhysEvaluate(I, ptsx);
154  ymap = m_xmap->PhysEvaluate(I, ptsy);
155 
156  F1 = coords[0] - xmap;
157  F2 = coords[1] - ymap;
158 
159  if(F1*F1 + F2*F2 < ScaledTol)
160  {
161  resid = sqrt(F1*F1 + F2*F2);
162  break;
163  }
164 
165  //Interpolate derivative metric at Lcoords
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);
170 
171  jac = dery_2*derx_1 - dery_1*derx_2;
172 
173  // use analytical inverse of derivitives which are
174  // also similar to those of metric factors.
175  Lcoords[0] = Lcoords[0] + (dery_2*(coords[0]-xmap) -
176  derx_2*(coords[1]-ymap))/jac;
177 
178  Lcoords[1] = Lcoords[1] + ( - dery_1*(coords[0]-xmap)
179  + derx_1*(coords[1]-ymap))/jac;
180 
181  if(fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv)
182  {
183  break; // lcoords have diverged so stop iteration
184  }
185  }
186 
187  resid = sqrt(F1*F1 + F2*F2);
188 
189  if(cnt >= MaxIterations)
190  {
191  Array<OneD, NekDouble> collCoords(2);
192  m_xmap->LocCoordToLocCollapsed(Lcoords,collCoords);
193 
194  // if coordinate is inside element dump error!
195  if((collCoords[0] >= -1.0 && collCoords[0] <= 1.0)&&
196  (collCoords[1] >= -1.0 && collCoords[1] <= 1.0))
197  {
198  std::ostringstream ss;
199 
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]
205  << "," << ") ";
206  ss << "Resid = " << resid << " Tolerance = "
207  << sqrt(ScaledTol) ;
208 
209  WARNINGL1(cnt < MaxIterations,ss.str());
210  }
211  }
212  }
213 
214  int Geometry2D::v_GetFid() const
215  {
217  "This function is only valid for shape type geometries");
218  return 0;
219  }
220 
222  {
224  "This function is only valid for shape type geometries");
226  return returnval;
227  }
228 
229 
230  int Geometry2D::v_GetEid(int i) const
231  {
233  "This function is only valid for shape type geometries");
234  return 0;
235  }
236 
238  {
240  "This function is only valid for shape type geometries");
241  SegGeomSharedPtr returnval;
242  return returnval;
243  }
244 
246  {
248  "This function is only valid for shape type geometries");
249  PointGeomSharedPtr returnval;
250  return returnval;
251  }
252 
254  {
256  "This function is only valid for shape type geometries");
257  return StdRegions::eForwards;
258  }
259 
261  {
263  "This function is only valid for shape type geometries");
264  Geometry2DSharedPtr returnval;
265  return returnval;
266  }
267 
269  {
271  "This function is only valid for shape type geometries");
273  }
274 
276  {
278  "This function is only valid for shape type geometries");
279  return StdRegions::eForwards;
280  }
281 
283  {
285  "This function is only valid for shape type geometries");
286  return 0;
287  }
288 
290  {
292  "This function is only valid for shape type geometries");
293  return 0;
294  }
295 
297  {
298  return 2;
299  }
300 
302  const Array<OneD, const NekDouble>& gloCoord,
303  NekDouble tol)
304  {
306  "This function has not been defined for this geometry");
307  return false;
308  }
309 
310 
311  }; //end of namespace
312 }; //end of namespace
313 
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:191
int WhichEdge(SegGeomSharedPtr edge)
Definition: Geometry2D.cpp:93
Base class for shape geometry information.
Definition: Geometry.h:76
GeomFactorsSharedPtr m_geomFactors
Definition: Geometry.h:170
virtual int v_WhichFace(Geometry2DSharedPtr face)
Definition: Geometry2D.cpp:289
const LibUtilities::BasisSharedPtr GetEdgeBasis(const int i)
Definition: Geometry2D.cpp:68
StdRegions::Orientation GetCartesianEorient(const int i) const
Definition: Geometry2D.cpp:88
STL namespace.
const Geometry2DSharedPtr GetFace(int i) const
Definition: Geometry2D.cpp:73
virtual int v_WhichEdge(SegGeomSharedPtr edge)
Definition: Geometry2D.cpp:282
virtual PointGeomSharedPtr v_GetVertex(int i) const
Definition: Geometry2D.cpp:245
virtual int v_GetShapeDim() const
Definition: Geometry2D.cpp:296
virtual StdRegions::Orientation v_GetCartesianEorient(const int i) const
Definition: Geometry2D.cpp:275
virtual const LibUtilities::BasisSharedPtr v_GetEdgeBasis(const int i)
Definition: Geometry2D.cpp:221
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
virtual StdRegions::Orientation v_GetEorient(const int i) const
Definition: Geometry2D.cpp:253
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)
Definition: Geometry2D.cpp:103
double NekDouble
virtual bool v_ContainsPoint(const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0)
Definition: Geometry2D.cpp:301
StdRegions::Orientation GetFaceOrient(const int i) const
Definition: Geometry2D.cpp:78
int WhichFace(Geometry2DSharedPtr face)
Definition: Geometry2D.cpp:98
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
virtual StdRegions::Orientation v_GetFaceOrient(const int i) const
Definition: Geometry2D.cpp:268
virtual const Geometry2DSharedPtr v_GetFace(int i) const
Definition: Geometry2D.cpp:260
#define WARNINGL1(condition, msg)
Definition: ErrorUtil.hpp:229
const Geometry1DSharedPtr GetEdge(int i) const
Definition: Geometry2D.cpp:83
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry1D.h:48
virtual const Geometry1DSharedPtr v_GetEdge(int i) const
Definition: Geometry2D.cpp:237
boost::shared_ptr< Basis > BasisSharedPtr
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:737
int m_coordim
coordinate dimension
Definition: Geometry.h:169
virtual int v_GetEid(int i) const
Definition: Geometry2D.cpp:230
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60