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