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