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