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