Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HexGeom.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: HexGeom.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: Hexahedral geometry definition.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <SpatialDomains/HexGeom.h>
38 #include <StdRegions/StdHexExp.h>
39 #include <SpatialDomains/SegGeom.h>
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46  namespace SpatialDomains
47  {
48  const unsigned int HexGeom::VertexEdgeConnectivity[8][3] = {
49  {0,3,4},{0,1,5},{1,2,6},{2,3,7},
50  {4,8,11},{5,8,9},{6,9,10},{7,10,11}};
51  const unsigned int HexGeom::VertexFaceConnectivity[8][3] = {
52  {0,1,4},{0,1,2},{0,2,3},{0,3,4},
53  {1,4,5},{1,2,5},{2,3,5},{3,4,5}};
54  const unsigned int HexGeom::EdgeFaceConnectivity[12][2] = {
55  {0,1},{0,2},{0,3},{0,4},{1,4},{1,2},{2,3},{3,4},
56  {1,5},{2,5},{3,5},{4,5}};
57 
58  HexGeom::HexGeom()
59  {
60  m_shapeType = LibUtilities::eHexahedron;
61  }
62 
63  HexGeom::HexGeom(const QuadGeomSharedPtr faces[]):
64  Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
65  {
67 
68  /// Copy the face shared pointers
69  m_faces.insert(m_faces.begin(), faces, faces+HexGeom::kNfaces);
70 
71  /// Set up orientation vectors with correct amount of elements.
72  m_eorient.resize(kNedges);
73  m_forient.resize(kNfaces);
74 
79  SetUpXmap();
80  SetUpCoeffs(m_xmap->GetNcoeffs());
81  }
82 
84  {
85 
86  }
87 
89  {
91  {
92  int i,f;
93  GeomType Gtype = eRegular;
94 
95  v_FillGeom();
96 
97  // check to see if expansions are linear
98  for(i = 0; i < m_coordim; ++i)
99  {
100  if (m_xmap->GetBasisNumModes(0) != 2 ||
101  m_xmap->GetBasisNumModes(1) != 2 ||
102  m_xmap->GetBasisNumModes(2) != 2 )
103  {
104  Gtype = eDeformed;
105  }
106  }
107 
108  // check to see if all faces are parallelograms
109  if(Gtype == eRegular)
110  {
111  const unsigned int faceVerts[kNfaces][QuadGeom::kNverts] =
112  { {0,1,2,3} ,
113  {0,1,5,4} ,
114  {1,2,6,5} ,
115  {3,2,6,7} ,
116  {0,3,7,4} ,
117  {4,5,6,7} };
118 
119  for(f = 0; f < kNfaces; f++)
120  {
121  // Ensure each face is a parallelogram? Check this.
122  for (i = 0; i < m_coordim; i++)
123  {
124  if( fabs( (*m_verts[ faceVerts[f][0] ])(i) -
125  (*m_verts[ faceVerts[f][1] ])(i) +
126  (*m_verts[ faceVerts[f][2] ])(i) -
127  (*m_verts[ faceVerts[f][3] ])(i) )
129  {
130  Gtype = eDeformed;
131  break;
132  }
133  }
134 
135  if (Gtype == eDeformed)
136  {
137  break;
138  }
139  }
140  }
141 
143  Gtype, m_coordim, m_xmap, m_coeffs);
145  }
146  }
147 
149  const Array<OneD, const NekDouble> &coords,
150  Array<OneD, NekDouble> &Lcoords)
151  {
152  NekDouble ptdist = 1e6;
153  int i;
154 
155  v_FillGeom();
156 
157  // calculate local coordinate for coord
158  if(GetMetricInfo()->GetGtype() == eRegular)
159  {
160  NekDouble len0 = 0.0 ;
161  NekDouble len1 = 0.0;
162  NekDouble len2 = 0.0;
163  NekDouble xi0 = 0.0;
164  NekDouble xi1 = 0.0;
165  NekDouble xi2 = 0.0;
166  Array<OneD, NekDouble> pts(m_xmap->GetTotPoints());
167  int nq0, nq1, nq2;
168 
169  // get points;
170  //find end points
171  for(i = 0; i < m_coordim; ++i)
172  {
173  nq0 = m_xmap->GetNumPoints(0);
174  nq1 = m_xmap->GetNumPoints(1);
175  nq2 = m_xmap->GetNumPoints(2);
176 
177  m_xmap->BwdTrans(m_coeffs[i], pts);
178 
179  // use projection to side 1 to determine xi_1 coordinate based on length
180  len0 += (pts[nq0-1]-pts[0])*(pts[nq0-1]-pts[0]);
181  xi0 += (coords[i] -pts[0])*(pts[nq0-1]-pts[0]);
182 
183  // use projection to side 4 to determine xi_2 coordinate based on length
184  len1 += (pts[nq0*(nq1-1)]-pts[0])*(pts[nq0*(nq1-1)]-pts[0]);
185  xi1 += (coords[i] -pts[0])*(pts[nq0*(nq1-1)]-pts[0]);
186 
187  // use projection to side 4 to determine xi_3 coordinate based on length
188  len2 += (pts[nq0*nq1*(nq2-1)]-pts[0])*(pts[nq0*nq1*(nq2-1)]-pts[0]);
189  xi2 += (coords[i] -pts[0])*(pts[nq0*nq1*(nq2-1)]-pts[0]);
190  }
191 
192  Lcoords[0] = 2*xi0/len0-1.0;
193  Lcoords[1] = 2*xi1/len1-1.0;
194  Lcoords[2] = 2*xi2/len2-1.0;
195 
196  // Set ptdist to distance to nearest vertex
197  // Point inside tetrahedron
198  PointGeom r(m_coordim, 0, coords[0], coords[1], coords[2]);
199  for(int i = 0; i < 8; ++i)
200  {
201  ptdist = min(ptdist,r.dist(*m_verts[i]));
202  }
203  }
204  else
205  {
206  // Determine nearest point of coords to values in m_xmap
207  int npts = m_xmap->GetTotPoints();
208  Array<OneD, NekDouble> ptsx(npts), ptsy(npts), ptsz(npts);
209  Array<OneD, NekDouble> tmp1(npts), tmp2(npts);
210 
211  m_xmap->BwdTrans(m_coeffs[0], ptsx);
212  m_xmap->BwdTrans(m_coeffs[1], ptsy);
213  m_xmap->BwdTrans(m_coeffs[2], ptsz);
214 
215  const Array<OneD, const NekDouble> za = m_xmap->GetPoints(0);
216  const Array<OneD, const NekDouble> zb = m_xmap->GetPoints(1);
217  const Array<OneD, const NekDouble> zc = m_xmap->GetPoints(2);
218 
219  //guess the first local coords based on nearest point
220  Vmath::Sadd(npts, -coords[0], ptsx,1,tmp1,1);
221  Vmath::Vmul (npts, tmp1,1,tmp1,1,tmp1,1);
222  Vmath::Sadd(npts, -coords[1], ptsy,1,tmp2,1);
223  Vmath::Vvtvp(npts, tmp2,1,tmp2,1,tmp1,1,tmp1,1);
224  Vmath::Sadd(npts, -coords[2], ptsz,1,tmp2,1);
225  Vmath::Vvtvp(npts, tmp2,1,tmp2,1,tmp1,1,tmp1,1);
226 
227  int min_i = Vmath::Imin(npts,tmp1,1);
228 
229  // distance from coordinate to nearest point for return value.
230  ptdist = sqrt(tmp1[min_i]);
231 
232  // Get Local coordinates
233  int qa = za.num_elements(), qb = zb.num_elements();
234  Lcoords[2] = zc[min_i/(qa*qb)];
235  min_i = min_i%(qa*qb);
236  Lcoords[1] = zb[min_i/qa];
237  Lcoords[0] = za[min_i%qa];
238 
239  // Perform newton iteration to find local coordinates
240  NekDouble resid = 0.0;
241  NewtonIterationForLocCoord(coords, ptsx, ptsy, ptsz, Lcoords,
242  resid);
243  }
244  return ptdist;
245  }
246 
247  /**
248  * @brief Determines if a point specified in global coordinates is
249  * located within this tetrahedral geometry.
250  */
252  const Array<OneD, const NekDouble> &gloCoord, NekDouble tol)
253  {
254  Array<OneD,NekDouble> locCoord(GetCoordim(),0.0);
255  return v_ContainsPoint(gloCoord,locCoord,tol);
256  }
257 
258 
260  const Array<OneD, const NekDouble> &gloCoord,
261  Array<OneD, NekDouble> &locCoord,
262  NekDouble tol)
263  {
264  NekDouble resid;
265  return v_ContainsPoint(gloCoord,locCoord,tol,resid);
266  }
267 
269  const Array<OneD, const NekDouble> &gloCoord,
270  Array<OneD, NekDouble> &locCoord,
271  NekDouble tol,
272  NekDouble &resid)
273  {
274  ASSERTL1(gloCoord.num_elements() == 3,
275  "Three dimensional geometry expects three coordinates.");
276 
277  // find min, max point and check if within twice this
278  // distance other false this is advisable since
279  // GetLocCoord is expensive for non regular elements.
280  if(GetMetricInfo()->GetGtype() != eRegular)
281  {
282  int i;
283  Array<OneD, NekDouble> mincoord(3), maxcoord(3);
284  NekDouble diff = 0.0;
285 
286  v_FillGeom();
287 
288  const int npts = m_xmap->GetTotPoints();
289  Array<OneD, NekDouble> pts(npts);
290 
291  for(i = 0; i < 3; ++i)
292  {
293  m_xmap->BwdTrans(m_coeffs[i], pts);
294 
295  mincoord[i] = Vmath::Vmin(pts.num_elements(),pts,1);
296  maxcoord[i] = Vmath::Vmax(pts.num_elements(),pts,1);
297 
298  diff = max(maxcoord[i] - mincoord[i],diff);
299  }
300 
301  for(i = 0; i < 3; ++i)
302  {
303  if((gloCoord[i] < mincoord[i] - 0.2*diff)||
304  (gloCoord[i] > maxcoord[i] + 0.2*diff))
305  {
306  return false;
307  }
308  }
309  }
310 
311  v_GetLocCoords(gloCoord, locCoord);
312 
313  if (locCoord[0] >= -(1+tol) && locCoord[0] <= 1+tol
314  && locCoord[1] >= -(1+tol) && locCoord[1] <= 1+tol
315  && locCoord[2] >= -(1+tol) && locCoord[2] <= 1+tol)
316  {
317  return true;
318  }
319 
320  // If out of range clamp locCoord to be within [-1,1]^3
321  // since any larger value will be very oscillatory if
322  // called by 'returnNearestElmt' option in
323  // ExpList::GetExpIndex
324  for(int i = 0; i < 3; ++i)
325  {
326  if(locCoord[i] <-(1+tol))
327  {
328  locCoord[i] = -(1+tol);
329  }
330 
331  if(locCoord[i] > (1+tol))
332  {
333  locCoord[i] = 1+tol;
334  }
335  }
336 
337  return false;
338  }
339 
341  {
342  return 8;
343  }
344 
346  {
347  return 12;
348  }
349 
351  {
352  return 6;
353  }
354 
355  int HexGeom::v_GetVertexEdgeMap(const int i, const int j) const
356  {
357  return VertexEdgeConnectivity[i][j];
358  }
359 
360  int HexGeom::v_GetVertexFaceMap(const int i, const int j) const
361  {
362  return VertexFaceConnectivity[i][j];
363  }
364 
365  int HexGeom::v_GetEdgeFaceMap(const int i, const int j) const
366  {
367  return EdgeFaceConnectivity[i][j];
368  }
369 
370  int HexGeom::v_GetDir(const int faceidx, const int facedir) const
371  {
372  if (faceidx == 0 || faceidx == 5)
373  {
374  return facedir;
375  }
376  else if (faceidx == 1 || faceidx == 3)
377  {
378  return 2 * facedir;
379  }
380  else
381  {
382  return 1 + facedir;
383  }
384  }
385 
387  {
388  // find edge 0
389  int i,j;
390  unsigned int check;
391 
392  SegGeomSharedPtr edge;
393 
394  // First set up the 4 bottom edges
395  int f;
396  for(f = 1; f < 5 ; f++)
397  {
398  check = 0;
399  for(i = 0; i < 4; i++)
400  {
401  for(j = 0; j < 4; j++)
402  {
403  if( (m_faces[0])->GetEid(i) == (m_faces[f])->GetEid(j) )
404  {
405  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[0])->GetEdge(i));
406  m_edges.push_back(edge);
407  check++;
408  }
409  }
410  }
411 
412  if( check < 1 )
413  {
414  std::ostringstream errstrm;
415  errstrm << "Connected faces do not share an edge. Faces ";
416  errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[f])->GetFid();
417  ASSERTL0(false, errstrm.str());
418  }
419  else if( check > 1)
420  {
421  std::ostringstream errstrm;
422  errstrm << "Connected faces share more than one edge. Faces ";
423  errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[f])->GetFid();
424  ASSERTL0(false, errstrm.str());
425  }
426  }
427 
428  // Then, set up the 4 vertical edges
429  check = 0;
430  for(i = 0; i < 4; i++)
431  {
432  for(j = 0; j < 4; j++)
433  {
434  if( (m_faces[1])->GetEid(i) == (m_faces[4])->GetEid(j) )
435  {
436  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[1])->GetEdge(i));
437  m_edges.push_back(edge);
438  check++;
439  }
440  }
441  }
442  if( check < 1 )
443  {
444  std::ostringstream errstrm;
445  errstrm << "Connected faces do not share an edge. Faces ";
446  errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[4])->GetFid();
447  ASSERTL0(false, errstrm.str());
448  }
449  else if( check > 1)
450  {
451  std::ostringstream errstrm;
452  errstrm << "Connected faces share more than one edge. Faces ";
453  errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[4])->GetFid();
454  ASSERTL0(false, errstrm.str());
455  }
456  for(f = 1; f < 4 ; f++)
457  {
458  check = 0;
459  for(i = 0; i < 4; i++)
460  {
461  for(j = 0; j < 4; j++)
462  {
463  if( (m_faces[f])->GetEid(i) == (m_faces[f+1])->GetEid(j) )
464  {
465  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[f])->GetEdge(i));
466  m_edges.push_back(edge);
467  check++;
468  }
469  }
470  }
471 
472  if( check < 1 )
473  {
474  std::ostringstream errstrm;
475  errstrm << "Connected faces do not share an edge. Faces ";
476  errstrm << (m_faces[f])->GetFid() << ", " << (m_faces[f+1])->GetFid();
477  ASSERTL0(false, errstrm.str());
478  }
479  else if( check > 1)
480  {
481  std::ostringstream errstrm;
482  errstrm << "Connected faces share more than one edge. Faces ";
483  errstrm << (m_faces[f])->GetFid() << ", " << (m_faces[f+1])->GetFid();
484  ASSERTL0(false, errstrm.str());
485  }
486  }
487 
488  // Finally, set up the 4 top vertices
489  for(f = 1; f < 5 ; f++)
490  {
491  check = 0;
492  for(i = 0; i < 4; i++)
493  {
494  for(j = 0; j < 4; j++)
495  {
496  if( (m_faces[5])->GetEid(i) == (m_faces[f])->GetEid(j) )
497  {
498  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[5])->GetEdge(i));
499  m_edges.push_back(edge);
500  check++;
501  }
502  }
503  }
504 
505  if( check < 1 )
506  {
507  std::ostringstream errstrm;
508  errstrm << "Connected faces do not share an edge. Faces ";
509  errstrm << (m_faces[5])->GetFid() << ", " << (m_faces[f])->GetFid();
510  ASSERTL0(false, errstrm.str());
511  }
512  else if( check > 1)
513  {
514  std::ostringstream errstrm;
515  errstrm << "Connected faces share more than one edge. Faces ";
516  errstrm << (m_faces[5])->GetFid() << ", " << (m_faces[f])->GetFid();
517  ASSERTL0(false, errstrm.str());
518  }
519  }
520  }
521 
523  {
524  // Set up the first 2 vertices (i.e. vertex 0,1)
525  if( ( m_edges[0]->GetVid(0) == m_edges[1]->GetVid(0) ) ||
526  ( m_edges[0]->GetVid(0) == m_edges[1]->GetVid(1) ) )
527  {
528  m_verts.push_back(m_edges[0]->GetVertex(1));
529  m_verts.push_back(m_edges[0]->GetVertex(0));
530  }
531  else if( ( m_edges[0]->GetVid(1) == m_edges[1]->GetVid(0) ) ||
532  ( m_edges[0]->GetVid(1) == m_edges[1]->GetVid(1) ) )
533  {
534  m_verts.push_back(m_edges[0]->GetVertex(0));
535  m_verts.push_back(m_edges[0]->GetVertex(1));
536  }
537  else
538  {
539  std::ostringstream errstrm;
540  errstrm << "Connected edges do not share a vertex. Edges ";
541  errstrm << m_edges[0]->GetEid() << ", " << m_edges[1]->GetEid();
542  ASSERTL0(false, errstrm.str());
543  }
544 
545  // set up the other bottom vertices (i.e. vertex 2,3)
546  int i;
547  for(i = 1; i < 3; i++)
548  {
549  if( m_edges[i]->GetVid(0) == m_verts[i]->GetVid() )
550  {
551  m_verts.push_back(m_edges[i]->GetVertex(1));
552  }
553  else if( m_edges[i]->GetVid(1) == m_verts[i]->GetVid() )
554  {
555  m_verts.push_back(m_edges[i]->GetVertex(0));
556  }
557  else
558  {
559  std::ostringstream errstrm;
560  errstrm << "Connected edges do not share a vertex. Edges ";
561  errstrm << m_edges[i]->GetEid() << ", " << m_edges[i-1]->GetEid();
562  ASSERTL0(false, errstrm.str());
563  }
564  }
565 
566  // set up top vertices
567  // First, set up vertices 4,5
568  if( ( m_edges[8]->GetVid(0) == m_edges[9]->GetVid(0) ) ||
569  ( m_edges[8]->GetVid(0) == m_edges[9]->GetVid(1) ) )
570  {
571  m_verts.push_back(m_edges[8]->GetVertex(1));
572  m_verts.push_back(m_edges[8]->GetVertex(0));
573  }
574  else if( ( m_edges[8]->GetVid(1) == m_edges[9]->GetVid(0) ) ||
575  ( m_edges[8]->GetVid(1) == m_edges[9]->GetVid(1) ) )
576  {
577  m_verts.push_back(m_edges[8]->GetVertex(0));
578  m_verts.push_back(m_edges[8]->GetVertex(1));
579  }
580  else
581  {
582  std::ostringstream errstrm;
583  errstrm << "Connected edges do not share a vertex. Edges ";
584  errstrm << m_edges[8]->GetEid() << ", " << m_edges[9]->GetEid();
585  ASSERTL0(false, errstrm.str());
586  }
587 
588  // set up the other top vertices (i.e. vertex 6,7)
589  for(i = 9; i < 11; i++)
590  {
591  if( m_edges[i]->GetVid(0) == m_verts[i-4]->GetVid() )
592  {
593  m_verts.push_back(m_edges[i]->GetVertex(1));
594  }
595  else if( m_edges[i]->GetVid(1) == m_verts[i-4]->GetVid() )
596  {
597  m_verts.push_back(m_edges[i]->GetVertex(0));
598  }
599  else
600  {
601  std::ostringstream errstrm;
602  errstrm << "Connected edges do not share a vertex. Edges ";
603  errstrm << m_edges[i]->GetEid() << ", " << m_edges[i-1]->GetEid();
604  ASSERTL0(false, errstrm.str());
605  }
606  }
607  }
608 
610  {
611  int f,i;
612 
613  // These arrays represent the vector of the A and B
614  // coordinate of the local elemental coordinate system
615  // where A corresponds with the coordinate direction xi_i
616  // with the lowest index i (for that particular face)
617  // Coordinate 'B' then corresponds to the other local
618  // coordinate (i.e. with the highest index)
619  Array<OneD,NekDouble> elementAaxis(m_coordim);
620  Array<OneD,NekDouble> elementBaxis(m_coordim);
621 
622  // These arrays correspond to the local coordinate
623  // system of the face itself (i.e. the Geometry2D)
624  // faceAaxis correspond to the xi_0 axis
625  // faceBaxis correspond to the xi_1 axis
626  Array<OneD,NekDouble> faceAaxis(m_coordim);
627  Array<OneD,NekDouble> faceBaxis(m_coordim);
628 
629  // This is the base vertex of the face (i.e. the Geometry2D)
630  // This corresponds to thevertex with local ID 0 of the
631  // Geometry2D
632  unsigned int baseVertex;
633 
634  // The lenght of the vectors above
635  NekDouble elementAaxis_length;
636  NekDouble elementBaxis_length;
637  NekDouble faceAaxis_length;
638  NekDouble faceBaxis_length;
639 
640  // This 2D array holds the local id's of all the vertices
641  // for every face. For every face, they are ordered in such
642  // a way that the implementation below allows a unified approach
643  // for all faces.
644  const unsigned int faceVerts[kNfaces][QuadGeom::kNverts] =
645  { {0,1,2,3} ,
646  {0,1,5,4} ,
647  {1,2,6,5} ,
648  {3,2,6,7} ,
649  {0,3,7,4} ,
650  {4,5,6,7} };
651 
652  NekDouble dotproduct1 = 0.0;
653  NekDouble dotproduct2 = 0.0;
654 
655  unsigned int orientation;
656 
657  // Loop over all the faces to set up the orientation
658  for(f = 0; f < kNqfaces + kNtfaces; f++)
659  {
660  // initialisation
661  elementAaxis_length = 0.0;
662  elementBaxis_length = 0.0;
663  faceAaxis_length = 0.0;
664  faceBaxis_length = 0.0;
665 
666  dotproduct1 = 0.0;
667  dotproduct2 = 0.0;
668 
669  baseVertex = m_faces[f]->GetVid(0);
670 
671  // We are going to construct the vectors representing the A and B axis
672  // of every face. These vectors will be constructed as a vector-representation
673  // of the edges of the face. However, for both coordinate directions, we can
674  // represent the vectors by two different edges. That's why we need to make sure that
675  // we pick the edge to which the baseVertex of the Geometry2D-representation of the face
676  // belongs...
677  if( baseVertex == m_verts[ faceVerts[f][0] ]->GetVid() )
678  {
679  for(i = 0; i < m_coordim; i++)
680  {
681  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
682  elementBaxis[i] = (*m_verts[ faceVerts[f][3] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
683  }
684  }
685  else if( baseVertex == m_verts[ faceVerts[f][1] ]->GetVid() )
686  {
687  for(i = 0; i < m_coordim; i++)
688  {
689  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
690  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
691  }
692  }
693  else if( baseVertex == m_verts[ faceVerts[f][2] ]->GetVid() )
694  {
695  for(i = 0; i < m_coordim; i++)
696  {
697  elementAaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][3] ])[i];
698  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
699  }
700  }
701  else if( baseVertex == m_verts[ faceVerts[f][3] ]->GetVid() )
702  {
703  for(i = 0; i < m_coordim; i++)
704  {
705  elementAaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][3] ])[i];
706  elementBaxis[i] = (*m_verts[ faceVerts[f][3] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
707  }
708  }
709  else
710  {
711  ASSERTL0(false, "Could not find matching vertex for the face");
712  }
713 
714  // Now, construct the edge-vectors of the local coordinates of
715  // the Geometry2D-representation of the face
716  for(i = 0; i < m_coordim; i++)
717  {
718  faceAaxis[i] = (*m_faces[f]->GetVertex(1))[i] - (*m_faces[f]->GetVertex(0))[i];
719  faceBaxis[i] = (*m_faces[f]->GetVertex(3))[i] - (*m_faces[f]->GetVertex(0))[i];
720 
721  elementAaxis_length += pow(elementAaxis[i],2);
722  elementBaxis_length += pow(elementBaxis[i],2);
723  faceAaxis_length += pow(faceAaxis[i],2);
724  faceBaxis_length += pow(faceBaxis[i],2);
725  }
726 
727  elementAaxis_length = sqrt(elementAaxis_length);
728  elementBaxis_length = sqrt(elementBaxis_length);
729  faceAaxis_length = sqrt(faceAaxis_length);
730  faceBaxis_length = sqrt(faceBaxis_length);
731 
732  // Calculate the inner product of both the A-axis
733  // (i.e. Elemental A axis and face A axis)
734  for(i = 0 ; i < m_coordim; i++)
735  {
736  dotproduct1 += elementAaxis[i]*faceAaxis[i];
737  }
738 
739  orientation = 0;
740  // if the innerproduct is equal to the (absolute value of the ) products of the lengths
741  // of both vectors, then, the coordinate systems will NOT be transposed
742  if( fabs(elementAaxis_length*faceAaxis_length - fabs(dotproduct1)) < NekConstants::kNekZeroTol )
743  {
744  // if the inner product is negative, both A-axis point
745  // in reverse direction
746  if(dotproduct1 < 0.0)
747  {
748  orientation += 2;
749  }
750 
751  // calculate the inner product of both B-axis
752  for(i = 0 ; i < m_coordim; i++)
753  {
754  dotproduct2 += elementBaxis[i]*faceBaxis[i];
755  }
756 
757  // check that both these axis are indeed parallel
758  ASSERTL1(fabs(elementBaxis_length*faceBaxis_length - fabs(dotproduct2)) <
760  "These vectors should be parallel");
761 
762  // if the inner product is negative, both B-axis point
763  // in reverse direction
764  if( dotproduct2 < 0.0 )
765  {
766  orientation++;
767  }
768  }
769  // The coordinate systems are transposed
770  else
771  {
772  orientation = 4;
773 
774  // Calculate the inner product between the elemental A-axis
775  // and the B-axis of the face (which are now the corresponding axis)
776  dotproduct1 = 0.0;
777  for(i = 0 ; i < m_coordim; i++)
778  {
779  dotproduct1 += elementAaxis[i]*faceBaxis[i];
780  }
781 
782  // check that both these axis are indeed parallel
783  ASSERTL1(fabs(elementAaxis_length*faceBaxis_length - fabs(dotproduct1)) <
785  "These vectors should be parallel");
786 
787  // if the result is negative, both axis point in reverse
788  // directions
789  if(dotproduct1 < 0.0)
790  {
791  orientation += 2;
792  }
793 
794  // Do the same for the other two corresponding axis
795  dotproduct2 = 0.0;
796  for(i = 0 ; i < m_coordim; i++)
797  {
798  dotproduct2 += elementBaxis[i]*faceAaxis[i];
799  }
800 
801  // check that both these axis are indeed parallel
802  ASSERTL1(fabs(elementBaxis_length*faceAaxis_length - fabs(dotproduct2)) <
804  "These vectors should be parallel");
805 
806  if( dotproduct2 < 0.0 )
807  {
808  orientation++;
809  }
810  }
811 
812  orientation = orientation + 5;
813  // Fill the m_forient array
814  m_forient[f] = (StdRegions::Orientation) orientation;
815  }
816  }
817 
819  {
820 
821  // This 2D array holds the local id's of all the vertices
822  // for every edge. For every edge, they are ordered to what we
823  // define as being Forwards
824  const unsigned int edgeVerts[kNedges][2] =
825  { {0,1} ,
826  {1,2} ,
827  {2,3} ,
828  {3,0} ,
829  {0,4} ,
830  {1,5} ,
831  {2,6} ,
832  {3,7} ,
833  {4,5} ,
834  {5,6} ,
835  {6,7} ,
836  {7,4} };
837 
838  int i;
839  for(i = 0; i < kNedges; i++)
840  {
841  if( m_edges[i]->GetVid(0) == m_verts[ edgeVerts[i][0] ]->GetVid() )
842  {
844  }
845  else if( m_edges[i]->GetVid(0) == m_verts[ edgeVerts[i][1] ]->GetVid() )
846  {
848  }
849  else
850  {
851  ASSERTL0(false,"Could not find matching vertex for the edge");
852  }
853  }
854  }
855 
857  CurveMap &curvedEdges,
858  CurveMap &curvedFaces)
859  {
860  Geometry::v_Reset(curvedEdges, curvedFaces);
861 
862  for (int i = 0; i < 6; ++i)
863  {
864  m_faces[i]->Reset(curvedEdges, curvedFaces);
865  }
866 
867  SetUpXmap();
868  SetUpCoeffs(m_xmap->GetNcoeffs());
869  }
870 
871  /**
872  * @brief Set up the #m_xmap object by determining the order of each
873  * direction from derived faces.
874  */
876  {
877  /// Determine necessary order for standard region. This can almost
878  /// certainly be simplified but works for now!
879  vector<int> tmp1, tmp2;
880 
881  if (m_forient[0] < 9)
882  {
883  tmp1.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs (0));
884  tmp1.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs (2));
885  tmp2.push_back(m_faces[0]->GetXmap()->GetEdgeNumPoints(0));
886  tmp2.push_back(m_faces[0]->GetXmap()->GetEdgeNumPoints(2));
887  }
888  else
889  {
890  tmp1.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs (1));
891  tmp1.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs (3));
892  tmp2.push_back(m_faces[0]->GetXmap()->GetEdgeNumPoints(1));
893  tmp2.push_back(m_faces[0]->GetXmap()->GetEdgeNumPoints(3));
894  }
895 
896  if (m_forient[5] < 9)
897  {
898  tmp1.push_back(m_faces[5]->GetXmap()->GetEdgeNcoeffs (0));
899  tmp1.push_back(m_faces[5]->GetXmap()->GetEdgeNcoeffs (2));
900  tmp2.push_back(m_faces[5]->GetXmap()->GetEdgeNumPoints(0));
901  tmp2.push_back(m_faces[5]->GetXmap()->GetEdgeNumPoints(2));
902  }
903  else
904  {
905  tmp1.push_back(m_faces[5]->GetXmap()->GetEdgeNcoeffs (1));
906  tmp1.push_back(m_faces[5]->GetXmap()->GetEdgeNcoeffs (3));
907  tmp2.push_back(m_faces[5]->GetXmap()->GetEdgeNumPoints(1));
908  tmp2.push_back(m_faces[5]->GetXmap()->GetEdgeNumPoints(3));
909  }
910 
911  int order0 = *max_element(tmp1.begin(), tmp1.end());
912 
913  tmp1.clear();
914  tmp2.clear();
915 
916  if (m_forient[0] < 9)
917  {
918  tmp1.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs (1));
919  tmp1.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs (3));
920  tmp2.push_back(m_faces[0]->GetXmap()->GetEdgeNumPoints(1));
921  tmp2.push_back(m_faces[0]->GetXmap()->GetEdgeNumPoints(3));
922  }
923  else
924  {
925  tmp1.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs (0));
926  tmp1.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs (2));
927  tmp2.push_back(m_faces[0]->GetXmap()->GetEdgeNumPoints(0));
928  tmp2.push_back(m_faces[0]->GetXmap()->GetEdgeNumPoints(2));
929  }
930 
931  if (m_forient[5] < 9)
932  {
933  tmp1.push_back(m_faces[5]->GetXmap()->GetEdgeNcoeffs (1));
934  tmp1.push_back(m_faces[5]->GetXmap()->GetEdgeNcoeffs (3));
935  tmp2.push_back(m_faces[5]->GetXmap()->GetEdgeNumPoints(1));
936  tmp2.push_back(m_faces[5]->GetXmap()->GetEdgeNumPoints(3));
937  }
938  else
939  {
940  tmp1.push_back(m_faces[5]->GetXmap()->GetEdgeNcoeffs (0));
941  tmp1.push_back(m_faces[5]->GetXmap()->GetEdgeNcoeffs (2));
942  tmp2.push_back(m_faces[5]->GetXmap()->GetEdgeNumPoints(0));
943  tmp2.push_back(m_faces[5]->GetXmap()->GetEdgeNumPoints(2));
944  }
945 
946  int order1 = *max_element(tmp1.begin(), tmp1.end());
947 
948  tmp1.clear();
949  tmp2.clear();
950 
951  if (m_forient[1] < 9)
952  {
953  tmp1.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs (1));
954  tmp1.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs (3));
955  tmp2.push_back(m_faces[1]->GetXmap()->GetEdgeNumPoints(1));
956  tmp2.push_back(m_faces[1]->GetXmap()->GetEdgeNumPoints(3));
957  }
958  else
959  {
960  tmp1.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs (0));
961  tmp1.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs (2));
962  tmp2.push_back(m_faces[1]->GetXmap()->GetEdgeNumPoints(0));
963  tmp2.push_back(m_faces[1]->GetXmap()->GetEdgeNumPoints(2));
964  }
965 
966  if (m_forient[3] < 9)
967  {
968  tmp1.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs (1));
969  tmp1.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs (3));
970  tmp2.push_back(m_faces[3]->GetXmap()->GetEdgeNumPoints(1));
971  tmp2.push_back(m_faces[3]->GetXmap()->GetEdgeNumPoints(3));
972  }
973  else
974  {
975  tmp1.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs (0));
976  tmp1.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs (2));
977  tmp2.push_back(m_faces[3]->GetXmap()->GetEdgeNumPoints(0));
978  tmp2.push_back(m_faces[3]->GetXmap()->GetEdgeNumPoints(2));
979  }
980 
981  int order2 = *max_element(tmp1.begin(), tmp1.end());
982 
983  const LibUtilities::BasisKey A(
987  const LibUtilities::BasisKey B(
991  const LibUtilities::BasisKey C(
995 
997  A, B, C);
998  }
999  }; //end of namespace
1000 }; //end of namespace
static const unsigned int VertexEdgeConnectivity[8][3]
Definition: HexGeom.h:105
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
void SetUpXmap()
Set up the m_xmap object by determining the order of each direction from derived faces.
Definition: HexGeom.cpp:875
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry3D.h:92
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual void v_GenGeomFactors()
Definition: HexGeom.cpp:88
GeomFactorsSharedPtr m_geomFactors
Definition: Geometry.h:170
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:765
static const int kNtfaces
Definition: HexGeom.h:63
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:857
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:833
int GetEid(int i) const
Return the ID of edge i in this element.
Definition: Geometry3D.cpp:74
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Principle Modified Functions .
Definition: BasisType.h:49
virtual int v_GetNumFaces() const
Definition: HexGeom.cpp:350
STL namespace.
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
Definition: Geometry.cpp:307
static const int kNfaces
Definition: HexGeom.h:64
virtual int v_GetNumEdges() const
Definition: HexGeom.cpp:345
int GetFid(int i) const
Definition: Geometry.h:329
StdRegions::StdExpansionSharedPtr GetXmap() const
Definition: Geometry.h:383
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:93
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
static const NekDouble kNekZeroTol
virtual bool v_ContainsPoint(const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0)
Determines if a point specified in global coordinates is located within this tetrahedral geometry...
Definition: HexGeom.cpp:251
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
virtual NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Definition: HexGeom.cpp:148
virtual int v_GetVertexEdgeMap(const int i, const int j) const
Definition: HexGeom.cpp:355
static std::string npts
Definition: InputFld.cpp:43
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
Definition: HexGeom.cpp:856
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
PointGeomSharedPtr GetVertex(int i) const
Definition: Geometry.h:348
virtual void v_FillGeom()
Put all quadrature information into face/edge structure and backward transform.
Definition: Geometry3D.cpp:247
virtual int v_GetEdgeFaceMap(const int i, const int j) const
Definition: HexGeom.cpp:365
const Geometry1DSharedPtr GetEdge(int i) const
3D geometry information
Definition: Geometry3D.h:70
virtual int v_GetVertexFaceMap(const int i, const int j) const
Definition: HexGeom.cpp:360
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Definition: Geometry.h:180
static const unsigned int EdgeFaceConnectivity[12][2]
Definition: HexGeom.h:107
virtual int v_GetNumVerts() const
Definition: HexGeom.cpp:340
Geometry is straight-sided with constant geometric factors.
int GetVid(int i) const
Definition: Geometry.h:319
virtual int v_GetDir(const int faceidx, const int facedir) const
Definition: HexGeom.cpp:370
LibUtilities::ShapeType m_shapeType
Definition: Geometry.h:177
static const unsigned int VertexFaceConnectivity[8][3]
Definition: HexGeom.h:106
NekDouble dist(PointGeom &a)
Definition: PointGeom.cpp:186
boost::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:63
Geometric information has been generated.
void SetUpCoeffs(const int nCoeffs)
Initialise the m_coeffs array.
Definition: Geometry.h:484
static const int kNqfaces
Definition: HexGeom.h:62
GeomFactorsSharedPtr GetMetricInfo()
Definition: Geometry.h:299
GeomType
Indicates the type of element geometry.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
static const int kNedges
Definition: HexGeom.h:61
Geometry is curved or has non-constant factors.
int m_coordim
coordinate dimension
Definition: Geometry.h:169
Describes the specification for a Basis.
Definition: Basis.h:50
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169