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