Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 namespace Nektar
42 {
43  namespace SpatialDomains
44  {
46  {
47  }
48 
49  Geometry2D::Geometry2D(const int coordim):
50  Geometry(coordim)
51  {
52  ASSERTL0(m_coordim > 1,
53  "Coordinate dimension should be at least 2 for a 2D geometry");
54  }
55 
57  {
58  }
59 
60  int Geometry2D::GetFid() const
61  {
62  return v_GetFid();
63  }
64 
66  {
67  return v_GetEdgeBasis(i);
68  }
69 
71  {
72  return v_GetFace(i);
73  }
74 
76  {
77  return v_GetFaceOrient(i);
78  }
79 
81  {
82  return v_GetEdge(i);
83  }
84 
86  {
87  return v_GetCartesianEorient(i);
88  }
89 
91  {
92  return v_WhichEdge(edge);
93  }
94 
96  {
97  return v_WhichFace(face);
98  }
99 
101  const Array<OneD, const NekDouble> &coords,
102  const Array<OneD, const NekDouble> &ptsx,
103  const Array<OneD, const NekDouble> &ptsy,
104  Array<OneD, NekDouble> &Lcoords,
105  NekDouble &resid)
106  {
107  // Maximum iterations for convergence
108  const int MaxIterations = 51;
109  // |x-xp|^2 < EPSILON error tolerance
110  const NekDouble Tol = 1.e-8;
111  // |r,s| > LcoordDIV stop the search
112  const NekDouble LcoordDiv = 15.0;
113 
114  Array<OneD, const NekDouble > Jac =
115  m_geomFactors->GetJac(m_xmap->GetPointsKeys());
116 
117  NekDouble ScaledTol = Vmath::Vsum(Jac.num_elements(),Jac,1)/
118  ((NekDouble)Jac.num_elements());
119  ScaledTol *= Tol;
120 
121  NekDouble xmap,ymap, F1,F2;
122  NekDouble derx_1, derx_2, dery_1, dery_2,jac;
123 
124  // save intiial guess for later reference if required.
125  NekDouble init0 = Lcoords[0], init1 = Lcoords[1];
126 
127  Array<OneD, NekDouble> DxD1(ptsx.num_elements());
128  Array<OneD, NekDouble> DxD2(ptsx.num_elements());
129  Array<OneD, NekDouble> DyD1(ptsx.num_elements());
130  Array<OneD, NekDouble> DyD2(ptsx.num_elements());
131 
132  // Ideally this will be stored in m_geomfactors
133  m_xmap->PhysDeriv(ptsx,DxD1,DxD2);
134  m_xmap->PhysDeriv(ptsy,DyD1,DyD2);
135 
136  int cnt=0;
137  Array<OneD, DNekMatSharedPtr > I(2);
138  Array<OneD, NekDouble> eta(2);
139 
140  F1 = F2 = 2000; // Starting value of Function
141 
142  while(cnt++ < MaxIterations)
143  {
144  // evaluate lagrange interpolant at Lcoords
145  m_xmap->LocCoordToLocCollapsed(Lcoords,eta);
146  I[0] = m_xmap->GetBasis(0)->GetI(eta);
147  I[1] = m_xmap->GetBasis(1)->GetI(eta+1);
148 
149  //calculate the global point `corresponding to Lcoords
150  xmap = m_xmap->PhysEvaluate(I, ptsx);
151  ymap = m_xmap->PhysEvaluate(I, ptsy);
152 
153  F1 = coords[0] - xmap;
154  F2 = coords[1] - ymap;
155 
156  if(F1*F1 + F2*F2 < ScaledTol)
157  {
158  resid = sqrt(F1*F1 + F2*F2);
159  break;
160  }
161 
162  //Interpolate derivative metric at Lcoords
163  derx_1 = m_xmap->PhysEvaluate(I, DxD1);
164  derx_2 = m_xmap->PhysEvaluate(I, DxD2);
165  dery_1 = m_xmap->PhysEvaluate(I, DyD1);
166  dery_2 = m_xmap->PhysEvaluate(I, DyD2);
167 
168  jac = dery_2*derx_1 - dery_1*derx_2;
169 
170  // use analytical inverse of derivitives which are
171  // also similar to those of metric factors.
172  Lcoords[0] = Lcoords[0] + (dery_2*(coords[0]-xmap) -
173  derx_2*(coords[1]-ymap))/jac;
174 
175  Lcoords[1] = Lcoords[1] + ( - dery_1*(coords[0]-xmap)
176  + derx_1*(coords[1]-ymap))/jac;
177 
178  if(fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv)
179  {
180  break; // lcoords have diverged so stop iteration
181  }
182  }
183 
184  resid = sqrt(F1*F1 + F2*F2);
185 
186  if(cnt >= MaxIterations)
187  {
188  Array<OneD, NekDouble> collCoords(2);
189  m_xmap->LocCoordToLocCollapsed(Lcoords,collCoords);
190 
191  // if coordinate is inside element dump error!
192  if((collCoords[0] >= -1.0 && collCoords[0] <= 1.0)&&
193  (collCoords[1] >= -1.0 && collCoords[1] <= 1.0))
194  {
195  std::ostringstream ss;
196 
197  ss << "Reached MaxIterations (" << MaxIterations
198  << ") in Newton iteration ";
199  ss << "Init value ("<< setprecision(4) << init0 << ","
200  << init1<< "," <<") ";
201  ss << "Fin value ("<<Lcoords[0] << "," << Lcoords[1]
202  << "," << ") ";
203  ss << "Resid = " << resid << " Tolerance = "
204  << sqrt(ScaledTol) ;
205 
206  WARNINGL1(cnt < MaxIterations,ss.str());
207  }
208  }
209  }
210 
211  int Geometry2D::v_GetFid() const
212  {
214  "This function is only valid for shape type geometries");
215  return 0;
216  }
217 
219  {
221  "This function is only valid for shape type geometries");
223  return returnval;
224  }
225 
226 
227  int Geometry2D::v_GetEid(int i) const
228  {
230  "This function is only valid for shape type geometries");
231  return 0;
232  }
233 
235  {
237  "This function is only valid for shape type geometries");
238  SegGeomSharedPtr returnval;
239  return returnval;
240  }
241 
243  {
245  "This function is only valid for shape type geometries");
246  PointGeomSharedPtr returnval;
247  return returnval;
248  }
249 
251  {
253  "This function is only valid for shape type geometries");
254  return StdRegions::eForwards;
255  }
256 
258  {
260  "This function is only valid for shape type geometries");
261  Geometry2DSharedPtr returnval;
262  return returnval;
263  }
264 
266  {
268  "This function is only valid for shape type geometries");
270  }
271 
273  {
275  "This function is only valid for shape type geometries");
276  return StdRegions::eForwards;
277  }
278 
280  {
282  "This function is only valid for shape type geometries");
283  return 0;
284  }
285 
287  {
289  "This function is only valid for shape type geometries");
290  return 0;
291  }
292 
294  {
295  return 2;
296  }
297 
299  const Array<OneD, const NekDouble>& gloCoord,
300  NekDouble tol)
301  {
303  "This function has not been defined for this geometry");
304  return false;
305  }
306 
307 
308  }; //end of namespace
309 }; //end of namespace
310