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