Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TetGeom.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: TetGeom.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: Tetrahedral geometry information.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <SpatialDomains/TetGeom.h>
37 
39 #include <StdRegions/StdTetExp.h>
40 #include <SpatialDomains/SegGeom.h>
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46  namespace SpatialDomains
47  {
48  const unsigned int TetGeom::VertexEdgeConnectivity[4][3] = {
49  {0,2,3},{0,1,4},{1,2,5},{3,4,5}};
50  const unsigned int TetGeom::VertexFaceConnectivity[4][3] = {
51  {0,1,3},{0,1,2},{0,2,3},{1,2,3}};
52  const unsigned int TetGeom::EdgeFaceConnectivity [6][2] = {
53  {0,1},{0,2},{0,3},{1,3},{1,2},{2,3}};
54 
55  TetGeom::TetGeom()
56  {
57  m_shapeType = LibUtilities::eTetrahedron;
58  }
59 
60  TetGeom::TetGeom(const TriGeomSharedPtr faces[]) :
61  Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
62  {
64 
65  /// Copy the face shared pointers
66  m_faces.insert(m_faces.begin(), faces, faces+TetGeom::kNfaces);
67 
68  /// Set up orientation vectors with correct amount of elements.
69  m_eorient.resize(kNedges);
70  m_forient.resize(kNfaces);
71 
76  SetUpXmap();
77  SetUpCoeffs(m_xmap->GetNcoeffs());
78  }
79 
81  {
82 
83  }
84 
85  /**
86  * @brief Determines if a point specified in global coordinates is
87  * located within this tetrahedral geometry.
88  */
90  const Array<OneD, const NekDouble> &gloCoord, NekDouble tol)
91  {
92  Array<OneD,NekDouble> locCoord(GetCoordim(),0.0);
93  return v_ContainsPoint(gloCoord,locCoord,tol);
94  }
95 
97  Array<OneD, NekDouble> &locCoord,
98  NekDouble tol)
99  {
100  NekDouble resid;
101  return v_ContainsPoint(gloCoord,locCoord,tol,resid);
102  }
103 
104  /**
105  * @brief Determines if a point specified in global coordinates is
106  * located within this tetrahedral geometry and return local caretsian coordinates
107  */
109  Array<OneD, NekDouble> &locCoord,
110  NekDouble tol,
111  NekDouble &resid)
112  {
113  // Validation checks
114  ASSERTL1(gloCoord.num_elements() == 3,
115  "Three dimensional geometry expects three coordinates.");
116 
117  // find min, max point and check if within twice this
118  // distance other false this is advisable since
119  // GetLocCoord is expensive for non regular elements.
120  if(GetMetricInfo()->GetGtype() != eRegular)
121  {
122  int i;
123  Array<OneD, NekDouble> mincoord(3), maxcoord(3);
124  NekDouble diff = 0.0;
125 
126  v_FillGeom();
127 
128  const int npts = m_xmap->GetTotPoints();
129  Array<OneD, NekDouble> pts(npts);
130 
131  for(i = 0; i < 3; ++i)
132  {
133  m_xmap->BwdTrans(m_coeffs[i], pts);
134 
135  mincoord[i] = Vmath::Vmin(pts.num_elements(),pts,1);
136  maxcoord[i] = Vmath::Vmax(pts.num_elements(),pts,1);
137 
138  diff = max(maxcoord[i] - mincoord[i],diff);
139  }
140 
141  for(i = 0; i < 3; ++i)
142  {
143  if((gloCoord[i] < mincoord[i] - 0.2*diff)||
144  (gloCoord[i] > maxcoord[i] + 0.2*diff))
145  {
146  return false;
147  }
148  }
149  }
150 
151  // Convert to the local (eta) coordinates.
152  resid = v_GetLocCoords(gloCoord, locCoord);
153 
154  // Check local coordinate is within cartesian bounds.
155  if (locCoord[0] >= -(1+tol) && locCoord[1] >= -(1+tol) &&
156  locCoord[2] >= -(1+tol) &&
157  locCoord[0] + locCoord[1] + locCoord[2] <= -1+tol)
158  {
159  return true;
160  }
161 
162  // If out of range clamp locCoord to be within [-1,1]^3
163  // since any larger value will be very oscillatory if
164  // called by 'returnNearestElmt' option in
165  // ExpList::GetExpIndex
166  for(int i = 0; i < 3; ++i)
167  {
168  if(locCoord[i] <-(1+tol))
169  {
170  locCoord[i] = -(1+tol);
171  }
172 
173  if(locCoord[i] > (1+tol))
174  {
175  locCoord[i] = 1+tol;
176  }
177  }
178 
179  return false;
180  }
181 
182 
183  /// Get Local cartesian points
185  const Array<OneD, const NekDouble>& coords,
186  Array<OneD, NekDouble>& Lcoords)
187  {
188  NekDouble ptdist = 1e6;
189 
190  // calculate local coordinates (eta) for coord
191  if(GetMetricInfo()->GetGtype() == eRegular)
192  {
193  // Point inside tetrahedron
194  PointGeom r(m_coordim, 0, coords[0], coords[1], coords[2]);
195 
196  // Edges
197  PointGeom er0, e10, e20, e30;
198  er0.Sub(r,*m_verts[0]);
199  e10.Sub(*m_verts[1],*m_verts[0]);
200  e20.Sub(*m_verts[2],*m_verts[0]);
201  e30.Sub(*m_verts[3],*m_verts[0]);
202 
203 
204  // Cross products (Normal times area)
205  PointGeom cp1020, cp2030, cp3010;
206  cp1020.Mult(e10,e20);
207  cp2030.Mult(e20,e30);
208  cp3010.Mult(e30,e10);
209 
210 
211  // Barycentric coordinates (relative volume)
212  NekDouble V = e30.dot(cp1020); //Tet Volume={(e30)dot(e10)x(e20)}/6
213  NekDouble beta = er0.dot(cp2030)/V; //volume1={(er0)dot(e20)x(e30)}/6
214  NekDouble gamma = er0.dot(cp3010)/V; //volume1={(er0)dot(e30)x(e10)}/6
215  NekDouble delta = er0.dot(cp1020)/V; //volume1={(er0)dot(e10)x(e20)}/6
216 
217  // Make tet bigger
218  Lcoords[0] = 2.0*beta - 1.0;
219  Lcoords[1] = 2.0*gamma - 1.0;
220  Lcoords[2] = 2.0*delta - 1.0;
221 
222  // Set ptdist to distance to nearest vertex
223  for(int i = 0; i < 4; ++i)
224  {
225  ptdist = min(ptdist,r.dist(*m_verts[i]));
226  }
227  }
228  else
229  {
230  v_FillGeom();
231 
232  // Determine nearest point of coords to values in m_xmap
233  int npts = m_xmap->GetTotPoints();
234  Array<OneD, NekDouble> ptsx(npts), ptsy(npts), ptsz(npts);
235  Array<OneD, NekDouble> tmp1(npts), tmp2(npts);
236 
237  m_xmap->BwdTrans(m_coeffs[0], ptsx);
238  m_xmap->BwdTrans(m_coeffs[1], ptsy);
239  m_xmap->BwdTrans(m_coeffs[2], ptsz);
240 
241  const Array<OneD, const NekDouble> za = m_xmap->GetPoints(0);
242  const Array<OneD, const NekDouble> zb = m_xmap->GetPoints(1);
243  const Array<OneD, const NekDouble> zc = m_xmap->GetPoints(2);
244 
245  //guess the first local coords based on nearest point
246  Vmath::Sadd(npts, -coords[0], ptsx,1,tmp1,1);
247  Vmath::Vmul (npts, tmp1,1,tmp1,1,tmp1,1);
248  Vmath::Sadd(npts, -coords[1], ptsy,1,tmp2,1);
249  Vmath::Vvtvp(npts, tmp2,1,tmp2,1,tmp1,1,tmp1,1);
250  Vmath::Sadd(npts, -coords[2], ptsz,1,tmp2,1);
251  Vmath::Vvtvp(npts, tmp2,1,tmp2,1,tmp1,1,tmp1,1);
252 
253  int min_i = Vmath::Imin(npts,tmp1,1);
254 
255  // distance from coordinate to nearest point for return value.
256  ptdist = sqrt(tmp1[min_i]);
257 
258  // Get collapsed coordinate
259  int qa = za.num_elements(), qb = zb.num_elements();
260  Lcoords[2] = zc[min_i/(qa*qb)];
261  min_i = min_i%(qa*qb);
262  Lcoords[1] = zb[min_i/qa];
263  Lcoords[0] = za[min_i%qa];
264 
265  // recover cartesian coordinate from collapsed coordinate.
266  Lcoords[1] = (1.0+Lcoords[0])*(1.0-Lcoords[2])/2 -1.0;
267  Lcoords[0] = (1.0+Lcoords[0])*(-Lcoords[1]-Lcoords[2])/2 -1.0;
268 
269  // Perform newton iteration to find local coordinates
270  NekDouble resid = 0.0;
271  NewtonIterationForLocCoord(coords, ptsx, ptsy, ptsz, Lcoords,resid);
272  }
273  return ptdist;
274  }
275 
277  {
278  return 4;
279  }
280 
282  {
283  return 6;
284  }
285 
287  {
288  return 4;
289  }
290 
291  int TetGeom::v_GetDir(const int faceidx, const int facedir) const
292  {
293  if (faceidx == 0)
294  {
295  return facedir;
296  }
297  else if (faceidx == 1)
298  {
299  return 2 * facedir;
300  }
301  else
302  {
303  return 1 + facedir;
304  }
305  }
306 
307  int TetGeom::v_GetVertexEdgeMap(const int i, const int j) const
308  {
309  return VertexEdgeConnectivity[i][j];
310  }
311 
312  int TetGeom::v_GetVertexFaceMap(const int i, const int j) const
313  {
314  return VertexFaceConnectivity[i][j];
315  }
316 
317  int TetGeom::v_GetEdgeFaceMap(const int i, const int j) const
318  {
319  return EdgeFaceConnectivity[i][j];
320  }
321 
323  {
324 
325  // find edge 0
326  int i,j;
327  unsigned int check;
328 
329  SegGeomSharedPtr edge;
330 
331  // First set up the 3 bottom edges
332 
333  if(m_faces[0]->GetEid(0) != m_faces[1]->GetEid(0))
334  {
335  std::ostringstream errstrm;
336  errstrm << "Local edge 0 (eid=" << m_faces[0]->GetEid(0);
337  errstrm << ") on face " << m_faces[0]->GetFid();
338  errstrm << " must be the same as local edge 0 (eid="<<m_faces[1]->GetEid(0);
339  errstrm << ") on face " << m_faces[1]->GetFid();
340  ASSERTL0(false, errstrm.str());
341  }
342 
343  int faceConnected;
344  for(faceConnected = 1; faceConnected < 4 ; faceConnected++)
345  {
346  check = 0;
347  for(i = 0; i < 3; i++)
348  {
349  if( (m_faces[0])->GetEid(i) == (m_faces[faceConnected])->GetEid(0) )
350  {
351  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[0])->GetEdge(i));
352  m_edges.push_back(edge);
353  check++;
354  }
355  }
356 
357  if( check < 1 )
358  {
359  std::ostringstream errstrm;
360  errstrm << "Face 0 does not share an edge with first edge of adjacent face. Faces ";
361  errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[faceConnected])->GetFid();
362  ASSERTL0(false, errstrm.str());
363  }
364  else if( check > 1)
365  {
366  std::ostringstream errstrm;
367  errstrm << "Connected faces share more than one edge. Faces ";
368  errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[faceConnected])->GetFid();
369  ASSERTL0(false, errstrm.str());
370  }
371  }
372 
373 
374  // Then, set up the 3 vertical edges
375  check = 0;
376  for(i = 0; i < 3; i++) //Set up the vertical edge :face(1) and face(3)
377  {
378  for(j = 0; j < 3; j++)
379  {
380  if( (m_faces[1])->GetEid(i) == (m_faces[3])->GetEid(j) )
381  {
382  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[1])->GetEdge(i));
383  m_edges.push_back(edge);
384  check++;
385  }
386  }
387  }
388  if( check < 1 )
389  {
390  std::ostringstream errstrm;
391  errstrm << "Connected faces do not share an edge. Faces ";
392  errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[3])->GetFid();
393  ASSERTL0(false, errstrm.str());
394  }
395  else if( check > 1)
396  {
397  std::ostringstream errstrm;
398  errstrm << "Connected faces share more than one edge. Faces ";
399  errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[3])->GetFid();
400  ASSERTL0(false, errstrm.str());
401  }
402  // Set up vertical edges: face(1) through face(3)
403  for(faceConnected = 1; faceConnected < 3 ; faceConnected++)
404  {
405  check = 0;
406  for(i = 0; i < 3; i++)
407  {
408  for(j = 0; j < 3; j++)
409  {
410  if( (m_faces[faceConnected])->GetEid(i) == (m_faces[faceConnected+1])->GetEid(j) )
411  {
412  edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[faceConnected])->GetEdge(i));
413  m_edges.push_back(edge);
414  check++;
415  }
416  }
417  }
418 
419  if( check < 1 )
420  {
421  std::ostringstream errstrm;
422  errstrm << "Connected faces do not share an edge. Faces ";
423  errstrm << (m_faces[faceConnected])->GetFid() << ", " << (m_faces[faceConnected+1])->GetFid();
424  ASSERTL0(false, errstrm.str());
425  }
426  else if( check > 1)
427  {
428  std::ostringstream errstrm;
429  errstrm << "Connected faces share more than one edge. Faces ";
430  errstrm << (m_faces[faceConnected])->GetFid() << ", " << (m_faces[faceConnected+1])->GetFid();
431  ASSERTL0(false, errstrm.str());
432  }
433  }
434  };
435 
437 
438  // Set up the first 2 vertices (i.e. vertex 0,1)
439  if( ( m_edges[0]->GetVid(0) == m_edges[1]->GetVid(0) ) ||
440  ( m_edges[0]->GetVid(0) == m_edges[1]->GetVid(1) ) )
441  {
442  m_verts.push_back(m_edges[0]->GetVertex(1));
443  m_verts.push_back(m_edges[0]->GetVertex(0));
444  }
445  else if( ( m_edges[0]->GetVid(1) == m_edges[1]->GetVid(0) ) ||
446  ( m_edges[0]->GetVid(1) == m_edges[1]->GetVid(1) ) )
447  {
448  m_verts.push_back(m_edges[0]->GetVertex(0));
449  m_verts.push_back(m_edges[0]->GetVertex(1));
450  }
451  else
452  {
453  std::ostringstream errstrm;
454  errstrm << "Connected edges do not share a vertex. Edges ";
455  errstrm << m_edges[0]->GetEid() << ", " << m_edges[1]->GetEid();
456  ASSERTL0(false, errstrm.str());
457  }
458 
459  // set up the other bottom vertices (i.e. vertex 2)
460  for(int i = 1; i < 2; i++)
461  {
462  if( m_edges[i]->GetVid(0) == m_verts[i]->GetVid() )
463  {
464  m_verts.push_back(m_edges[i]->GetVertex(1));
465  }
466  else if( m_edges[i]->GetVid(1) == m_verts[i]->GetVid() )
467  {
468  m_verts.push_back(m_edges[i]->GetVertex(0));
469  }
470  else
471  {
472  std::ostringstream errstrm;
473  errstrm << "Connected edges do not share a vertex. Edges ";
474  errstrm << m_edges[i]->GetEid() << ", " << m_edges[i-1]->GetEid();
475  ASSERTL0(false, errstrm.str());
476  }
477  }
478 
479  // set up top vertex
480  if (m_edges[3]->GetVid(0) == m_verts[0]->GetVid())
481  {
482  m_verts.push_back(m_edges[3]->GetVertex(1));
483  }
484  else {
485  m_verts.push_back(m_edges[3]->GetVertex(0));
486  }
487 
488  // Check the other edges match up.
489  int check = 0;
490  for (int i = 4; i < 6; ++i)
491  {
492  if( (m_edges[i]->GetVid(0) == m_verts[i-3]->GetVid()
493  && m_edges[i]->GetVid(1) == m_verts[3]->GetVid())
494  ||(m_edges[i]->GetVid(1) == m_verts[i-3]->GetVid()
495  && m_edges[i]->GetVid(0) == m_verts[3]->GetVid()))
496  {
497  check++;
498  }
499  }
500  if (check != 2) {
501  std::ostringstream errstrm;
502  errstrm << "Connected edges do not share a vertex. Edges ";
503  errstrm << m_edges[3]->GetEid() << ", " << m_edges[2]->GetEid();
504  ASSERTL0(false, errstrm.str());
505  }
506  };
507 
509 
510  // This 2D array holds the local id's of all the vertices
511  // for every edge. For every edge, they are ordered to what we
512  // define as being Forwards
513  const unsigned int edgeVerts[kNedges][2] =
514  { {0,1},
515  {1,2},
516  {0,2},
517  {0,3},
518  {1,3},
519  {2,3} };
520 
521  int i;
522  for(i = 0; i < kNedges; i++)
523  {
524  if( m_edges[i]->GetVid(0) == m_verts[ edgeVerts[i][0] ]->GetVid() )
525  {
527  }
528  else if( m_edges[i]->GetVid(0) == m_verts[ edgeVerts[i][1] ]->GetVid() )
529  {
531  }
532  else
533  {
534  ASSERTL0(false,"Could not find matching vertex for the edge");
535  }
536  }
537 
538  };
539 
541 
542  int f,i;
543 
544  // These arrays represent the vector of the A and B
545  // coordinate of the local elemental coordinate system
546  // where A corresponds with the coordinate direction xi_i
547  // with the lowest index i (for that particular face)
548  // Coordinate 'B' then corresponds to the other local
549  // coordinate (i.e. with the highest index)
550  Array<OneD,NekDouble> elementAaxis(m_coordim);
551  Array<OneD,NekDouble> elementBaxis(m_coordim);
552 
553  // These arrays correspond to the local coordinate
554  // system of the face itself (i.e. the Geometry2D)
555  // faceAaxis correspond to the xi_0 axis
556  // faceBaxis correspond to the xi_1 axis
557  Array<OneD,NekDouble> faceAaxis(m_coordim);
558  Array<OneD,NekDouble> faceBaxis(m_coordim);
559 
560  // This is the base vertex of the face (i.e. the Geometry2D)
561  // This corresponds to thevertex with local ID 0 of the
562  // Geometry2D
563  unsigned int baseVertex;
564 
565  // The lenght of the vectors above
566  NekDouble elementAaxis_length;
567  NekDouble elementBaxis_length;
568  NekDouble faceAaxis_length;
569  NekDouble faceBaxis_length;
570 
571  // This 2D array holds the local id's of all the vertices
572  // for every face. For every face, they are ordered in such
573  // a way that the implementation below allows a unified approach
574  // for all faces.
575  const unsigned int faceVerts[kNfaces][TriGeom::kNverts] =
576  { {0,1,2} ,
577  {0,1,3} ,
578  {1,2,3} ,
579  {0,2,3} };
580 
581  NekDouble dotproduct1 = 0.0;
582  NekDouble dotproduct2 = 0.0;
583 
584  unsigned int orientation;
585 
586  // Loop over all the faces to set up the orientation
587  for(f = 0; f < kNqfaces + kNtfaces; f++)
588  {
589  // initialisation
590  elementAaxis_length = 0.0;
591  elementBaxis_length = 0.0;
592  faceAaxis_length = 0.0;
593  faceBaxis_length = 0.0;
594 
595  dotproduct1 = 0.0;
596  dotproduct2 = 0.0;
597 
598  baseVertex = m_faces[f]->GetVid(0);
599 
600  // We are going to construct the vectors representing the A and B axis
601  // of every face. These vectors will be constructed as a vector-representation
602  // of the edges of the face. However, for both coordinate directions, we can
603  // represent the vectors by two different edges. That's why we need to make sure that
604  // we pick the edge to which the baseVertex of the Geometry2D-representation of the face
605  // belongs...
606 
607  // Compute the length of edges on a base-face
608  if( baseVertex == m_verts[ faceVerts[f][0] ]->GetVid() )
609  {
610  for(i = 0; i < m_coordim; i++)
611  {
612  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
613  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
614  }
615  }
616  else if( baseVertex == m_verts[ faceVerts[f][1] ]->GetVid() )
617  {
618  for(i = 0; i < m_coordim; i++)
619  {
620  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
621  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
622  }
623  }
624  else if( baseVertex == m_verts[ faceVerts[f][2] ]->GetVid() )
625  {
626  for(i = 0; i < m_coordim; i++)
627  {
628  elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][2] ])[i];
629  elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
630  }
631  }
632  else
633  {
634  ASSERTL0(false, "Could not find matching vertex for the face");
635  }
636 
637  // Now, construct the edge-vectors of the local coordinates of
638  // the Geometry2D-representation of the face
639  for(i = 0; i < m_coordim; i++)
640  {
641  faceAaxis[i] = (*m_faces[f]->GetVertex(1))[i] - (*m_faces[f]->GetVertex(0))[i];
642  faceBaxis[i] = (*m_faces[f]->GetVertex(2))[i] - (*m_faces[f]->GetVertex(0))[i];
643 
644  elementAaxis_length += pow(elementAaxis[i],2);
645  elementBaxis_length += pow(elementBaxis[i],2);
646  faceAaxis_length += pow(faceAaxis[i],2);
647  faceBaxis_length += pow(faceBaxis[i],2);
648  }
649 
650  elementAaxis_length = sqrt(elementAaxis_length);
651  elementBaxis_length = sqrt(elementBaxis_length);
652  faceAaxis_length = sqrt(faceAaxis_length);
653  faceBaxis_length = sqrt(faceBaxis_length);
654 
655  // Calculate the inner product of both the A-axis
656  // (i.e. Elemental A axis and face A axis)
657  for(i = 0 ; i < m_coordim; i++)
658  {
659  dotproduct1 += elementAaxis[i]*faceAaxis[i];
660  }
661 
662  orientation = 0;
663  // if the innerproduct is equal to the (absolute value of the ) products of the lengths
664  // of both vectors, then, the coordinate systems will NOT be transposed
665  if( fabs(elementAaxis_length*faceAaxis_length - fabs(dotproduct1)) < NekConstants::kNekZeroTol )
666  {
667  // if the inner product is negative, both A-axis point
668  // in reverse direction
669  if(dotproduct1 < 0.0)
670  {
671  orientation += 2;
672  }
673 
674  // calculate the inner product of both B-axis
675  for(i = 0 ; i < m_coordim; i++)
676  {
677  dotproduct2 += elementBaxis[i]*faceBaxis[i];
678  }
679 
680  // check that both these axis are indeed parallel
681  ASSERTL1(fabs(elementBaxis_length*faceBaxis_length - fabs(dotproduct2)) <
683  "These vectors should be parallel");
684 
685  // if the inner product is negative, both B-axis point
686  // in reverse direction
687  if( dotproduct2 < 0.0 )
688  {
689  orientation++;
690  }
691  }
692  // The coordinate systems are transposed
693  else
694  {
695  orientation = 4;
696 
697  // Calculate the inner product between the elemental A-axis
698  // and the B-axis of the face (which are now the corresponding axis)
699  dotproduct1 = 0.0;
700  for(i = 0 ; i < m_coordim; i++)
701  {
702  dotproduct1 += elementAaxis[i]*faceBaxis[i];
703  }
704 
705  // check that both these axis are indeed parallel
706  ASSERTL1(fabs(elementAaxis_length*faceBaxis_length - fabs(dotproduct1)) <
708  "These vectors should be parallel");
709 
710  // if the result is negative, both axis point in reverse
711  // directions
712  if(dotproduct1 < 0.0)
713  {
714  orientation += 2;
715  }
716 
717  // Do the same for the other two corresponding axis
718  dotproduct2 = 0.0;
719  for(i = 0 ; i < m_coordim; i++)
720  {
721  dotproduct2 += elementBaxis[i]*faceAaxis[i];
722  }
723 
724  // check that both these axis are indeed parallel
725  ASSERTL1(fabs(elementBaxis_length*faceAaxis_length - fabs(dotproduct2)) <
727  "These vectors should be parallel");
728 
729  if( dotproduct2 < 0.0 )
730  {
731  orientation++;
732  }
733  }
734 
735  orientation = orientation + 5;
736 
737  // Fill the m_forient array
738  m_forient[f] = (StdRegions::Orientation) orientation;
739  }
740  }
741 
743  CurveMap &curvedEdges,
744  CurveMap &curvedFaces)
745  {
746  Geometry::v_Reset(curvedEdges, curvedFaces);
747 
748  for (int i = 0; i < 4; ++i)
749  {
750  m_faces[i]->Reset(curvedEdges, curvedFaces);
751  }
752 
753  SetUpXmap();
754  SetUpCoeffs(m_xmap->GetNcoeffs());
755  }
756 
757  /**
758  * @brief Set up the #m_xmap object by determining the order of each
759  * direction from derived faces.
760  */
762  {
763  vector<int> tmp;
764  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
765  int order0 = *max_element(tmp.begin(), tmp.end());
766 
767  tmp.clear();
768  tmp.push_back(order0);
769  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
770  tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
771  int order1 = *max_element(tmp.begin(), tmp.end());
772 
773  tmp.clear();
774  tmp.push_back(order0);
775  tmp.push_back(order1);
776  tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(1));
777  tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(2));
778  tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(1));
779  int order2 = *max_element(tmp.begin(), tmp.end());
780 
785  const LibUtilities::BasisKey B(
789  const LibUtilities::BasisKey C(
793 
795  A, B, C);
796  }
797  }; //end of namespace
798 }; //end of namespace
virtual int v_GetEdgeFaceMap(const int i, const int j) const
Definition: TetGeom.cpp:317
StdRegions::StdExpansionSharedPtr m_xmap
Definition: Geometry.h:172
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
std::vector< StdRegions::Orientation > m_eorient
Definition: Geometry3D.h:92
Principle Modified Functions .
Definition: BasisType.h:51
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
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: TetGeom.cpp:89
static const int kNverts
Definition: TriGeom.h:99
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:779
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:871
NekDouble dot(PointGeom &a)
Definition: PointGeom.cpp:192
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:847
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:442
Principle Modified Functions .
Definition: BasisType.h:49
virtual int v_GetDir(const int faceidx, const int facedir) const
Definition: TetGeom.cpp:291
void Sub(PointGeom &a, PointGeom &b)
Definition: PointGeom.cpp:168
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
int GetFid(int i) const
Definition: Geometry.h:329
static const int kNtfaces
Definition: TetGeom.h:59
StdRegions::StdExpansionSharedPtr GetXmap() const
Definition: Geometry.h:383
static const int kNfaces
Definition: TetGeom.h:60
std::vector< StdRegions::Orientation > m_forient
Definition: Geometry3D.h:93
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
static const unsigned int EdgeFaceConnectivity[6][2]
Definition: TetGeom.h:103
virtual int v_GetNumFaces() const
Definition: TetGeom.cpp:286
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
static const int kNedges
Definition: TetGeom.h:57
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Principle Modified Functions .
Definition: BasisType.h:50
virtual int v_GetNumEdges() const
Definition: TetGeom.cpp:281
static const unsigned int VertexFaceConnectivity[4][3]
Definition: TetGeom.h:102
static std::string npts
Definition: InputFld.cpp:43
void SetUpXmap()
Set up the m_xmap object by determining the order of each direction from derived faces.
Definition: TetGeom.cpp:761
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
virtual int v_GetVertexEdgeMap(const int i, const int j) const
Definition: TetGeom.cpp:307
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:315
static const int kNqfaces
Definition: TetGeom.h:58
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 NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Get Local cartesian points.
Definition: TetGeom.cpp:184
void Mult(PointGeom &a, PointGeom &b)
Definition: PointGeom.cpp:177
const Geometry1DSharedPtr GetEdge(int i) const
virtual void v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
Reset this geometry object: unset the current state and remove allocated GeomFactors.
Definition: TetGeom.cpp:742
3D geometry information
Definition: Geometry3D.h:70
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Definition: Geometry.h:180
Geometry is straight-sided with constant geometric factors.
virtual int v_GetVertexFaceMap(const int i, const int j) const
Definition: TetGeom.cpp:312
int GetVid(int i) const
Definition: Geometry.h:319
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:60
virtual int v_GetNumVerts() const
Definition: TetGeom.cpp:276
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
LibUtilities::ShapeType m_shapeType
Definition: Geometry.h:177
NekDouble dist(PointGeom &a)
Definition: PointGeom.cpp:186
boost::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:63
void SetUpCoeffs(const int nCoeffs)
Initialise the m_coeffs array.
Definition: Geometry.h:484
static const unsigned int VertexEdgeConnectivity[4][3]
Definition: TetGeom.h:101
GeomFactorsSharedPtr GetMetricInfo()
Definition: Geometry.h:299
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
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:52
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:183