Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MeshElements.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MeshElements.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: Mesh manipulation objects.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <vector>
38 #include <loki/Singleton.h>
39 
40 #include <StdRegions/StdNodalTriExp.h>
41 #include <StdRegions/StdNodalTetExp.h>
42 #include <StdRegions/StdNodalPrismExp.h>
43 #include <LocalRegions/TriExp.h>
44 #include <LocalRegions/QuadExp.h>
45 #include <LocalRegions/TetExp.h>
46 #include <LocalRegions/PrismExp.h>
47 #include <SpatialDomains/SegGeom.h>
48 #include <SpatialDomains/TetGeom.h>
49 #include <SpatialDomains/PyrGeom.h>
50 #include <SpatialDomains/PrismGeom.h>
51 #include <SpatialDomains/HexGeom.h>
52 
53 #include "MeshElements.h"
54 using namespace std;
55 
56 namespace Nektar
57 {
58  namespace Utilities
59  {
60  ElementFactory& GetElementFactory()
61  {
62  typedef Loki::SingletonHolder<ElementFactory,
63  Loki::CreateUsingNew,
64  Loki::NoDestroy > Type;
65  return Type::Instance();
66  }
67 
68  Element::Element(ElmtConfig pConf,
69  unsigned int pNumNodes,
70  unsigned int pGotNodes)
71  : m_conf(pConf),
72  m_curveType(LibUtilities::ePolyEvenlySpaced),
73  m_geom()
74  {
75  if (pNumNodes != pGotNodes)
76  {
77  cerr << "Number of modes mismatch for type "
78  << pConf.m_e << "! Should be " << pNumNodes
79  << " but got " << pGotNodes << " nodes." << endl;
80  abort();
81  }
82  }
83 
84  /**
85  * @brief Return the number of elements of the expansion dimension.
86  */
87  unsigned int Mesh::GetNumElements()
88  {
89  return m_element[m_expDim].size();
90  }
91 
92  /**
93  * @brief Return the number of boundary elements (i.e. one below the
94  * expansion dimension).
95  */
96  unsigned int Mesh::GetNumBndryElements()
97  {
98  unsigned int i, nElmt = 0;
99 
100  for (i = 0; i < m_expDim; ++i)
101  nElmt += m_element[i].size();
102 
103  return nElmt;
104  }
105 
106  /**
107  * @brief Return the total number of entities in the mesh (i.e. all
108  * elements, regardless of dimension).
109  */
110  unsigned int Mesh::GetNumEntities()
111  {
112  unsigned int nEnt = 0;
113 
114  for (unsigned int d = 0; d <= m_expDim; ++d)
115  {
116  nEnt += m_element[d].size();
117  }
118 
119  return nEnt;
120  }
121 
122  /**
123  * @brief Test equality of two conditions - i.e. compare types, fields
124  * and values but _not_ composite ids.
125  */
126  bool operator==(ConditionSharedPtr const &c1, ConditionSharedPtr const &c2)
127  {
128  int i, n = c1->type.size();
129 
130  if (n != c2->type.size())
131  {
132  return false;
133  }
134 
135  for (i = 0; i < n; ++i)
136  {
137  if (c1->type[i] != c2->type[i])
138  {
139  return false;
140  }
141 
142  if (c1->field[i] != c2->field[i] || c1->value[i] != c2->value[i])
143  {
144  return false;
145  }
146  }
147 
148  return true;
149  }
150 
151  /**
152  * @brief Defines equality between two #NodeSharedPtr objects.
153  */
154  bool operator==(NodeSharedPtr const &p1, NodeSharedPtr const &p2)
155  {
156  return *p1 == *p2;
157  }
158 
159  /**
160  * @brief Defines ordering between two #NodeSharedPtr objects.
161  */
162  bool operator< (NodeSharedPtr const &p1, NodeSharedPtr const &p2)
163  {
164  return *p1 < *p2;
165  }
166 
167  std::ostream &operator<<(std::ostream &os, const NodeSharedPtr &n)
168  {
169  os << n->m_x << " " << n->m_y << " " << n->m_z;
170  return os;
171  }
172 
173  /**
174  * @brief Defines equality of two edges (equal if IDs of end nodes
175  * match in either ordering).
176  */
177  bool operator==(EdgeSharedPtr const &p1, EdgeSharedPtr const &p2)
178  {
179  return ( ((*(p1->m_n1) == *(p2->m_n1)) && (*(p1->m_n2) == *(p2->m_n2)))
180  || ((*(p1->m_n2) == *(p2->m_n1)) && (*(p1->m_n1) == *(p2->m_n2))));
181  }
182 
183  /**
184  * @brief Defines ordering between two edges (based on ID of edges).
185  */
186  bool operator< (EdgeSharedPtr const &p1, EdgeSharedPtr const &p2)
187  {
188  return p1->m_id < p2->m_id;
189  }
190 
191  /**
192  * @brief Defines equality of two faces (equal if IDs of vertices are
193  * the same.)
194  */
195  bool operator==(FaceSharedPtr const &p1, FaceSharedPtr const &p2)
196  {
197  std::vector<NodeSharedPtr>::iterator it1, it2;
198  for (it1 = p1->m_vertexList.begin(); it1 != p1->m_vertexList.end(); ++it1)
199  {
200  if (find(p2->m_vertexList.begin(), p2->m_vertexList.end(), *it1)
201  == p2->m_vertexList.end())
202  {
203  return false;
204  }
205  }
206  return true;
207  }
208 
209  /**
210  * @brief Defines ordering between two faces (depending on ID of
211  * faces).
212  */
213  bool operator< (FaceSharedPtr const &p1, FaceSharedPtr const &p2)
214  {
215  return p1->m_id < p2->m_id;
216  }
217 
218  /**
219  * @brief Replace a vertex in the element.
220  *
221  * When a vertex is replaced, the element edges and faces are also
222  * searched and the corresponding edge/face nodes are updated to
223  * maintain consistency.
224  *
225  * @param p Index of the vertex to replace.
226  * @param pNew New vertex.
227  */
228  void Element::SetVertex(unsigned int p, NodeSharedPtr pNew)
229  {
230  NodeSharedPtr vOld = m_vertex[p];
231  m_vertex[p] = pNew;
232  for (unsigned int i = 0; i < m_edge.size(); ++i)
233  {
234  if (m_edge[i]->m_n1 == vOld)
235  {
236  m_edge[i]->m_n1 = pNew;
237  }
238  else if (m_edge[i]->m_n2 == vOld)
239  {
240  m_edge[i]->m_n2 = pNew;
241  }
242  }
243  for (unsigned int i = 0; i < m_face.size(); ++i)
244  {
245  // Replace vertices in faces
246  for (unsigned int j = 0; j < m_face[i]->m_vertexList.size(); ++j)
247  {
248  if (m_face[i]->m_vertexList[j] == vOld)
249  {
250  m_face[i]->m_vertexList[j] = pNew;
251  }
252  }
253  for (unsigned int j = 0; j < m_face[i]->m_edgeList.size(); ++j)
254  {
255  if (m_face[i]->m_edgeList[j]->m_n1 == vOld)
256  {
257  m_face[i]->m_edgeList[j]->m_n1 = pNew;
258  }
259  else if (m_face[i]->m_edgeList[j]->m_n2 == vOld)
260  {
261  m_face[i]->m_edgeList[j]->m_n2 = pNew;
262  }
263  }
264  }
265  }
266 
267  /**
268  * @brief Replace an edge in the element.
269  *
270  * When an edge is replaced, the element faces are also searched and
271  * the corresponding face edges are updated to maintain consistency.
272  *
273  * @param p Index of the edge to replace.
274  * @param pNew New edge.
275  */
276  void Element::SetEdge(unsigned int p, EdgeSharedPtr pNew)
277  {
278  EdgeSharedPtr vOld = m_edge[p];
279  m_edge[p] = pNew;
280  for (unsigned int i = 0; i < m_face.size(); ++i)
281  {
282  for (unsigned int j = 0; j < m_face[i]->m_edgeList.size(); ++j)
283  {
284  if (m_face[i]->m_edgeList[j] == vOld)
285  {
286  m_face[i]->m_edgeList[j] = pNew;
287  }
288  }
289  }
290  }
291 
292  /**
293  * @brief Replace a face in the element.
294  *
295  * When a face is replaced, no other consistency checks are required.
296  *
297  * @param p Index of the face to replace.
298  * @param pNew New face.
299  */
300  void Element::SetFace(unsigned int p, FaceSharedPtr pNew)
301  {
302  m_face[p] = pNew;
303  }
304 
305  /**
306  * @brief Obtain the order of an element by looking at edges.
307  */
308  int Element::GetMaxOrder()
309  {
310  int i, ret = 1;
311 
312  for (i = 0; i < m_edge.size(); ++i)
313  {
314  int edgeOrder = m_edge[i]->GetNodeCount()-1;
315  if (edgeOrder > ret)
316  {
317  ret = edgeOrder;
318  }
319  }
320 
321  return ret;
322  }
323 
324  /**
325  * @brief Generate a Nektar++ string describing the composite.
326  *
327  * The list of composites may include individual element IDs or ranges
328  * of element IDs.
329  */
330  string Composite::GetXmlString(bool doSort)
331  {
332 
333 #if 0 // turn this option off since causes problem with InputNekpp.cpp
334  if (doSort)
335  {
336  element_id_less_than sortOperator;
337  sort(m_items.begin(), m_items.end(), sortOperator);
338  }
339 #endif
340 
341  stringstream st;
342  vector<ElementSharedPtr>::iterator it;
343  bool range = false;
344  int vId = m_items[0]->GetId();
345  int prevId = vId;
346 
347  st << " " << m_tag << "[" << vId;
348 
349  for (it = m_items.begin()+1; it != m_items.end(); ++it){
350  // store previous element ID and get current one
351  prevId = vId;
352  vId = (*it)->GetId();
353 
354  // continue an already started range
355  if (prevId > -1 && vId == prevId + 1)
356  {
357  range = true;
358  // if this is the last element, it's the end of a range, so write
359  if (*it == m_items.back())
360  {
361  st << "-" << vId;
362  }
363  continue;
364  }
365 
366  // terminate a range, if present
367  if (range)
368  {
369  st << "-" << prevId;
370  range = false;
371  }
372 
373  // write what will be either a single entry or start of new range
374  st << "," << vId;
375  }
376  // terminate
377  st << "] ";
378  return st.str();
379  }
380 
381  LibUtilities::ShapeType Point::m_type = GetElementFactory().
382  RegisterCreatorFunction(LibUtilities::ePoint, Point::create, "Point");
383 
384  /**
385  * @brief Create a point element.
386  */
387  Point::Point(ElmtConfig pConf,
388  vector<NodeSharedPtr> pNodeList,
389  vector<int> pTagList)
390  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
391  {
392  m_tag = "";
393  m_dim = 0;
394  m_taglist = pTagList;
395  m_vertex.push_back(pNodeList[0]);
396  }
397 
398  /**
399  * @brief Return the number of nodes defining a point (i.e. return 1).
400  */
401  unsigned int Point::GetNumNodes(ElmtConfig pConf)
402  {
403  return 1;
404  }
405 
406 
407  LibUtilities::ShapeType Line::m_type = GetElementFactory().
408  RegisterCreatorFunction(LibUtilities::eSegment, Line::create, "Line");
409 
410  /**
411  * @brief Create a line element.
412  */
413  Line::Line(ElmtConfig pConf,
414  vector<NodeSharedPtr> pNodeList,
415  vector<int> pTagList)
416  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
417  {
418  m_tag = "S";
419  m_dim = 1;
420  m_taglist = pTagList;
421  int n = m_conf.m_order-1;
422 
423  // Add vertices
424  for (int i = 0; i < 2; ++i) {
425  m_vertex.push_back(pNodeList[i]);
426  }
427  vector<NodeSharedPtr> edgeNodes;
428  if (m_conf.m_order > 1) {
429  for (int j = 0; j<n; ++j) {
430  edgeNodes.push_back(pNodeList[2+j]);
431  }
432  }
433  m_edge.push_back(boost::shared_ptr<Edge>(
434  new Edge(pNodeList[0], pNodeList[1], edgeNodes, m_conf.m_edgeCurveType)));
435  }
436 
437  SpatialDomains::GeometrySharedPtr Line::GetGeom(int coordDim)
438  {
439  // Create edge vertices.
440  SpatialDomains::PointGeomSharedPtr p[2];
441  SpatialDomains::SegGeomSharedPtr ret;
442 
443  p[0] = m_vertex[0]->GetGeom(coordDim);
444  p[1] = m_vertex[1]->GetGeom(coordDim);
445 
446  if (m_edge[0]->m_edgeNodes.size() > 0)
447  {
448  SpatialDomains::CurveSharedPtr c =
449  MemoryManager<SpatialDomains::Curve>::
450  AllocateSharedPtr(m_id, m_edge[0]->m_curveType);
451 
452  c->m_points.push_back(p[0]);
453  for (int i = 0; i < m_edge[0]->m_edgeNodes.size(); ++i)
454  {
455  c->m_points.push_back(m_edge[0]->m_edgeNodes[i]->GetGeom(coordDim));
456  }
457  c->m_points.push_back(p[1]);
458 
459  ret = MemoryManager<SpatialDomains::SegGeom>::
460  AllocateSharedPtr(m_id, 2, p, c);
461  }
462  else
463  {
464  ret = MemoryManager<SpatialDomains::SegGeom>::
465  AllocateSharedPtr(m_id, 2, p);
466  }
467 
468  return ret;
469  }
470 
471  /**
472  * @brief Return the number of nodes defining a line.
473  */
474  unsigned int Line::GetNumNodes(ElmtConfig pConf)
475  {
476  return pConf.m_order+1;
477  }
478 
479 
480  LibUtilities::ShapeType Triangle::m_type = GetElementFactory().
481  RegisterCreatorFunction(LibUtilities::eTriangle, Triangle::create, "Triangle");
482 
483  /**
484  * @brief Create a triangle element.
485  */
486  Triangle::Triangle(ElmtConfig pConf,
487  vector<NodeSharedPtr> pNodeList,
488  vector<int> pTagList)
489  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
490  {
491  m_tag = "T";
492  m_dim = 2;
493  m_taglist = pTagList;
494  m_curveType = LibUtilities::eNodalTriEvenlySpaced;
495  int n = m_conf.m_order-1;
496 
497  // Create a map to relate edge nodes to a pair of vertices
498  // defining an edge. This is based on the ordering produced by
499  // gmsh.
500  map<pair<int,int>, int> edgeNodeMap;
501  map<pair<int,int>, int>::iterator it;
502  edgeNodeMap[pair<int,int>(1,2)] = 4;
503  edgeNodeMap[pair<int,int>(2,3)] = 4 + n;
504  edgeNodeMap[pair<int,int>(3,1)] = 4 + 2*n;
505 
506  // Add vertices. This logic will determine (in 2D) whether the
507  // element is clockwise (sum > 0) or counter-clockwise (sum < 0).
508  NekDouble sum = 0.0;
509  for (int i = 0; i < 3; ++i) {
510  int o = (i+1) % 3;
511  m_vertex.push_back(pNodeList[i]);
512  sum += (pNodeList[o]->m_x - pNodeList[i]->m_x) *
513  (pNodeList[o]->m_y + pNodeList[i]->m_y);
514  }
515 
516  // Create edges (with corresponding set of edge points)
517  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
518  {
519  vector<NodeSharedPtr> edgeNodes;
520  if (m_conf.m_order > 1) {
521  for (int j = it->second; j < it->second + n; ++j) {
522  edgeNodes.push_back(pNodeList[j-1]);
523  }
524  }
525  m_edge.push_back(EdgeSharedPtr(new Edge(pNodeList[it->first.first-1],
526  pNodeList[it->first.second-1],
527  edgeNodes,
528  m_conf.m_edgeCurveType)));
529  }
530 
531  if (pConf.m_reorient)
532  {
533  if (sum > 0.0)
534  {
535  reverse(m_edge.begin(), m_edge.end());
536  }
537  }
538 
539  if (m_conf.m_faceNodes)
540  {
541  m_volumeNodes.insert(m_volumeNodes.begin(),
542  pNodeList.begin()+3*m_conf.m_order,
543  pNodeList.end());
544  }
545  }
546 
547  SpatialDomains::GeometrySharedPtr Triangle::GetGeom(int coordDim)
548  {
549  SpatialDomains::SegGeomSharedPtr edges[3];
550  SpatialDomains::PointGeomSharedPtr verts[3];
551  SpatialDomains::TriGeomSharedPtr ret;
552 
553  for (int i = 0; i < 3; ++i)
554  {
555  edges[i] = m_edge [i]->GetGeom(coordDim);
556  verts[i] = m_vertex[i]->GetGeom(coordDim);
557  }
558 
559  StdRegions::Orientation edgeorient[3] = {
560  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
561  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
562  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[2], *edges[0])
563  };
564 
565  ret = MemoryManager<SpatialDomains::TriGeom>::
566  AllocateSharedPtr(m_id, verts, edges, edgeorient);
567 
568  return ret;
569  }
570 
571  /**
572  * @brief Return the number of nodes defining a triangle.
573  */
574  unsigned int Triangle::GetNumNodes(ElmtConfig pConf)
575  {
576  int n = pConf.m_order;
577  if (!pConf.m_faceNodes)
578  return (n+1)+2*(n-1)+1;
579  else
580  return (n+1)*(n+2)/2;
581  }
582 
583  void Triangle::Complete(int order)
584  {
585  int i, j;
586 
587  // Create basis key for a nodal tetrahedron.
588  LibUtilities::BasisKey B0(
589  LibUtilities::eOrtho_A, order+1,
590  LibUtilities::PointsKey(
591  order+1,LibUtilities::eGaussLobattoLegendre));
592  LibUtilities::BasisKey B1(
593  LibUtilities::eOrtho_B, order+1,
594  LibUtilities::PointsKey(
595  order+1,LibUtilities::eGaussRadauMAlpha1Beta0));
596 
597  // Create a standard nodal triangle in order to get the
598  // Vandermonde matrix to perform interpolation to nodal points.
599  StdRegions::StdNodalTriExpSharedPtr nodalTri =
600  MemoryManager<StdRegions::StdNodalTriExp>::AllocateSharedPtr(
601  B0, B1, LibUtilities::eNodalTriEvenlySpaced);
602 
603  SpatialDomains::TriGeomSharedPtr geom =
604  boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
605  this->GetGeom(3));
606 
607  // Create basis key for a triangle.
608  LibUtilities::BasisKey C0(
609  LibUtilities::eOrtho_A, order+1,
610  LibUtilities::PointsKey(
611  order+1,LibUtilities::eGaussLobattoLegendre));
612  LibUtilities::BasisKey C1(
613  LibUtilities::eOrtho_B, order+1,
614  LibUtilities::PointsKey(
615  order+1,LibUtilities::eGaussRadauMAlpha1Beta0));
616 
617  // Create a triangle.
618  LocalRegions::TriExpSharedPtr tri =
619  MemoryManager<LocalRegions::TriExp>::AllocateSharedPtr(
620  C0, C1, geom);
621 
622  // Get coordinate array for tetrahedron.
623  int nqtot = tri->GetTotPoints();
624  Array<OneD, NekDouble> alloc(6*nqtot);
625  Array<OneD, NekDouble> xi(alloc);
626  Array<OneD, NekDouble> yi(alloc+ nqtot);
627  Array<OneD, NekDouble> zi(alloc+2*nqtot);
628  Array<OneD, NekDouble> xo(alloc+3*nqtot);
629  Array<OneD, NekDouble> yo(alloc+4*nqtot);
630  Array<OneD, NekDouble> zo(alloc+5*nqtot);
631  Array<OneD, NekDouble> tmp;
632 
633  tri->GetCoords(xi, yi, zi);
634 
635  for (i = 0; i < 3; ++i)
636  {
637  Array<OneD, NekDouble> coeffs(nodalTri->GetNcoeffs());
638  tri->FwdTrans(alloc+i*nqtot, coeffs);
639  // Apply Vandermonde matrix to project onto nodal space.
640  nodalTri->ModalToNodal(coeffs, tmp=alloc+(i+3)*nqtot);
641  }
642 
643  // Now extract points from the co-ordinate arrays into the
644  // edge/face/volume nodes. First, extract edge-interior nodes.
645  for (i = 0; i < 3; ++i)
646  {
647  int pos = 3 + i*(order-1);
648  m_edge[i]->m_edgeNodes.clear();
649  for (j = 0; j < order-1; ++j)
650  {
651  m_edge[i]->m_edgeNodes.push_back(
652  NodeSharedPtr(new Node(0, xo[pos+j], yo[pos+j], zo[pos+j])));
653  }
654  }
655 
656  // Extract face-interior nodes.
657  int pos = 3 + 3*(order-1);
658  for (i = pos; i < (order+1)*(order+2)/2; ++i)
659  {
660  m_volumeNodes.push_back(
661  NodeSharedPtr(new Node(0, xo[i], yo[i], zo[i])));
662  }
663 
664  m_conf.m_order = order;
665  m_conf.m_faceNodes = true;
666  m_conf.m_volumeNodes = true;
667  }
668 
669 
670  LibUtilities::ShapeType Quadrilateral::m_type = GetElementFactory().
671  RegisterCreatorFunction(LibUtilities::eQuadrilateral, Quadrilateral::create,
672  "Quadrilateral");
673 
674  /**
675  * @brief Create a quadrilateral element.
676  */
677  Quadrilateral::Quadrilateral(ElmtConfig pConf,
678  vector<NodeSharedPtr> pNodeList,
679  vector<int> pTagList)
680  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
681  {
682  m_tag = "Q";
683  m_dim = 2;
684  m_taglist = pTagList;
685  int n = m_conf.m_order-1;
686 
687  // Create a map to relate edge nodes to a pair of vertices
688  // defining an edge. This is based on the ordering produced by
689  // gmsh.
690  map<pair<int,int>, int> edgeNodeMap;
691  map<pair<int,int>, int>::iterator it;
692  edgeNodeMap[pair<int,int>(1,2)] = 5;
693  edgeNodeMap[pair<int,int>(2,3)] = 5 + n;
694  edgeNodeMap[pair<int,int>(3,4)] = 5 + 2*n;
695  edgeNodeMap[pair<int,int>(4,1)] = 5 + 3*n;
696 
697  // Add vertices. This logic will determine (in 2D) whether the
698  // element is clockwise (sum > 0) or counter-clockwise (sum < 0).
699  NekDouble sum = 0.0;
700  for (int i = 0; i < 4; ++i) {
701  int o = (i+1) % 4;
702  m_vertex.push_back(pNodeList[i]);
703  sum += (pNodeList[o]->m_x - pNodeList[i]->m_x) *
704  (pNodeList[o]->m_y + pNodeList[i]->m_y);
705  }
706 
707  // Create edges (with corresponding set of edge points)
708  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
709  {
710  vector<NodeSharedPtr> edgeNodes;
711  if (m_conf.m_order > 1) {
712  for (int j = it->second; j < it->second + n; ++j) {
713  edgeNodes.push_back(pNodeList[j-1]);
714  }
715  }
716  m_edge.push_back(EdgeSharedPtr(new Edge(pNodeList[it->first.first-1],
717  pNodeList[it->first.second-1],
718  edgeNodes,
719  m_conf.m_edgeCurveType)));
720  }
721 
722  if (pConf.m_reorient)
723  {
724  if (sum > 0.0)
725  {
726  reverse(m_edge.begin(), m_edge.end());
727  }
728  }
729 
730  if (m_conf.m_faceNodes)
731  {
732  m_volumeNodes.insert(m_volumeNodes.begin(),
733  pNodeList.begin()+4*m_conf.m_order,
734  pNodeList.end());
735  }
736  }
737 
738  void Quadrilateral::Complete(int order)
739  {
740  LibUtilities::BasisKey C0(
741  LibUtilities::eOrtho_A, order+1,
742  LibUtilities::PointsKey(
743  order+1,LibUtilities::eGaussLobattoLegendre));
744 
745  SpatialDomains::QuadGeomSharedPtr geom =
746  boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
747  this->GetGeom(3));
748 
749  // Create a quad.
750  LocalRegions::QuadExpSharedPtr quad =
751  MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(
752  C0, C0, geom);
753 
754  // Get coordinate array for quadrilateral.
755  int nqtot = quad->GetTotPoints();
756  Array<OneD, NekDouble> alloc(3*nqtot);
757  Array<OneD, NekDouble> x (alloc );
758  Array<OneD, NekDouble> y (alloc+1*nqtot);
759  Array<OneD, NekDouble> z (alloc+2*nqtot);
760 
761  quad->GetCoords(x, y, z);
762 
763  // Now extract points from the co-ordinate arrays into the edge
764  // and face nodes. First, extract edge-interior nodes.
765  int edgeMap[4][2] = {{0,1},{order,order+1},
766  {nqtot-1,-1},{order*(order+1),-order-1}};
767 
768  for (int i = 0; i < 4; ++i)
769  {
770  int pos = edgeMap[i][0] + edgeMap[i][1];
771  m_edge[i]->m_edgeNodes.clear();
772  for (int j = 1; j < order; ++j, pos += edgeMap[i][1])
773  {
774  m_edge[i]->m_edgeNodes.push_back(
775  NodeSharedPtr(new Node(0, x[pos], y[pos], z[pos])));
776  }
777  }
778 
779  // Extract face-interior nodes.
780  m_volumeNodes.clear();
781  for (int i = 1; i < order; ++i)
782  {
783  int pos = i*(order+1);
784  for (int j = 1; j < order; ++j)
785  {
786  m_volumeNodes.push_back(
787  NodeSharedPtr(new Node(0, x[pos+j], y[pos+j], z[pos+j])));
788  }
789  }
790 
791  m_conf.m_order = order;
792  m_conf.m_faceNodes = true;
793  m_conf.m_volumeNodes = true;
794  }
795 
796  SpatialDomains::GeometrySharedPtr Quadrilateral::GetGeom(int coordDim)
797  {
798  SpatialDomains::SegGeomSharedPtr edges[4];
799  SpatialDomains::PointGeomSharedPtr verts[4];
800  SpatialDomains::QuadGeomSharedPtr ret;
801 
802  for (int i = 0; i < 4; ++i)
803  {
804  edges[i] = m_edge [i]->GetGeom(coordDim);
805  verts[i] = m_vertex[i]->GetGeom(coordDim);
806  }
807 
808  StdRegions::Orientation edgeorient[4] = {
809  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
810  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
811  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
812  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
813  };
814 
815  ret = MemoryManager<SpatialDomains::QuadGeom>::
816  AllocateSharedPtr(m_id, verts, edges, edgeorient);
817 
818  return ret;
819  }
820 
821  /**
822  * @brief Return the number of nodes defining a quadrilateral.
823  */
824  unsigned int Quadrilateral::GetNumNodes(ElmtConfig pConf)
825  {
826  int n = pConf.m_order;
827  if (!pConf.m_faceNodes)
828  return 4*n;
829  else
830  return (n+1)*(n+1);
831  }
832 
833 
834  LibUtilities::ShapeType Tetrahedron::m_type = GetElementFactory().
835  RegisterCreatorFunction(LibUtilities::eTetrahedron, Tetrahedron::create, "Tetrahedron");
836 
837  /**
838  * @brief Create a tetrahedron element.
839  */
840  Tetrahedron::Tetrahedron(ElmtConfig pConf,
841  vector<NodeSharedPtr> pNodeList,
842  vector<int> pTagList)
843  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
844  {
845  m_tag = "A";
846  m_dim = 3;
847  m_taglist = pTagList;
848  int n = m_conf.m_order-1;
849 
850  // Create a map to relate edge nodes to a pair of vertices
851  // defining an edge.
852  map<pair<int,int>, int> edgeNodeMap;
853  map<pair<int,int>, int>::iterator it;
854  edgeNodeMap[pair<int,int>(1,2)] = 5;
855  edgeNodeMap[pair<int,int>(2,3)] = 5 + n;
856  edgeNodeMap[pair<int,int>(1,3)] = 5 + 2*n;
857  edgeNodeMap[pair<int,int>(1,4)] = 5 + 3*n;
858  edgeNodeMap[pair<int,int>(2,4)] = 5 + 4*n;
859  edgeNodeMap[pair<int,int>(3,4)] = 5 + 5*n;
860 
861  // Add vertices
862  for (int i = 0; i < 4; ++i)
863  {
864  m_vertex.push_back(pNodeList[i]);
865  }
866 
867  // Create edges (with corresponding set of edge points)
868  int eid = 0;
869  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
870  {
871  vector<NodeSharedPtr> edgeNodes;
872  if (m_conf.m_order > 1) {
873  for (int j = it->second; j < it->second + n; ++j) {
874  edgeNodes.push_back(pNodeList[j-1]);
875  }
876  }
877  m_edge.push_back(
878  EdgeSharedPtr(new Edge(pNodeList[it->first.first-1],
879  pNodeList[it->first.second-1],
880  edgeNodes,
881  m_conf.m_edgeCurveType)));
882  m_edge.back()->m_id = eid++;
883  }
884 
885  // Reorient the tet to ensure collapsed coordinates align between
886  // adjacent elements.
887  if (m_conf.m_reorient)
888  {
889  OrientTet();
890  }
891  else
892  {
893  // If we didn't need to orient the tet then set up the
894  // orientation map as the identity mapping.
895  for (int i = 0; i < 4; ++i)
896  {
897  m_orientationMap[i] = i;
898  }
899  }
900 
901  // Face-vertex IDs
902  int face_ids[4][3] = {
903  {0,1,2}, {0,1,3}, {1,2,3}, {0,2,3}
904  };
905 
906  // Face-edge IDs
907  int face_edges[4][3];
908 
909  // Create faces
910  for (int j = 0; j < 4; ++j)
911  {
912  vector<NodeSharedPtr> faceVertices;
913  vector<EdgeSharedPtr> faceEdges;
914  vector<NodeSharedPtr> faceNodes;
915 
916  // Extract the edges for this face. We need to do this because
917  // of the reorientation which might have been applied (see the
918  // additional note below).
919  for (int k = 0; k < 3; ++k)
920  {
921  faceVertices.push_back(m_vertex[face_ids[j][k]]);
922  NodeSharedPtr a = m_vertex[face_ids[j][k]];
923  NodeSharedPtr b = m_vertex[face_ids[j][(k+1)%3]];
924  for (unsigned int i = 0; i < m_edge.size(); ++i)
925  {
926  if ( ((*(m_edge[i]->m_n1)==*a) && (*(m_edge[i]->m_n2)==*b))
927  || ((*(m_edge[i]->m_n1)==*b) && (*(m_edge[i]->m_n2) == *a)) )
928  {
929  face_edges[j][k] = i;
930  faceEdges.push_back(m_edge[i]);
931  break;
932  }
933  }
934  }
935 
936  // When face curvature is supplied, it may have been the case
937  // that our tetrahedron was reoriented. In this case, faces have
938  // different vertex IDs and so we have to rotate the face
939  // curvature so that the two align appropriately.
940  if (m_conf.m_faceNodes)
941  {
942  const int nFaceNodes = n*(n-1)/2;
943 
944  // Get the vertex IDs of whatever face we are processing.
945  vector<int> faceIds(3);
946  faceIds[0] = faceVertices[0]->m_id;
947  faceIds[1] = faceVertices[1]->m_id;
948  faceIds[2] = faceVertices[2]->m_id;
949 
950  // Find out the original face number as we were given it in
951  // the constructor using the orientation map.
952  int origFace = -1;
953  for (int i = 0; i < 4; ++i)
954  {
955  if (m_orientationMap[i] == j)
956  {
957  origFace = i;
958  break;
959  }
960  }
961 
962  ASSERTL0(origFace >= 0, "Couldn't find face");
963 
964  // Now get the face nodes for the original face.
965  int N = 4 + 6*n + origFace * nFaceNodes;
966  for (int i = 0; i < nFaceNodes; ++i)
967  {
968  faceNodes.push_back(pNodeList[N+i]);
969  }
970 
971  // Find the original face vertex IDs.
972  vector<int> origFaceIds(3);
973  origFaceIds[0] = pNodeList[face_ids[origFace][0]]->m_id;
974  origFaceIds[1] = pNodeList[face_ids[origFace][1]]->m_id;
975  origFaceIds[2] = pNodeList[face_ids[origFace][2]]->m_id;
976 
977  // Construct a HOTriangle object which performs the
978  // orientation magically for us.
979  HOTriangle<NodeSharedPtr> hoTri(origFaceIds, faceNodes);
980  hoTri.Align(faceIds);
981 
982  // Copy the face nodes back again.
983  faceNodes = hoTri.surfVerts;
984  }
985  m_face.push_back(FaceSharedPtr(
986  new Face(faceVertices, faceNodes, faceEdges, m_conf.m_faceCurveType)));
987  }
988 
989  vector<EdgeSharedPtr> tmp(6);
990  tmp[0] = m_edge[face_edges[0][0]];
991  tmp[1] = m_edge[face_edges[0][1]];
992  tmp[2] = m_edge[face_edges[0][2]];
993  tmp[3] = m_edge[face_edges[1][2]];
994  tmp[4] = m_edge[face_edges[1][1]];
995  tmp[5] = m_edge[face_edges[2][1]];
996  m_edge = tmp;
997  }
998 
999  SpatialDomains::GeometrySharedPtr Tetrahedron::GetGeom(int coordDim)
1000  {
1001  SpatialDomains::TriGeomSharedPtr tfaces[4];
1002  SpatialDomains::TetGeomSharedPtr ret;
1003 
1004  for (int i = 0; i < 4; ++i)
1005  {
1006  tfaces[i] = boost::dynamic_pointer_cast
1007  <SpatialDomains::TriGeom>(m_face[i]->GetGeom(coordDim));
1008  }
1009 
1010  ret = MemoryManager<SpatialDomains::TetGeom>::
1011  AllocateSharedPtr(tfaces);
1012 
1013  return ret;
1014  }
1015 
1016  /**
1017  * @brief Return the number of nodes defining a tetrahedron.
1018  */
1019  unsigned int Tetrahedron::GetNumNodes(ElmtConfig pConf)
1020  {
1021  int n = pConf.m_order;
1022  if (pConf.m_volumeNodes && pConf.m_faceNodes)
1023  return (n+1)*(n+2)*(n+3)/6;
1024  else if (!pConf.m_volumeNodes && pConf.m_faceNodes)
1025  return 4*(n+1)*(n+2)/2-6*(n+1)+4;
1026  else
1027  return 6*(n+1)-8;
1028  }
1029 
1030  /**
1031  * @brief .
1032  */
1033  void Tetrahedron::Complete(int order)
1034  {
1035  int i, j;
1036 
1037  // Create basis key for a nodal tetrahedron.
1038  LibUtilities::BasisKey B0(
1039  LibUtilities::eOrtho_A, order+1,
1040  LibUtilities::PointsKey(
1041  order+1,LibUtilities::eGaussLobattoLegendre));
1042  LibUtilities::BasisKey B1(
1043  LibUtilities::eOrtho_B, order+1,
1044  LibUtilities::PointsKey(
1045  order+1,LibUtilities::eGaussRadauMAlpha1Beta0));
1046  LibUtilities::BasisKey B2(
1047  LibUtilities::eOrtho_C, order+1,
1048  LibUtilities::PointsKey(
1049  order+1,LibUtilities::eGaussRadauMAlpha2Beta0));
1050 
1051  // Create a standard nodal tetrahedron in order to get the
1052  // Vandermonde matrix to perform interpolation to nodal points.
1053  StdRegions::StdNodalTetExpSharedPtr nodalTet =
1054  MemoryManager<StdRegions::StdNodalTetExp>::AllocateSharedPtr(
1055  B0, B1, B2, LibUtilities::eNodalTetEvenlySpaced);
1056 
1057  Array<OneD, NekDouble> x, y, z;
1058 
1059  nodalTet->GetNodalPoints(x,y,z);
1060 
1061  SpatialDomains::TetGeomSharedPtr geom =
1062  boost::dynamic_pointer_cast<SpatialDomains::TetGeom>(
1063  this->GetGeom(3));
1064 
1065  // Create basis key for a tetrahedron.
1066  LibUtilities::BasisKey C0(
1067  LibUtilities::eOrtho_A, order+1,
1068  LibUtilities::PointsKey(
1069  order+1,LibUtilities::eGaussLobattoLegendre));
1070  LibUtilities::BasisKey C1(
1071  LibUtilities::eOrtho_B, order+1,
1072  LibUtilities::PointsKey(
1073  order+1,LibUtilities::eGaussRadauMAlpha1Beta0));
1074  LibUtilities::BasisKey C2(
1075  LibUtilities::eOrtho_C, order+1,
1076  LibUtilities::PointsKey(
1077  order+1,LibUtilities::eGaussRadauMAlpha2Beta0));
1078 
1079  // Create a tet.
1080  LocalRegions::TetExpSharedPtr tet =
1081  MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(
1082  C0, C1, C2, geom);
1083 
1084  // Get coordinate array for tetrahedron.
1085  int nqtot = tet->GetTotPoints();
1086  Array<OneD, NekDouble> alloc(6*nqtot);
1087  Array<OneD, NekDouble> xi(alloc);
1088  Array<OneD, NekDouble> yi(alloc+ nqtot);
1089  Array<OneD, NekDouble> zi(alloc+2*nqtot);
1090  Array<OneD, NekDouble> xo(alloc+3*nqtot);
1091  Array<OneD, NekDouble> yo(alloc+4*nqtot);
1092  Array<OneD, NekDouble> zo(alloc+5*nqtot);
1093  Array<OneD, NekDouble> tmp;
1094 
1095  tet->GetCoords(xi, yi, zi);
1096 
1097  for (i = 0; i < 3; ++i)
1098  {
1099  Array<OneD, NekDouble> coeffs(nodalTet->GetNcoeffs());
1100  tet->FwdTrans(alloc+i*nqtot, coeffs);
1101  // Apply Vandermonde matrix to project onto nodal space.
1102  nodalTet->ModalToNodal(coeffs, tmp=alloc+(i+3)*nqtot);
1103  }
1104 
1105  // Now extract points from the co-ordinate arrays into the
1106  // edge/face/volume nodes. First, extract edge-interior nodes.
1107  for (i = 0; i < 6; ++i)
1108  {
1109  int pos = 4 + i*(order-1);
1110  m_edge[i]->m_edgeNodes.clear();
1111  for (j = 0; j < order-1; ++j)
1112  {
1113  m_edge[i]->m_edgeNodes.push_back(
1114  NodeSharedPtr(new Node(0, xo[pos+j], yo[pos+j], zo[pos+j])));
1115  }
1116  }
1117 
1118  // Now extract face-interior nodes.
1119  for (i = 0; i < 4; ++i)
1120  {
1121  int pos = 4 + 6*(order-1) + i*(order-2)*(order-1)/2;
1122  m_face[i]->m_faceNodes.clear();
1123  for (j = 0; j < (order-2)*(order-1)/2; ++j)
1124  {
1125  m_face[i]->m_faceNodes.push_back(
1126  NodeSharedPtr(new Node(0, xo[pos+j], yo[pos+j], zo[pos+j])));
1127  }
1128  }
1129 
1130  // Finally extract volume nodes.
1131  int pos = 4 + 6*(order-1) + 4*(order-2)*(order-1)/2;
1132  for (i = pos; i < (order+1)*(order+2)*(order+3)/6; ++i)
1133  {
1134  m_volumeNodes.push_back(
1135  NodeSharedPtr(new Node(0, xo[i], yo[i], zo[i])));
1136  }
1137 
1138  m_conf.m_order = order;
1139  m_conf.m_faceNodes = true;
1140  m_conf.m_volumeNodes = true;
1141  }
1142 
1143  struct TetOrient
1144  {
1145  TetOrient(vector<int> nid, int fid) : nid(nid), fid(fid) {}
1146  vector<int> nid;
1147  int fid;
1148  };
1149 
1150  struct TetOrientHash : std::unary_function<struct TetOrient, std::size_t>
1151  {
1152  std::size_t operator()(struct TetOrient const& p) const
1153  {
1154  return boost::hash_range(p.nid.begin(), p.nid.end());
1155  }
1156  };
1157  typedef boost::unordered_set<struct TetOrient, TetOrientHash> TetOrientSet;
1158 
1159  bool operator==(const struct TetOrient &a, const struct TetOrient &b)
1160  {
1161  if (a.nid.size() != b.nid.size())
1162  {
1163  return false;
1164  }
1165 
1166  for (int i = 0; i < a.nid.size(); ++i)
1167  {
1168  if (a.nid[i] != b.nid[i])
1169  {
1170  return false;
1171  }
1172  }
1173 
1174  return true;
1175  }
1176 
1177  /**
1178  * @brief Orient tetrahedron to align degenerate vertices.
1179  *
1180  * Orientation of tetrahedral elements is required so that the
1181  * singular vertices of triangular faces (which occur as a part of the
1182  * collapsed co-ordinate system) align. The algorithm is based on that
1183  * used in T. Warburton's thesis and in the original Nektar source.
1184  *
1185  * First the vertices are ordered with the highest global vertex at
1186  * the top degenerate point, and the base degenerate point has second
1187  * lowest ID. These vertices are swapped if the element is incorrectly
1188  * oriented.
1189  */
1190  void Tetrahedron::OrientTet()
1191  {
1192  static int face_ids[4][3] = {
1193  {0,1,2},{0,1,3},{1,2,3},{0,2,3}};
1194 
1195  TetOrientSet faces;
1196 
1197  // Create a copy of the original vertex ordering. This is used to
1198  // construct a mapping, #orientationMap, which maps the original
1199  // face ordering to the new face ordering.
1200  for (int i = 0; i < 4; ++i)
1201  {
1202  vector<int> nodes(3);
1203 
1204  nodes[0] = m_vertex[face_ids[i][0]]->m_id;
1205  nodes[1] = m_vertex[face_ids[i][1]]->m_id;
1206  nodes[2] = m_vertex[face_ids[i][2]]->m_id;
1207 
1208  sort(nodes.begin(), nodes.end());
1209  struct TetOrient faceNodes(nodes, i);
1210  faces.insert(faceNodes);
1211  }
1212 
1213  // Store a copy of the original vertex ordering so we can create a
1214  // permutation map later.
1215  vector<NodeSharedPtr> origVert = m_vertex;
1216 
1217  // Order vertices with highest global vertex at top degenerate
1218  // point. Place second highest global vertex at base degenerate
1219  // point.
1220  sort(m_vertex.begin(), m_vertex.end());
1221 
1222  // Calculate a.(b x c) to determine tet volume; if negative,
1223  // reverse order of non-degenerate points to correctly orientate
1224  // the tet.
1225  double ax = m_vertex[1]->m_x-m_vertex[0]->m_x;
1226  double ay = m_vertex[1]->m_y-m_vertex[0]->m_y;
1227  double az = m_vertex[1]->m_z-m_vertex[0]->m_z;
1228  double bx = m_vertex[2]->m_x-m_vertex[0]->m_x;
1229  double by = m_vertex[2]->m_y-m_vertex[0]->m_y;
1230  double bz = m_vertex[2]->m_z-m_vertex[0]->m_z;
1231  double cx = m_vertex[3]->m_x-m_vertex[0]->m_x;
1232  double cy = m_vertex[3]->m_y-m_vertex[0]->m_y;
1233  double cz = m_vertex[3]->m_z-m_vertex[0]->m_z;
1234  double vol = cx*(ay*bz-az*by)+cy*(az*bx-ax*bz)+cz*(ax*by-ay*bx);
1235  vol /= 6.0;
1236 
1237  if (fabs(vol) <= 1e-10)
1238  {
1239  cerr << "Warning: degenerate tetrahedron, volume = " << vol << endl;
1240  }
1241 
1242  if (vol < 0)
1243  {
1244  swap(m_vertex[0], m_vertex[1]);
1245  }
1246 
1247  TetOrientSet::iterator it;
1248 
1249  // Search for the face in the original set of face nodes. Then use
1250  // this to construct the #orientationMap.
1251  for (int i = 0; i < 4; ++i)
1252  {
1253  vector<int> nodes(3);
1254 
1255  nodes[0] = m_vertex[face_ids[i][0]]->m_id;
1256  nodes[1] = m_vertex[face_ids[i][1]]->m_id;
1257  nodes[2] = m_vertex[face_ids[i][2]]->m_id;
1258  sort(nodes.begin(), nodes.end());
1259 
1260  struct TetOrient faceNodes(nodes, 0);
1261 
1262  it = faces.find(faceNodes);
1263  m_orientationMap[it->fid] = i;
1264 
1265  for (int j = 0; j < 4; ++j)
1266  {
1267  if (m_vertex[i]->m_id == origVert[j]->m_id)
1268  {
1269  m_origVertMap[j] = i;
1270  break;
1271  }
1272  }
1273  }
1274  }
1275 
1276  LibUtilities::ShapeType Pyramid::type = GetElementFactory().
1277  RegisterCreatorFunction(LibUtilities::ePyramid, Pyramid::create, "Pyramid");
1278 
1279  /**
1280  * @brief Create a pyramidic element.
1281  */
1282  Pyramid::Pyramid(ElmtConfig pConf,
1283  vector<NodeSharedPtr> pNodeList,
1284  vector<int> pTagList)
1285  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
1286  {
1287  m_tag = "P";
1288  m_dim = 3;
1289  m_taglist = pTagList;
1290  int n = m_conf.m_order-1;
1291 
1292  // This edge-node map is based on Nektar++ ordering.
1293  map<pair<int,int>, int> edgeNodeMap;
1294  map<pair<int,int>, int>::iterator it;
1295  edgeNodeMap[pair<int,int>(1,2)] = 6;
1296  edgeNodeMap[pair<int,int>(2,3)] = 6 + n;
1297  edgeNodeMap[pair<int,int>(4,3)] = 6 + 2*n;
1298  edgeNodeMap[pair<int,int>(1,4)] = 6 + 3*n;
1299  edgeNodeMap[pair<int,int>(1,5)] = 6 + 4*n;
1300  edgeNodeMap[pair<int,int>(2,5)] = 6 + 5*n;
1301  edgeNodeMap[pair<int,int>(3,5)] = 6 + 6*n;
1302  edgeNodeMap[pair<int,int>(4,5)] = 6 + 7*n;
1303 
1304  // Add vertices
1305  for (int i = 0; i < 5; ++i)
1306  {
1307  m_vertex.push_back(pNodeList[i]);
1308  }
1309 
1310  // Create edges (with corresponding set of edge points)
1311  int eid = 0;
1312  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
1313  {
1314  vector<NodeSharedPtr> edgeNodes;
1315  if (m_conf.m_order > 1)
1316  {
1317  for (int j = it->second; j < it->second + n; ++j)
1318  {
1319  edgeNodes.push_back(pNodeList[j-1]);
1320  }
1321  }
1322  m_edge.push_back(
1323  EdgeSharedPtr(new Edge(pNodeList[it->first.first-1],
1324  pNodeList[it->first.second-1],
1325  edgeNodes,
1326  m_conf.m_edgeCurveType)));
1327  m_edge.back()->m_id = eid++;
1328  }
1329 
1330  // Create faces
1331  int face_ids[5][4] = {
1332  {0,1,2,3}, {0,1,4,-1}, {1,2,4,-1}, {3,2,4,-1}, {0,3,4,-1}
1333  };
1334  int face_edges[5][4];
1335  int faceoffset = 5 + 8*n;
1336  for (int j = 0; j < 5; ++j)
1337  {
1338  vector<NodeSharedPtr> faceVertices;
1339  vector<EdgeSharedPtr> faceEdges;
1340  vector<NodeSharedPtr> faceNodes;
1341  int nEdge = j > 0 ? 3 : 4;
1342 
1343  for (int k = 0; k < nEdge; ++k)
1344  {
1345  faceVertices.push_back(m_vertex[face_ids[j][k]]);
1346  NodeSharedPtr a = m_vertex[face_ids[j][k]];
1347  NodeSharedPtr b = m_vertex[face_ids[j][(k+1) % nEdge]];
1348  for (unsigned int i = 0; i < m_edge.size(); ++i)
1349  {
1350  if ((m_edge[i]->m_n1 == a && m_edge[i]->m_n2 == b) ||
1351  (m_edge[i]->m_n1 == b && m_edge[i]->m_n2 == a))
1352  {
1353  faceEdges.push_back(m_edge[i]);
1354  face_edges[j][k] = i;
1355  break;
1356  }
1357  }
1358  }
1359 
1360  if (m_conf.m_faceNodes)
1361  {
1362  int facenodes = j == 0 ? n*n : n*(n-1)/2;
1363  for (int i = 0; i < facenodes; ++i)
1364  {
1365  faceNodes.push_back(pNodeList[faceoffset+i]);
1366  }
1367  faceoffset += facenodes;
1368  }
1369  m_face.push_back(FaceSharedPtr(
1370  new Face(faceVertices, faceNodes, faceEdges, m_conf.m_faceCurveType)));
1371  }
1372 
1373  // Reorder edges to align with Nektar++ order.
1374  vector<EdgeSharedPtr> tmp(8);
1375  tmp[0] = m_edge[face_edges[0][0]];
1376  tmp[1] = m_edge[face_edges[0][1]];
1377  tmp[2] = m_edge[face_edges[0][2]];
1378  tmp[3] = m_edge[face_edges[0][3]];
1379  tmp[4] = m_edge[face_edges[1][2]];
1380  tmp[5] = m_edge[face_edges[1][1]];
1381  tmp[6] = m_edge[face_edges[3][1]];
1382  tmp[7] = m_edge[face_edges[3][2]];
1383  m_edge = tmp;
1384  }
1385 
1386  SpatialDomains::GeometrySharedPtr Pyramid::GetGeom(int coordDim)
1387  {
1388  SpatialDomains::Geometry2DSharedPtr faces[5];
1389 
1390  for (int i = 0; i < 5; ++i)
1391  {
1392  faces[i] = m_face[i]->GetGeom(coordDim);
1393  }
1394 
1395  m_geom = MemoryManager<SpatialDomains::PyrGeom>::
1396  AllocateSharedPtr(faces);
1397 
1398  return m_geom;
1399  }
1400 
1401  /**
1402  * @brief Return the number of nodes defining a pyramid.
1403  */
1404  unsigned int Pyramid::GetNumNodes(ElmtConfig pConf)
1405  {
1406  int n = pConf.m_order;
1407  return 5 + 8*(n-1);
1408  }
1409 
1410  LibUtilities::ShapeType Prism::m_type = GetElementFactory().
1411  RegisterCreatorFunction(LibUtilities::ePrism, Prism::create, "Prism");
1412 
1413  /**
1414  * @brief Create a prism element.
1415  */
1416  Prism::Prism(ElmtConfig pConf,
1417  vector<NodeSharedPtr> pNodeList,
1418  vector<int> pTagList)
1419  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
1420  {
1421  m_tag = "R";
1422  m_dim = 3;
1423  m_taglist = pTagList;
1424  int n = m_conf.m_order-1;
1425 
1426  // Create a map to relate edge nodes to a pair of vertices
1427  // defining an edge. This is based on the ordering produced by
1428  // gmsh.
1429  map<pair<int,int>, int> edgeNodeMap;
1430  map<pair<int,int>, int>::iterator it;
1431 
1432  // This edge-node map is based on Nektar++ ordering.
1433  edgeNodeMap[pair<int,int>(1,2)] = 7;
1434  edgeNodeMap[pair<int,int>(2,3)] = 7 + n;
1435  edgeNodeMap[pair<int,int>(4,3)] = 7 + 2*n;
1436  edgeNodeMap[pair<int,int>(1,4)] = 7 + 3*n;
1437  edgeNodeMap[pair<int,int>(1,5)] = 7 + 4*n;
1438  edgeNodeMap[pair<int,int>(2,5)] = 7 + 5*n;
1439  edgeNodeMap[pair<int,int>(3,6)] = 7 + 6*n;
1440  edgeNodeMap[pair<int,int>(4,6)] = 7 + 7*n;
1441  edgeNodeMap[pair<int,int>(5,6)] = 7 + 8*n;
1442 
1443  // Add vertices
1444  for (int i = 0; i < 6; ++i)
1445  {
1446  m_vertex.push_back(pNodeList[i]);
1447  }
1448 
1449  int eid = 0;
1450  // Create edges (with corresponding set of edge points)
1451  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
1452  {
1453  vector<NodeSharedPtr> edgeNodes;
1454  if (m_conf.m_order > 1) {
1455  for (int j = it->second; j < it->second + n; ++j) {
1456  edgeNodes.push_back(pNodeList[j-1]);
1457  }
1458  }
1459  m_edge.push_back(EdgeSharedPtr(
1460  new Edge(pNodeList[it->first.first-1],
1461  pNodeList[it->first.second-1],
1462  edgeNodes,
1463  m_conf.m_edgeCurveType)));
1464  m_edge.back()->m_id = eid++;
1465  }
1466 
1467  if (m_conf.m_reorient)
1468  {
1469  OrientPrism();
1470  }
1471  else
1472  {
1473  m_orientation = 0;
1474  }
1475 
1476  // Create faces
1477  int face_ids[5][4] = {
1478  {0,1,2,3},{0,1,4,-1},{1,2,5,4},{3,2,5,-1},{0,3,5,4}};
1479  int face_edges[5][4];
1480 
1481  int face_offset[5];
1482  face_offset[0] = 6 + 9*n;
1483  for (int j = 0; j < 4; ++j)
1484  {
1485  int facenodes = j % 2 == 0 ? n*n : n*(n-1)/2;
1486  face_offset[j+1] = face_offset[j] + facenodes;
1487  }
1488 
1489  for (int j = 0; j < 5; ++j)
1490  {
1491  vector<NodeSharedPtr> faceVertices;
1492  vector<EdgeSharedPtr> faceEdges;
1493  vector<NodeSharedPtr> faceNodes;
1494  int nEdge = 3 - (j % 2 - 1);
1495 
1496  for (int k = 0; k < nEdge; ++k)
1497  {
1498  faceVertices.push_back(m_vertex[face_ids[j][k]]);
1499  NodeSharedPtr a = m_vertex[face_ids[j][k]];
1500  NodeSharedPtr b = m_vertex[face_ids[j][(k+1) % nEdge]];
1501  for (unsigned int i = 0; i < m_edge.size(); ++i)
1502  {
1503  if ((m_edge[i]->m_n1 == a && m_edge[i]->m_n2 == b) ||
1504  (m_edge[i]->m_n1 == b && m_edge[i]->m_n2 == a))
1505  {
1506  faceEdges.push_back(m_edge[i]);
1507  face_edges[j][k] = i;
1508  break;
1509  }
1510  }
1511  }
1512 
1513  if (m_conf.m_faceNodes)
1514  {
1515  int face = j, facenodes;
1516 
1517  if (j % 2 == 0)
1518  {
1519  facenodes = n*n;
1520  if (m_orientation == 1)
1521  {
1522  face = (face+4) % 6;
1523  }
1524  else if (m_orientation == 2)
1525  {
1526  face = (face+2) % 6;
1527  }
1528  }
1529  else
1530  {
1531  // TODO: need to rotate these too.
1532  facenodes = n*(n-1)/2;
1533  }
1534 
1535  for (int i = 0; i < facenodes; ++i)
1536  {
1537  faceNodes.push_back(pNodeList[face_offset[face]+i]);
1538  }
1539  }
1540  m_face.push_back(FaceSharedPtr(
1541  new Face(faceVertices, faceNodes, faceEdges, m_conf.m_faceCurveType)));
1542  }
1543 
1544  // Re-order edge array to be consistent with Nektar++ ordering.
1545  vector<EdgeSharedPtr> tmp(9);
1546  tmp[0] = m_edge[face_edges[0][0]];
1547  tmp[1] = m_edge[face_edges[0][1]];
1548  tmp[2] = m_edge[face_edges[0][2]];
1549  tmp[3] = m_edge[face_edges[0][3]];
1550  tmp[4] = m_edge[face_edges[1][2]];
1551  tmp[5] = m_edge[face_edges[1][1]];
1552  tmp[6] = m_edge[face_edges[2][1]];
1553  tmp[7] = m_edge[face_edges[3][2]];
1554  tmp[8] = m_edge[face_edges[4][2]];
1555  m_edge = tmp;
1556  }
1557 
1558  /**
1559  * @brief Return the number of nodes defining a prism.
1560  */
1561  unsigned int Prism::GetNumNodes(ElmtConfig pConf)
1562  {
1563  int n = pConf.m_order;
1564  if (pConf.m_faceNodes && pConf.m_volumeNodes)
1565  return (n+1)*(n+1)*(n+2)/2;
1566  else if (pConf.m_faceNodes && !pConf.m_volumeNodes)
1567  return 3*(n+1)*(n+1)+2*(n+1)*(n+2)/2-9*(n+1)+6;
1568  else
1569  return 9*(n+1)-12;
1570  }
1571 
1572  SpatialDomains::GeometrySharedPtr Prism::GetGeom(int coordDim)
1573  {
1574  SpatialDomains::Geometry2DSharedPtr faces[5];
1575  SpatialDomains::PrismGeomSharedPtr ret;
1576 
1577  for (int i = 0; i < 5; ++i)
1578  {
1579  faces[i] = m_face[i]->GetGeom(coordDim);
1580  }
1581 
1582  ret = MemoryManager<SpatialDomains::PrismGeom>::
1583  AllocateSharedPtr(faces);
1584 
1585  return ret;
1586  }
1587 
1588  /**
1589  * @brief .
1590  */
1591  void Prism::Complete(int order)
1592  {
1593  int i, j, pos;
1594 
1595  // Create basis key for a nodal tetrahedron.
1596  LibUtilities::BasisKey B0(
1597  LibUtilities::eOrtho_A, order+1,
1598  LibUtilities::PointsKey(
1599  order+1,LibUtilities::eGaussLobattoLegendre));
1600  LibUtilities::BasisKey B1(
1601  LibUtilities::eOrtho_A, order+1,
1602  LibUtilities::PointsKey(
1603  order+1,LibUtilities::eGaussLobattoLegendre));
1604  LibUtilities::BasisKey B2(
1605  LibUtilities::eOrtho_B, order+1,
1606  LibUtilities::PointsKey(
1607  order+1,LibUtilities::eGaussRadauMAlpha1Beta0));
1608 
1609  // Create a standard nodal prism in order to get the Vandermonde
1610  // matrix to perform interpolation to nodal points.
1611  StdRegions::StdNodalPrismExpSharedPtr nodalPrism =
1612  MemoryManager<StdRegions::StdNodalPrismExp>::AllocateSharedPtr(
1613  B0, B1, B2, LibUtilities::eNodalPrismEvenlySpaced);
1614 
1615  Array<OneD, NekDouble> x, y, z;
1616  nodalPrism->GetNodalPoints(x,y,z);
1617 
1618  SpatialDomains::PrismGeomSharedPtr geom =
1619  boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(
1620  this->GetGeom(3));
1621 
1622  // Create basis key for a prism.
1623  LibUtilities::BasisKey C0(
1624  LibUtilities::eOrtho_A, order+1,
1625  LibUtilities::PointsKey(
1626  order+1,LibUtilities::eGaussLobattoLegendre));
1627  LibUtilities::BasisKey C1(
1628  LibUtilities::eOrtho_A, order+1,
1629  LibUtilities::PointsKey(
1630  order+1,LibUtilities::eGaussLobattoLegendre));
1631  LibUtilities::BasisKey C2(
1632  LibUtilities::eOrtho_B, order+1,
1633  LibUtilities::PointsKey(
1634  order+1,LibUtilities::eGaussRadauMAlpha1Beta0));
1635 
1636  // Create a prism.
1637  LocalRegions::PrismExpSharedPtr prism =
1638  MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(
1639  C0, C1, C2, geom);
1640 
1641  // Get coordinate array for tetrahedron.
1642  int nqtot = prism->GetTotPoints();
1643  Array<OneD, NekDouble> alloc(6*nqtot);
1644  Array<OneD, NekDouble> xi(alloc);
1645  Array<OneD, NekDouble> yi(alloc+ nqtot);
1646  Array<OneD, NekDouble> zi(alloc+2*nqtot);
1647  Array<OneD, NekDouble> xo(alloc+3*nqtot);
1648  Array<OneD, NekDouble> yo(alloc+4*nqtot);
1649  Array<OneD, NekDouble> zo(alloc+5*nqtot);
1650  Array<OneD, NekDouble> tmp;
1651 
1652  prism->GetCoords(xi, yi, zi);
1653 
1654  for (i = 0; i < 3; ++i)
1655  {
1656  Array<OneD, NekDouble> coeffs(nodalPrism->GetNcoeffs());
1657  prism->FwdTrans(alloc+i*nqtot, coeffs);
1658  // Apply Vandermonde matrix to project onto nodal space.
1659  nodalPrism->ModalToNodal(coeffs, tmp=alloc+(i+3)*nqtot);
1660  }
1661 
1662  // Now extract points from the co-ordinate arrays into the
1663  // edge/face/volume nodes. First, extract edge-interior nodes.
1664  for (i = 0; i < 9; ++i)
1665  {
1666  pos = 6 + i*(order-1);
1667  m_edge[i]->m_edgeNodes.clear();
1668  for (j = 0; j < order-1; ++j)
1669  {
1670  m_edge[i]->m_edgeNodes.push_back(
1671  NodeSharedPtr(new Node(0, xo[pos+j], yo[pos+j], zo[pos+j])));
1672  }
1673  }
1674 
1675  // Now extract face-interior nodes.
1676  pos = 6 + 9*(order-1);
1677  for (i = 0; i < 5; ++i)
1678  {
1679  int facesize = i % 2 ? (order-2)*(order-1)/2 : (order-1)*(order-1);
1680  m_face[i]->m_faceNodes.clear();
1681  for (j = 0; j < facesize; ++j)
1682  {
1683  m_face[i]->m_faceNodes.push_back(
1684  NodeSharedPtr(new Node(0, xo[pos+j], yo[pos+j], zo[pos+j])));
1685  }
1686  pos += facesize;
1687  }
1688 
1689  // Finally extract volume nodes.
1690  for (i = pos; i < (order+1)*(order+1)*(order+2)/2; ++i)
1691  {
1692  m_volumeNodes.push_back(
1693  NodeSharedPtr(new Node(0, xo[i], yo[i], zo[i])));
1694  }
1695 
1696  m_conf.m_order = order;
1697  m_conf.m_faceNodes = true;
1698  m_conf.m_volumeNodes = true;
1699  }
1700 
1701  /**
1702  * @brief Orient prism to align degenerate vertices.
1703  *
1704  * Orientation of prismatric elements is required so that the singular
1705  * vertices of triangular faces (which occur as a part of the
1706  * collapsed co-ordinate system) align. The algorithm is based on that
1707  * used in T. Warburton's thesis and in the original Nektar source.
1708  *
1709  * First the points are re-ordered so that the highest global IDs
1710  * represent the two singular points of the prism. Then, if necessary,
1711  * the nodes are rotated either clockwise or counter-clockwise (w.r.t
1712  * to the p-r plane) to correctly align the prism. The #orientation
1713  * variable is set to:
1714  *
1715  * - 0 if the prism is not rotated;
1716  * - 1 if the prism is rotated clockwise;
1717  * - 2 if the prism is rotated counter-clockwise.
1718  *
1719  * This is necessary for some input modules (e.g. #InputNek) which add
1720  * high-order information or bounary conditions to faces.
1721  */
1722  void Prism::OrientPrism()
1723  {
1724  int lid[6], gid[6];
1725 
1726  // Re-order vertices.
1727  for (int i = 0; i < 6; ++i)
1728  {
1729  lid[i] = i;
1730  gid[i] = m_vertex[i]->m_id;
1731  }
1732 
1733  gid[0] = gid[3] = max(gid[0], gid[3]);
1734  gid[1] = gid[2] = max(gid[1], gid[2]);
1735  gid[4] = gid[5] = max(gid[4], gid[5]);
1736 
1737  for (int i = 1; i < 6; ++i)
1738  {
1739  if (gid[0] < gid[i])
1740  {
1741  swap(gid[i], gid[0]);
1742  swap(lid[i], lid[0]);
1743  }
1744  }
1745 
1746  if (lid[0] == 4 || lid[0] == 5)
1747  {
1748  m_orientation = 0;
1749  }
1750  else if (lid[0] == 1 || lid[0] == 2)
1751  {
1752  // Rotate prism clockwise in p-r plane
1753  vector<NodeSharedPtr> vertexmap(6);
1754  vertexmap[0] = m_vertex[4];
1755  vertexmap[1] = m_vertex[0];
1756  vertexmap[2] = m_vertex[3];
1757  vertexmap[3] = m_vertex[5];
1758  vertexmap[4] = m_vertex[1];
1759  vertexmap[5] = m_vertex[2];
1760  m_vertex = vertexmap;
1761  m_orientation = 1;
1762  }
1763  else if (lid[0] == 0 || lid[0] == 3)
1764  {
1765  // Rotate prism counter-clockwise in p-r plane
1766  vector<NodeSharedPtr> vertexmap(6);
1767  vertexmap[0] = m_vertex[1];
1768  vertexmap[1] = m_vertex[4];
1769  vertexmap[2] = m_vertex[5];
1770  vertexmap[3] = m_vertex[2];
1771  vertexmap[4] = m_vertex[0];
1772  vertexmap[5] = m_vertex[3];
1773  m_vertex = vertexmap;
1774  m_orientation = 2;
1775  }
1776  else
1777  {
1778  cerr << "Warning: possible prism orientation problem." << endl;
1779  }
1780  }
1781 
1782 
1783  LibUtilities::ShapeType Hexahedron::m_type = GetElementFactory().
1784  RegisterCreatorFunction(LibUtilities::eHexahedron, Hexahedron::create, "Hexahedron");
1785 
1786  /**
1787  * @brief Create a hexahedral element.
1788  */
1789  Hexahedron::Hexahedron(ElmtConfig pConf,
1790  vector<NodeSharedPtr> pNodeList,
1791  vector<int> pTagList)
1792  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
1793  {
1794  m_tag = "H";
1795  m_dim = 3;
1796  m_taglist = pTagList;
1797  int n = m_conf.m_order-1;
1798 
1799  // Create a map to relate edge nodes to a pair of vertices defining an edge
1800  // This is based on the ordering produced by gmsh.
1801  map<pair<int,int>, int> edgeNodeMap;
1802  map<pair<int,int>, int>::iterator it;
1803  edgeNodeMap[pair<int,int>(1,2)] = 9;
1804  edgeNodeMap[pair<int,int>(2,3)] = 9 + n;
1805  edgeNodeMap[pair<int,int>(3,4)] = 9 + 2*n;
1806  edgeNodeMap[pair<int,int>(4,1)] = 9 + 3*n;
1807  edgeNodeMap[pair<int,int>(1,5)] = 9 + 4*n;
1808  edgeNodeMap[pair<int,int>(2,6)] = 9 + 5*n;
1809  edgeNodeMap[pair<int,int>(3,7)] = 9 + 6*n;
1810  edgeNodeMap[pair<int,int>(4,8)] = 9 + 7*n;
1811  edgeNodeMap[pair<int,int>(5,6)] = 9 + 8*n;
1812  edgeNodeMap[pair<int,int>(6,7)] = 9 + 9*n;
1813  edgeNodeMap[pair<int,int>(7,8)] = 9 + 10*n;
1814  edgeNodeMap[pair<int,int>(8,5)] = 9 + 11*n;
1815 
1816  // Add vertices
1817  for (int i = 0; i < 8; ++i) {
1818  m_vertex.push_back(pNodeList[i]);
1819  }
1820 
1821  // Create edges (with corresponding set of edge points)
1822  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
1823  {
1824  vector<NodeSharedPtr> edgeNodes;
1825  if (m_conf.m_order > 1) {
1826  for (int j = it->second; j < it->second + n; ++j) {
1827  edgeNodes.push_back(pNodeList[j-1]);
1828  }
1829  }
1830  m_edge.push_back(EdgeSharedPtr(new Edge(pNodeList[it->first.first-1],
1831  pNodeList[it->first.second-1],
1832  edgeNodes,
1833  m_conf.m_edgeCurveType)));
1834  }
1835 
1836  // Create faces
1837  int face_edges[6][4];
1838  int face_ids[6][4] = {
1839  {0,1,2,3},{0,1,5,4},{1,2,6,5},{3,2,6,7},{0,3,7,4},{4,5,6,7}};
1840  for (int j = 0; j < 6; ++j)
1841  {
1842  vector<NodeSharedPtr> faceVertices;
1843  vector<EdgeSharedPtr> faceEdges;
1844  vector<NodeSharedPtr> faceNodes;
1845  for (int k = 0; k < 4; ++k)
1846  {
1847  faceVertices.push_back(m_vertex[face_ids[j][k]]);
1848  NodeSharedPtr a = m_vertex[face_ids[j][k]];
1849  NodeSharedPtr b = m_vertex[face_ids[j][(k+1)%4]];
1850  for (unsigned int i = 0; i < m_edge.size(); ++i)
1851  {
1852  if ( ((*(m_edge[i]->m_n1)==*a) && (*(m_edge[i]->m_n2)==*b))
1853  || ((*(m_edge[i]->m_n1)==*b) && (*(m_edge[i]->m_n2) == *a)) )
1854  {
1855  face_edges[j][k] = i;
1856  faceEdges.push_back(m_edge[i]);
1857  break;
1858  }
1859  }
1860  }
1861 
1862  if (m_conf.m_faceNodes)
1863  {
1864  int N = 8 + 12*n + j*n*n;
1865  for (int i = 0; i < n*n; ++i)
1866  {
1867  faceNodes.push_back(pNodeList[N+i]);
1868  }
1869  }
1870  m_face.push_back(FaceSharedPtr(
1871  new Face(faceVertices, faceNodes, faceEdges, m_conf.m_faceCurveType)));
1872  }
1873 
1874  // Reorder edges to be consistent with Nektar++ ordering.
1875  vector<EdgeSharedPtr> tmp(12);
1876  tmp[ 0] = m_edge[face_edges[0][0]];
1877  tmp[ 1] = m_edge[face_edges[0][1]];
1878  tmp[ 2] = m_edge[face_edges[0][2]];
1879  tmp[ 3] = m_edge[face_edges[0][3]];
1880  tmp[ 4] = m_edge[face_edges[1][3]];
1881  tmp[ 5] = m_edge[face_edges[1][1]];
1882  tmp[ 6] = m_edge[face_edges[2][1]];
1883  tmp[ 7] = m_edge[face_edges[4][1]];
1884  tmp[ 8] = m_edge[face_edges[5][0]];
1885  tmp[ 9] = m_edge[face_edges[5][1]];
1886  tmp[10] = m_edge[face_edges[5][2]];
1887  tmp[11] = m_edge[face_edges[5][3]];
1888  m_edge = tmp;
1889  }
1890 
1891  SpatialDomains::GeometrySharedPtr Hexahedron::GetGeom(int coordDim)
1892  {
1893  SpatialDomains::QuadGeomSharedPtr faces[6];
1894  SpatialDomains::HexGeomSharedPtr ret;
1895 
1896  for (int i = 0; i < 6; ++i)
1897  {
1898  faces[i] = boost::dynamic_pointer_cast
1899  <SpatialDomains::QuadGeom>(m_face[i]->GetGeom(coordDim));
1900  }
1901 
1902  ret = MemoryManager<SpatialDomains::HexGeom>::
1903  AllocateSharedPtr(faces);
1904 
1905  return ret;
1906  }
1907 
1908  /**
1909  * @brief Return the number of nodes defining a hexahedron.
1910  */
1911  unsigned int Hexahedron::GetNumNodes(ElmtConfig pConf)
1912  {
1913  int n = pConf.m_order;
1914  if (pConf.m_faceNodes && pConf.m_volumeNodes)
1915  return (n+1)*(n+1)*(n+1);
1916  else if (pConf.m_faceNodes && !pConf.m_volumeNodes)
1917  return 6*(n+1)*(n+1)-12*(n+1)+8;
1918  else
1919  return 12*(n+1)-16;
1920  }
1921  }
1922 }
Generated on Wed Mar 4 2015 16:33:15 for Nektar++ by   doxygen 1.8.1.2