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