Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Geometry3D.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Geometry3D.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: 3D geometry information.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 
41 #include <SpatialDomains/SegGeom.h>
42 
43 #include <iomanip>
44 
45 using namespace std;
46 
47 namespace Nektar
48 {
49 namespace SpatialDomains
50 {
51  Geometry3D::Geometry3D()
52  {
53  }
54 
55  Geometry3D::Geometry3D(const int coordim):
56  Geometry(coordim)
57  {
58  ASSERTL0(m_coordim > 2,
59  "Coordinate dimension should be at least 3 for a 3D geometry.");
60  }
61 
63  {
64  }
65 
66 
67  //---------------------------------------
68  // Helper functions
69  //---------------------------------------
70 
71  /**
72  * @brief Return the ID of edge i in this element.
73  */
74  int Geometry3D::GetEid(int i) const
75  {
76  return v_GetEid(i);
77  }
78 
79  /**
80  * @brief Return face i in this element.
81  */
83  {
84  return v_GetFace(i);
85  }
86 
87  /**
88  * @brief Returns the element coordinate direction corresponding to a
89  * given face coordinate direction
90  */
91  int Geometry3D::GetDir(const int faceidx, const int facedir) const
92  {
93  return v_GetDir(faceidx, facedir);
94  }
95 
96  //---------------------------------------
97  // 3D Geometry Methods
98  //---------------------------------------
100  const Array<OneD, const NekDouble> &coords,
101  const Array<OneD, const NekDouble> &ptsx,
102  const Array<OneD, const NekDouble> &ptsy,
103  const Array<OneD, const NekDouble> &ptsz,
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 
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,zmap, F1,F2, F3;
122 
123  NekDouble derx_1, derx_2, derx_3, dery_1, dery_2, dery_3,
124  derz_1, derz_2, derz_3, jac;
125 
126  // save intiial guess for later reference if required.
127  NekDouble init0 = Lcoords[0], init1 = Lcoords[1], init2 = Lcoords[2];
128 
129  Array<OneD, NekDouble> DxD1(ptsx.num_elements());
130  Array<OneD, NekDouble> DxD2(ptsx.num_elements());
131  Array<OneD, NekDouble> DxD3(ptsx.num_elements());
132  Array<OneD, NekDouble> DyD1(ptsx.num_elements());
133  Array<OneD, NekDouble> DyD2(ptsx.num_elements());
134  Array<OneD, NekDouble> DyD3(ptsx.num_elements());
135  Array<OneD, NekDouble> DzD1(ptsx.num_elements());
136  Array<OneD, NekDouble> DzD2(ptsx.num_elements());
137  Array<OneD, NekDouble> DzD3(ptsx.num_elements());
138 
139  // Ideally this will be stored in m_geomfactors
140  m_xmap->PhysDeriv(ptsx,DxD1,DxD2,DxD3);
141  m_xmap->PhysDeriv(ptsy,DyD1,DyD2,DyD3);
142  m_xmap->PhysDeriv(ptsz,DzD1,DzD2,DzD3);
143 
144  int cnt=0;
146  Array<OneD, NekDouble> eta(3);
147 
148  F1 = F2 = F3 = 2000; // Starting value of Function
149 
150  while(cnt++ < MaxIterations)
151  {
152  // evaluate lagrange interpolant at Lcoords
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);
157 
158  //calculate the global point `corresponding to Lcoords
159  xmap = m_xmap->PhysEvaluate(I, ptsx);
160  ymap = m_xmap->PhysEvaluate(I, ptsy);
161  zmap = m_xmap->PhysEvaluate(I, ptsz);
162 
163  F1 = coords[0] - xmap;
164  F2 = coords[1] - ymap;
165  F3 = coords[2] - zmap;
166 
167  if(F1*F1 + F2*F2 + F3*F3 < ScaledTol)
168  {
169  resid = sqrt(F1*F1 + F2*F2 + F3*F3);
170  break;
171  }
172 
173  //Interpolate derivative metric at Lcoords
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);
183 
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);
187 
188  // use analytical inverse of derivitives which are also similar to
189  // those of metric factors.
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)
194  )/jac;
195 
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)
200  )/jac;
201 
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)
206  )/jac;
207 
208  if (fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv ||
209  fabs(Lcoords[0]) > LcoordDiv)
210  {
211  break; // lcoords have diverged so stop iteration
212  }
213  }
214 
215  resid = sqrt(F1*F1 + F2*F2 + F3*F3);
216 
217  if(cnt >= MaxIterations)
218  {
219  Array<OneD, NekDouble> collCoords(3);
220  m_xmap->LocCoordToLocCollapsed(Lcoords,collCoords);
221 
222  // if coordinate is inside element dump error!
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))
226  {
227  std::ostringstream ss;
228 
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) ;
233 
234  WARNINGL1(cnt < MaxIterations,ss.str());
235  }
236  }
237  }
238 
239  /**
240  * @brief Put all quadrature information into face/edge structure and
241  * backward transform.
242  *
243  * Note verts, edges, and faces are listed according to anticlockwise
244  * convention but points in _coeffs have to be in array format from left
245  * to right.
246  */
248  {
249  if (m_state == ePtsFilled)
250  return;
251 
252  int i,j,k;
253 
254  for(i = 0; i < m_forient.size(); i++)
255  {
256  m_faces[i]->FillGeom();
257 
258  int nFaceCoeffs = m_faces[i]->GetXmap()->GetNcoeffs();
259 
260  Array<OneD, unsigned int> mapArray (nFaceCoeffs);
261  Array<OneD, int> signArray(nFaceCoeffs);
262 
263  if (m_forient[i] < 9)
264  {
265  m_xmap->GetFaceToElementMap(
266  i, m_forient[i], mapArray, signArray,
267  m_faces[i]->GetXmap()->GetEdgeNcoeffs(0),
268  m_faces[i]->GetXmap()->GetEdgeNcoeffs(1));
269  }
270  else
271  {
272  m_xmap->GetFaceToElementMap(
273  i, m_forient[i], mapArray, signArray,
274  m_faces[i]->GetXmap()->GetEdgeNcoeffs(1),
275  m_faces[i]->GetXmap()->GetEdgeNcoeffs(0));
276  }
277 
278  for (j = 0; j < m_coordim; j++)
279  {
280  const Array<OneD, const NekDouble> &coeffs =
281  m_faces[i]->GetCoeffs(j);
282 
283  for (k = 0; k < nFaceCoeffs; k++)
284  {
285  NekDouble v = signArray[k] * coeffs[k];
286  m_coeffs[j][mapArray[k]] = v;
287  }
288  }
289  }
290 
292  }
293 
294  /**
295  * Generate the geometry factors for this element.
296  */
298  {
300  {
301  GeomType Gtype = eRegular;
302 
303  v_FillGeom();
304 
305  // check to see if expansions are linear
306  for(int i = 0; i < m_coordim; ++i)
307  {
308  if (m_xmap->GetBasisNumModes(0) != 2 ||
309  m_xmap->GetBasisNumModes(1) != 2 ||
310  m_xmap->GetBasisNumModes(2) != 2)
311  {
312  Gtype = eDeformed;
313  }
314  }
315 
317  Gtype, m_coordim, m_xmap, m_coeffs);
319  }
320  }
321 
322  /**
323  * @brief Given local collapsed coordinate Lcoord return the value of
324  * physical coordinate in direction i.
325  */
327  const int i, const Array<OneD, const NekDouble> &Lcoord)
328  {
330  "Geometry is not in physical space");
331 
332  Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
333  m_xmap->BwdTrans(m_coeffs[i], tmp);
334 
335  return m_xmap->PhysEvaluate(Lcoord, tmp);
336  }
337 
338 
339  //---------------------------------------
340  // Helper functions
341  //---------------------------------------
342 
343  /**
344  * @brief Return the dimension of this element.
345  */
347  {
348  return 3;
349  }
350 
351  /**
352  * @brief Return the vertex ID of vertex i.
353  */
354  int Geometry3D::v_GetVid(const int i) const
355  {
356  ASSERTL2(i >= 0 && i <= m_verts.size() - 1,
357  "Vertex ID must be between 0 and "+
358  boost::lexical_cast<string>(m_verts.size() - 1));
359  return m_verts[i]->GetVid();
360  }
361 
362  /**
363  * @brief Return vertex i in this element.
364  */
366  {
367  return m_verts[i];
368  }
369 
370  /**
371  * @brief Return edge i of this element.
372  */
374  {
375  ASSERTL2(i >= 0 && i <= m_edges.size() - 1,
376  "Edge ID must be between 0 and "+
377  boost::lexical_cast<string>(m_edges.size() - 1));
378  return m_edges[i];
379  }
380 
381  /**
382  * @brief Return the orientation of edge i in this element.
383  */
385  const int i) const
386  {
387  ASSERTL2(i >= 0 && i <= m_edges.size() - 1,
388  "Edge ID must be between 0 and "+
389  boost::lexical_cast<string>(m_edges.size() - 1));
390  return m_eorient[i];
391  }
392 
393  /**
394  * @brief Return the ID of edge i in this element.
395  */
396  int Geometry3D::v_GetEid(int i) const
397  {
398  ASSERTL2(i >= 0 && i <= m_edges.size() - 1,
399  "Edge ID must be between 0 and "+
400  boost::lexical_cast<string>(m_edges.size() - 1));
401  return m_edges[i]->GetEid();
402  }
403 
404  /**
405  * @brief Return face i in this element.
406  */
408  {
409  ASSERTL2((i >=0) && (i <= 5),"Edge id must be between 0 and 4");
410  return m_faces[i];
411  }
412 
413  /**
414  * @brief Return the orientation of face i in this element.
415  */
417  {
418  ASSERTL2(i >= 0 && i <= m_faces.size() - 1,
419  "Face ID must be between 0 and "+
420  boost::lexical_cast<string>(m_faces.size() - 1));
421  return m_forient[i];
422  }
423 
424  /**
425  * @brief Return the ID of face i in this element.
426  */
427  int Geometry3D::v_GetFid(int i) const
428  {
429  ASSERTL2(i >= 0 && i <= m_faces.size() - 1,
430  "Face ID must be between 0 and "+
431  boost::lexical_cast<string>(m_faces.size() - 1));
432  return m_faces[i]->GetFid();
433  }
434 
435  /**
436  * @brief Return the ID of this element.
437  */
439  {
440  return m_eid;
441  }
442 
443  /**
444  * @brief Return the j-th basis of the i-th co-ordinate dimension.
445  */
447  const int i)
448  {
449  return m_xmap->GetBasis(i);
450  }
451 
452  /**
453  * @brief Return the local ID of a given edge.
454  *
455  * The local ID of an edge is a number between 0 and the number of edges
456  * in this element. If the edge is not found, this function returns -1.
457  */
459  {
460  int returnval = -1;
461 
462  SegGeomVector::iterator edgeIter;
463  int i;
464 
465  for (i=0,edgeIter = m_edges.begin(); edgeIter != m_edges.end(); ++edgeIter,++i)
466  {
467  if (*edgeIter == edge)
468  {
469  returnval = i;
470  break;
471  }
472  }
473 
474  return returnval;
475  }
476 
477  /**
478  * @brief Return the local ID of a given face.
479  *
480  * The local ID of a face is a number between 0 and the number of faces
481  * in this element. If the face is not found, this function returns -1.
482  */
484  {
485  int i = 0;
486 
488  for (i = 0, f = m_faces.begin(); f != m_faces.end(); ++f,++i)
489  {
490  if (*f == face)
491  {
492  break;
493  }
494  }
495  return i;
496  }
497 
498  //---------------------------------------
499  // Element connection functions
500  //---------------------------------------
501 
502  void Geometry3D::v_AddElmtConnected(int gvo_id, int locid)
503  {
504  CompToElmt ee(gvo_id,locid);
505  m_elmtmap.push_back(ee);
506  }
507 
509  {
510  return int(m_elmtmap.size());
511  }
512 
513  bool Geometry3D::v_IsElmtConnected(int gvo_id, int locid) const
514  {
515  std::list<CompToElmt>::const_iterator def;
516  CompToElmt ee(gvo_id,locid);
517 
518  def = find(m_elmtmap.begin(),m_elmtmap.end(),ee);
519 
520  // Found the element connectivity object in the list
521  return (def != m_elmtmap.end());
522  }
523 
525  {
526  m_owndata = true;
527  }
528 
529 }; //end of namespace
530 }; //end of namespace
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry3D.h:92
Base class for shape geometry information.
Definition: Geometry.h:76
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...
Definition: Geometry3D.cpp:91
GeomFactorsSharedPtr m_geomFactors
Definition: Geometry.h:170
virtual PointGeomSharedPtr v_GetVertex(int i) const
Return vertex i in this element.
Definition: Geometry3D.cpp:365
std::list< CompToElmt > m_elmtmap
Definition: Geometry3D.h:94
virtual int v_GetShapeDim() const
Return the dimension of this element.
Definition: Geometry3D.cpp:346
Structure holding graphvertexobject id and local element facet id.
int GetEid(int i) const
Return the ID of edge i in this element.
Definition: Geometry3D.cpp:74
Geometry2DSharedPtr GetFace(int i)
Return face i in this element.
Definition: Geometry3D.cpp:82
STL namespace.
StdRegions::StdExpansionSharedPtr GetXmap() const
Definition: Geometry.h:383
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:93
virtual int v_WhichEdge(SegGeomSharedPtr edge)
Return the local ID of a given edge.
Definition: Geometry3D.cpp:458
virtual int v_GetFid(int i) const
Return the ID of face i in this element.
Definition: Geometry3D.cpp:427
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)
Definition: Geometry3D.cpp:99
virtual int v_WhichFace(Geometry2DSharedPtr face)
Return the local ID of a given face.
Definition: Geometry3D.cpp:483
virtual void v_AddElmtConnected(int gvo_id, int locid)
Definition: Geometry3D.cpp:502
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
virtual bool v_IsElmtConnected(int gvo_id, int locid) const
Definition: Geometry3D.cpp:513
virtual int v_GetDir(const int faceidx, const int facedir) const =0
double NekDouble
virtual int v_GetVid(int i) const
Return the vertex ID of vertex i.
Definition: Geometry3D.cpp:354
virtual StdRegions::Orientation v_GetEorient(const int i) const
Return the orientation of edge i in this element.
Definition: Geometry3D.cpp:384
virtual const LibUtilities::BasisSharedPtr v_GetBasis(const int i)
Return the j-th basis of the i-th co-ordinate dimension.
Definition: Geometry3D.cpp:446
virtual int v_GetEid() const
Return the ID of this element.
Definition: Geometry3D.cpp:438
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
Definition: Geometry3D.cpp:247
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
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.
Definition: Geometry3D.cpp:416
virtual const Geometry2DSharedPtr v_GetFace(int i) const
Return face i in this element.
Definition: Geometry3D.cpp:407
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Definition: Geometry.h:180
#define WARNINGL1(condition, msg)
Definition: ErrorUtil.hpp:192
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
Geometry is straight-sided with constant geometric factors.
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:315
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...
Definition: Geometry3D.cpp:326
virtual const SegGeomSharedPtr v_GetEdge(int i) const
Return edge i of this element.
Definition: Geometry3D.cpp:373
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
Definition: Geometry.h:175
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:723
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
Geometry is curved or has non-constant factors.
int m_coordim
coordinate dimension
Definition: Geometry.h:169
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
virtual int v_NumElmtConnected() const
Definition: Geometry3D.cpp:508