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  /*
773  cout << "EDGE " << i << " = "
774  << m_edge[i]->m_n1->m_x << "," << m_edge[i]->m_n1->m_y << " "
775  << m_edge[i]->m_n2->m_x << "," << m_edge[i]->m_n2->m_y << endl;
776  */
777  for (int j = 1; j < order; ++j, pos += edgeMap[i][1])
778  {
779  //cout << "INSERTING: " << x[pos] << " " << y[pos] << endl;
780  m_edge[i]->m_edgeNodes.push_back(
781  NodeSharedPtr(new Node(0, x[pos], y[pos], z[pos])));
782  }
783  }
784 
785  // Extract face-interior nodes.
786  m_volumeNodes.clear();
787  for (int i = 1; i < order; ++i)
788  {
789  int pos = i*(order+1);
790  for (int j = 1; j < order; ++j)
791  {
792  m_volumeNodes.push_back(
793  NodeSharedPtr(new Node(0, x[pos+j], y[pos+j], z[pos+j])));
794  //cout << "here" << endl;
795  }
796  }
797 
798  m_conf.m_order = order;
799  m_conf.m_faceNodes = true;
800  m_conf.m_volumeNodes = true;
801  }
802 
803  SpatialDomains::GeometrySharedPtr Quadrilateral::GetGeom(int coordDim)
804  {
805  SpatialDomains::SegGeomSharedPtr edges[4];
806  SpatialDomains::PointGeomSharedPtr verts[4];
807  SpatialDomains::QuadGeomSharedPtr ret;
808 
809  for (int i = 0; i < 4; ++i)
810  {
811  edges[i] = m_edge [i]->GetGeom(coordDim);
812  verts[i] = m_vertex[i]->GetGeom(coordDim);
813  }
814 
815  StdRegions::Orientation edgeorient[4] = {
816  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[0], *edges[1]),
817  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[1], *edges[2]),
818  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[2], *edges[3]),
819  SpatialDomains::SegGeom::GetEdgeOrientation(*edges[3], *edges[0])
820  };
821 
822  ret = MemoryManager<SpatialDomains::QuadGeom>::
823  AllocateSharedPtr(m_id, verts, edges, edgeorient);
824 
825  return ret;
826  }
827 
828  /**
829  * @brief Return the number of nodes defining a quadrilateral.
830  */
831  unsigned int Quadrilateral::GetNumNodes(ElmtConfig pConf)
832  {
833  int n = pConf.m_order;
834  if (!pConf.m_faceNodes)
835  return 4*n;
836  else
837  return (n+1)*(n+1);
838  }
839 
840 
841  LibUtilities::ShapeType Tetrahedron::m_type = GetElementFactory().
842  RegisterCreatorFunction(LibUtilities::eTetrahedron, Tetrahedron::create, "Tetrahedron");
843 
844  /**
845  * @brief Create a tetrahedron element.
846  */
847  Tetrahedron::Tetrahedron(ElmtConfig pConf,
848  vector<NodeSharedPtr> pNodeList,
849  vector<int> pTagList)
850  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
851  {
852  m_tag = "A";
853  m_dim = 3;
854  m_taglist = pTagList;
855  int n = m_conf.m_order-1;
856 
857  // Create a map to relate edge nodes to a pair of vertices
858  // defining an edge.
859  map<pair<int,int>, int> edgeNodeMap;
860  map<pair<int,int>, int>::iterator it;
861  edgeNodeMap[pair<int,int>(1,2)] = 5;
862  edgeNodeMap[pair<int,int>(2,3)] = 5 + n;
863  edgeNodeMap[pair<int,int>(1,3)] = 5 + 2*n;
864  edgeNodeMap[pair<int,int>(1,4)] = 5 + 3*n;
865  edgeNodeMap[pair<int,int>(2,4)] = 5 + 4*n;
866  edgeNodeMap[pair<int,int>(3,4)] = 5 + 5*n;
867 
868  // Add vertices
869  for (int i = 0; i < 4; ++i) {
870  m_vertex.push_back(pNodeList[i]);
871  }
872 
873  // Create edges (with corresponding set of edge points)
874  int eid = 0;
875  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
876  {
877  vector<NodeSharedPtr> edgeNodes;
878  if (m_conf.m_order > 1) {
879  for (int j = it->second; j < it->second + n; ++j) {
880  edgeNodes.push_back(pNodeList[j-1]);
881  }
882  }
883  m_edge.push_back(EdgeSharedPtr(new Edge(pNodeList[it->first.first-1],
884  pNodeList[it->first.second-1],
885  edgeNodes,
886  m_conf.m_edgeCurveType)));
887  m_edge.back()->m_id = eid++;
888  }
889 
890  // Reorient the tet to ensure collapsed coordinates align between adjacent
891  // elements.
892  if (m_conf.m_reorient)
893  {
894  OrientTet();
895  }
896 
897  // Create faces
898  int face_ids[4][3] = {
899  {0,1,2},{0,1,3},{1,2,3},{0,2,3}};
900  int face_edges[4][3];
901 
902  for (int j = 0; j < 4; ++j)
903  {
904  vector<NodeSharedPtr> faceVertices;
905  vector<EdgeSharedPtr> faceEdges;
906  vector<NodeSharedPtr> faceNodes;
907  for (int k = 0; k < 3; ++k)
908  {
909  faceVertices.push_back(m_vertex[face_ids[j][k]]);
910  NodeSharedPtr a = m_vertex[face_ids[j][k]];
911  NodeSharedPtr b = m_vertex[face_ids[j][(k+1)%3]];
912  for (unsigned int i = 0; i < m_edge.size(); ++i)
913  {
914  if ( ((*(m_edge[i]->m_n1)==*a) && (*(m_edge[i]->m_n2)==*b))
915  || ((*(m_edge[i]->m_n1)==*b) && (*(m_edge[i]->m_n2) == *a)) )
916  {
917  face_edges[j][k] = i;
918  faceEdges.push_back(m_edge[i]);
919  break;
920  }
921  }
922  }
923 
924  if (m_conf.m_faceNodes)
925  {
926  int N = 4 + 6*n + j*n*(n-1)/2;
927  for (int i = 0; i < n*(n-1)/2; ++i)
928  {
929  faceNodes.push_back(pNodeList[N+i]);
930  }
931  }
932  m_face.push_back(FaceSharedPtr(
933  new Face(faceVertices, faceNodes, faceEdges, m_conf.m_faceCurveType)));
934  }
935 
936  vector<EdgeSharedPtr> tmp(6);
937  tmp[0] = m_edge[face_edges[0][0]];
938  tmp[1] = m_edge[face_edges[0][1]];
939  tmp[2] = m_edge[face_edges[0][2]];
940  tmp[3] = m_edge[face_edges[1][2]];
941  tmp[4] = m_edge[face_edges[1][1]];
942  tmp[5] = m_edge[face_edges[2][1]];
943  m_edge = tmp;
944  }
945 
946  SpatialDomains::GeometrySharedPtr Tetrahedron::GetGeom(int coordDim)
947  {
948  SpatialDomains::TriGeomSharedPtr tfaces[4];
949  SpatialDomains::TetGeomSharedPtr ret;
950 
951  for (int i = 0; i < 4; ++i)
952  {
953  tfaces[i] = boost::dynamic_pointer_cast
954  <SpatialDomains::TriGeom>(m_face[i]->GetGeom(coordDim));
955  }
956 
957  ret = MemoryManager<SpatialDomains::TetGeom>::
958  AllocateSharedPtr(tfaces);
959 
960  return ret;
961  }
962 
963  /**
964  * @brief Return the number of nodes defining a tetrahedron.
965  */
966  unsigned int Tetrahedron::GetNumNodes(ElmtConfig pConf)
967  {
968  int n = pConf.m_order;
969  if (pConf.m_volumeNodes && pConf.m_faceNodes)
970  return (n+1)*(n+2)*(n+3)/6;
971  else if (!pConf.m_volumeNodes && pConf.m_faceNodes)
972  return 4*(n+1)*(n+2)/2-6*(n+1)+4;
973  else
974  return 6*(n+1)-8;
975  }
976 
977  /**
978  * @brief .
979  */
980  void Tetrahedron::Complete(int order)
981  {
982  int i, j;
983 
984  // Create basis key for a nodal tetrahedron.
985  LibUtilities::BasisKey B0(
986  LibUtilities::eOrtho_A, order+1,
987  LibUtilities::PointsKey(
988  order+1,LibUtilities::eGaussLobattoLegendre));
989  LibUtilities::BasisKey B1(
990  LibUtilities::eOrtho_B, order+1,
991  LibUtilities::PointsKey(
992  order+1,LibUtilities::eGaussRadauMAlpha1Beta0));
993  LibUtilities::BasisKey B2(
994  LibUtilities::eOrtho_C, order+1,
995  LibUtilities::PointsKey(
996  order+1,LibUtilities::eGaussRadauMAlpha2Beta0));
997 
998  // Create a standard nodal tetrahedron in order to get the
999  // Vandermonde matrix to perform interpolation to nodal points.
1000  StdRegions::StdNodalTetExpSharedPtr nodalTet =
1001  MemoryManager<StdRegions::StdNodalTetExp>::AllocateSharedPtr(
1002  B0, B1, B2, LibUtilities::eNodalTetEvenlySpaced);
1003 
1004  Array<OneD, NekDouble> x, y, z;
1005 
1006  nodalTet->GetNodalPoints(x,y,z);
1007 
1008  SpatialDomains::TetGeomSharedPtr geom =
1009  boost::dynamic_pointer_cast<SpatialDomains::TetGeom>(
1010  this->GetGeom(3));
1011 
1012  // Create basis key for a tetrahedron.
1013  LibUtilities::BasisKey C0(
1014  LibUtilities::eOrtho_A, order+1,
1015  LibUtilities::PointsKey(
1016  order+1,LibUtilities::eGaussLobattoLegendre));
1017  LibUtilities::BasisKey C1(
1018  LibUtilities::eOrtho_B, order+1,
1019  LibUtilities::PointsKey(
1020  order+1,LibUtilities::eGaussRadauMAlpha1Beta0));
1021  LibUtilities::BasisKey C2(
1022  LibUtilities::eOrtho_C, order+1,
1023  LibUtilities::PointsKey(
1024  order+1,LibUtilities::eGaussRadauMAlpha2Beta0));
1025 
1026  // Create a tet.
1027  LocalRegions::TetExpSharedPtr tet =
1028  MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(
1029  C0, C1, C2, geom);
1030 
1031  // Get coordinate array for tetrahedron.
1032  int nqtot = tet->GetTotPoints();
1033  Array<OneD, NekDouble> alloc(6*nqtot);
1034  Array<OneD, NekDouble> xi(alloc);
1035  Array<OneD, NekDouble> yi(alloc+ nqtot);
1036  Array<OneD, NekDouble> zi(alloc+2*nqtot);
1037  Array<OneD, NekDouble> xo(alloc+3*nqtot);
1038  Array<OneD, NekDouble> yo(alloc+4*nqtot);
1039  Array<OneD, NekDouble> zo(alloc+5*nqtot);
1040  Array<OneD, NekDouble> tmp;
1041 
1042  tet->GetCoords(xi, yi, zi);
1043 
1044  for (i = 0; i < 3; ++i)
1045  {
1046  Array<OneD, NekDouble> coeffs(nodalTet->GetNcoeffs());
1047  tet->FwdTrans(alloc+i*nqtot, coeffs);
1048  // Apply Vandermonde matrix to project onto nodal space.
1049  nodalTet->ModalToNodal(coeffs, tmp=alloc+(i+3)*nqtot);
1050  }
1051 
1052  // Now extract points from the co-ordinate arrays into the
1053  // edge/face/volume nodes. First, extract edge-interior nodes.
1054  for (i = 0; i < 6; ++i)
1055  {
1056  int pos = 4 + i*(order-1);
1057  m_edge[i]->m_edgeNodes.clear();
1058  for (j = 0; j < order-1; ++j)
1059  {
1060  m_edge[i]->m_edgeNodes.push_back(
1061  NodeSharedPtr(new Node(0, xo[pos+j], yo[pos+j], zo[pos+j])));
1062  }
1063  }
1064 
1065  // Now extract face-interior nodes.
1066  for (i = 0; i < 4; ++i)
1067  {
1068  int pos = 4 + 6*(order-1) + i*(order-2)*(order-1)/2;
1069  m_face[i]->m_faceNodes.clear();
1070  for (j = 0; j < (order-2)*(order-1)/2; ++j)
1071  {
1072  m_face[i]->m_faceNodes.push_back(
1073  NodeSharedPtr(new Node(0, xo[pos+j], yo[pos+j], zo[pos+j])));
1074  }
1075  }
1076 
1077  // Finally extract volume nodes.
1078  int pos = 4 + 6*(order-1) + 4*(order-2)*(order-1)/2;
1079  for (i = pos; i < (order+1)*(order+2)*(order+3)/6; ++i)
1080  {
1081  m_volumeNodes.push_back(
1082  NodeSharedPtr(new Node(0, xo[i], yo[i], zo[i])));
1083  }
1084 
1085  m_conf.m_order = order;
1086  m_conf.m_faceNodes = true;
1087  m_conf.m_volumeNodes = true;
1088  }
1089 
1090  struct TetOrient
1091  {
1092  TetOrient(vector<int> nid, int fid) : nid(nid), fid(fid) {}
1093  vector<int> nid;
1094  int fid;
1095  };
1096 
1097  struct TetOrientHash : std::unary_function<struct TetOrient, std::size_t>
1098  {
1099  std::size_t operator()(struct TetOrient const& p) const
1100  {
1101  return boost::hash_range(p.nid.begin(), p.nid.end());
1102  }
1103  };
1104  typedef boost::unordered_set<struct TetOrient, TetOrientHash> TetOrientSet;
1105 
1106  bool operator==(const struct TetOrient &a, const struct TetOrient &b)
1107  {
1108  if (a.nid.size() != b.nid.size())
1109  {
1110  return false;
1111  }
1112 
1113  for (int i = 0; i < a.nid.size(); ++i)
1114  {
1115  if (a.nid[i] != b.nid[i])
1116  {
1117  return false;
1118  }
1119  }
1120 
1121  return true;
1122  }
1123 
1124  /**
1125  * @brief Orient tetrahedron to align degenerate vertices.
1126  *
1127  * Orientation of tetrahedral elements is required so that the
1128  * singular vertices of triangular faces (which occur as a part of the
1129  * collapsed co-ordinate system) align. The algorithm is based on that
1130  * used in T. Warburton's thesis and in the original Nektar source.
1131  *
1132  * First the vertices are ordered with the highest global vertex at
1133  * the top degenerate point, and the base degenerate point has second
1134  * lowest ID. These vertices are swapped if the element is incorrectly
1135  * oriented.
1136  */
1137  void Tetrahedron::OrientTet()
1138  {
1139  static int face_ids[4][3] = {
1140  {0,1,2},{0,1,3},{1,2,3},{0,2,3}};
1141 
1142  TetOrientSet faces;
1143 
1144  // Create a copy of the original vertex ordering. This is used to
1145  // construct a mapping, #orientationMap, which maps the original
1146  // face ordering to the new face ordering.
1147  for (int i = 0; i < 4; ++i)
1148  {
1149  vector<int> nodes(3);
1150 
1151  nodes[0] = m_vertex[face_ids[i][0]]->m_id;
1152  nodes[1] = m_vertex[face_ids[i][1]]->m_id;
1153  nodes[2] = m_vertex[face_ids[i][2]]->m_id;
1154 
1155  sort(nodes.begin(), nodes.end());
1156  struct TetOrient faceNodes(nodes, i);
1157  faces.insert(faceNodes);
1158  }
1159 
1160  // Order vertices with highest global vertex at top degenerate
1161  // point. Place second highest global vertex at base degenerate
1162  // point.
1163  sort(m_vertex.begin(), m_vertex.end());
1164 
1165  // Calculate a.(b x c) to determine tet volume; if negative,
1166  // reverse order of non-degenerate points to correctly orientate
1167  // the tet.
1168  double ax = m_vertex[1]->m_x-m_vertex[0]->m_x;
1169  double ay = m_vertex[1]->m_y-m_vertex[0]->m_y;
1170  double az = m_vertex[1]->m_z-m_vertex[0]->m_z;
1171  double bx = m_vertex[2]->m_x-m_vertex[0]->m_x;
1172  double by = m_vertex[2]->m_y-m_vertex[0]->m_y;
1173  double bz = m_vertex[2]->m_z-m_vertex[0]->m_z;
1174  double cx = m_vertex[3]->m_x-m_vertex[0]->m_x;
1175  double cy = m_vertex[3]->m_y-m_vertex[0]->m_y;
1176  double cz = m_vertex[3]->m_z-m_vertex[0]->m_z;
1177  double vol = cx*(ay*bz-az*by)+cy*(az*bx-ax*bz)+cz*(ax*by-ay*bx);
1178  vol /= 6.0;
1179 
1180  if (fabs(vol) <= 1e-10)
1181  {
1182  cerr << "Warning: degenerate tetrahedron, volume = " << vol << endl;
1183  }
1184 
1185  if (vol < 0)
1186  {
1187  swap(m_vertex[0], m_vertex[1]);
1188  }
1189 
1190  TetOrientSet::iterator it;
1191 
1192  // Search for the face in the original set of face nodes. Then use
1193  // this to construct the #orientationMap.
1194  for (int i = 0; i < 4; ++i)
1195  {
1196  vector<int> nodes(3);
1197 
1198  nodes[0] = m_vertex[face_ids[i][0]]->m_id;
1199  nodes[1] = m_vertex[face_ids[i][1]]->m_id;
1200  nodes[2] = m_vertex[face_ids[i][2]]->m_id;
1201  sort(nodes.begin(), nodes.end());
1202 
1203  struct TetOrient faceNodes(nodes, 0);
1204 
1205  it = faces.find(faceNodes);
1206  m_orientationMap[it->fid] = i;
1207  }
1208  }
1209 
1210  LibUtilities::ShapeType Pyramid::type = GetElementFactory().
1211  RegisterCreatorFunction(LibUtilities::ePyramid, Pyramid::create, "Pyramid");
1212 
1213  /**
1214  * @brief Create a pyramidic element.
1215  */
1216  Pyramid::Pyramid(ElmtConfig pConf,
1217  vector<NodeSharedPtr> pNodeList,
1218  vector<int> pTagList)
1219  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
1220  {
1221  m_tag = "P";
1222  m_dim = 3;
1223  m_taglist = pTagList;
1224  int n = m_conf.m_order-1;
1225 
1226  // This edge-node map is based on Nektar++ ordering.
1227  map<pair<int,int>, int> edgeNodeMap;
1228  map<pair<int,int>, int>::iterator it;
1229  edgeNodeMap[pair<int,int>(1,2)] = 6;
1230  edgeNodeMap[pair<int,int>(2,3)] = 6 + n;
1231  edgeNodeMap[pair<int,int>(4,3)] = 6 + 2*n;
1232  edgeNodeMap[pair<int,int>(1,4)] = 6 + 3*n;
1233  edgeNodeMap[pair<int,int>(1,5)] = 6 + 4*n;
1234  edgeNodeMap[pair<int,int>(2,5)] = 6 + 5*n;
1235  edgeNodeMap[pair<int,int>(3,5)] = 6 + 6*n;
1236  edgeNodeMap[pair<int,int>(4,5)] = 6 + 7*n;
1237 
1238  // Add vertices
1239  for (int i = 0; i < 5; ++i)
1240  {
1241  m_vertex.push_back(pNodeList[i]);
1242  }
1243 
1244  // Create edges (with corresponding set of edge points)
1245  int eid = 0;
1246  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
1247  {
1248  vector<NodeSharedPtr> edgeNodes;
1249  if (m_conf.m_order > 1)
1250  {
1251  for (int j = it->second; j < it->second + n; ++j)
1252  {
1253  edgeNodes.push_back(pNodeList[j-1]);
1254  }
1255  }
1256  m_edge.push_back(
1257  EdgeSharedPtr(new Edge(pNodeList[it->first.first-1],
1258  pNodeList[it->first.second-1],
1259  edgeNodes,
1260  m_conf.m_edgeCurveType)));
1261  m_edge.back()->m_id = eid++;
1262  }
1263 
1264  // Create faces
1265  int face_ids[5][4] = {
1266  {0,1,2,3}, {0,1,4,-1}, {1,2,4,-1}, {3,2,4,-1}, {0,3,4,-1}
1267  };
1268  int face_edges[5][4];
1269  int faceoffset = 0;
1270  for (int j = 0; j < 5; ++j)
1271  {
1272  vector<NodeSharedPtr> faceVertices;
1273  vector<EdgeSharedPtr> faceEdges;
1274  vector<NodeSharedPtr> faceNodes;
1275  int nEdge = j > 0 ? 3 : 4;
1276 
1277  for (int k = 0; k < nEdge; ++k)
1278  {
1279  faceVertices.push_back(m_vertex[face_ids[j][k]]);
1280  NodeSharedPtr a = m_vertex[face_ids[j][k]];
1281  NodeSharedPtr b = m_vertex[face_ids[j][(k+1) % nEdge]];
1282  for (unsigned int i = 0; i < m_edge.size(); ++i)
1283  {
1284  if ((m_edge[i]->m_n1 == a && m_edge[i]->m_n2 == b) ||
1285  (m_edge[i]->m_n1 == b && m_edge[i]->m_n2 == a))
1286  {
1287  faceEdges.push_back(m_edge[i]);
1288  face_edges[j][k] = i;
1289  break;
1290  }
1291  }
1292  }
1293 
1294  if (m_conf.m_faceNodes)
1295  {
1296  int facenodes = j == 0 ? n*n : n*(n-1)/2;
1297  faceoffset += facenodes;
1298  int N = 6 + 9*n + faceoffset;
1299  for (int i = 0; i < facenodes; ++i)
1300  {
1301  faceNodes.push_back(pNodeList[N+i]);
1302  }
1303  }
1304  m_face.push_back(FaceSharedPtr(
1305  new Face(faceVertices, faceNodes, faceEdges, m_conf.m_faceCurveType)));
1306  }
1307 
1308  vector<EdgeSharedPtr> tmp(8);
1309  tmp[0] = m_edge[face_edges[0][0]];
1310  tmp[1] = m_edge[face_edges[0][1]];
1311  tmp[2] = m_edge[face_edges[0][2]];
1312  tmp[3] = m_edge[face_edges[0][3]];
1313  tmp[4] = m_edge[face_edges[1][2]];
1314  tmp[5] = m_edge[face_edges[1][1]];
1315  tmp[6] = m_edge[face_edges[3][1]];
1316  tmp[7] = m_edge[face_edges[3][2]];
1317  m_edge = tmp;
1318  }
1319 
1320  SpatialDomains::GeometrySharedPtr Pyramid::GetGeom(int coordDim)
1321  {
1322  SpatialDomains::Geometry2DSharedPtr faces[5];
1323 
1324  for (int i = 0; i < 5; ++i)
1325  {
1326  faces[i] = m_face[i]->GetGeom(coordDim);
1327  }
1328 
1329  m_geom = MemoryManager<SpatialDomains::PyrGeom>::
1330  AllocateSharedPtr(faces);
1331 
1332  return m_geom;
1333  }
1334 
1335  /**
1336  * @brief Return the number of nodes defining a pyramid.
1337  */
1338  unsigned int Pyramid::GetNumNodes(ElmtConfig pConf)
1339  {
1340  int n = pConf.m_order;
1341  return 5 + 8*(n-1);
1342  }
1343 
1344  LibUtilities::ShapeType Prism::m_type = GetElementFactory().
1345  RegisterCreatorFunction(LibUtilities::ePrism, Prism::create, "Prism");
1346 
1347  /**
1348  * @brief Create a prism element.
1349  */
1350  Prism::Prism(ElmtConfig pConf,
1351  vector<NodeSharedPtr> pNodeList,
1352  vector<int> pTagList)
1353  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
1354  {
1355  m_tag = "R";
1356  m_dim = 3;
1357  m_taglist = pTagList;
1358  int n = m_conf.m_order-1;
1359 
1360  // Create a map to relate edge nodes to a pair of vertices
1361  // defining an edge. This is based on the ordering produced by
1362  // gmsh.
1363  map<pair<int,int>, int> edgeNodeMap;
1364  map<pair<int,int>, int>::iterator it;
1365 
1366  // This edge-node map is based on Nektar++ ordering.
1367  edgeNodeMap[pair<int,int>(1,2)] = 7;
1368  edgeNodeMap[pair<int,int>(2,3)] = 7 + n;
1369  edgeNodeMap[pair<int,int>(4,3)] = 7 + 2*n;
1370  edgeNodeMap[pair<int,int>(1,4)] = 7 + 3*n;
1371  edgeNodeMap[pair<int,int>(1,5)] = 7 + 4*n;
1372  edgeNodeMap[pair<int,int>(2,5)] = 7 + 5*n;
1373  edgeNodeMap[pair<int,int>(3,6)] = 7 + 6*n;
1374  edgeNodeMap[pair<int,int>(4,6)] = 7 + 7*n;
1375  edgeNodeMap[pair<int,int>(5,6)] = 7 + 8*n;
1376 
1377  // Add vertices
1378  for (int i = 0; i < 6; ++i)
1379  {
1380  m_vertex.push_back(pNodeList[i]);
1381  }
1382 
1383  int eid = 0;
1384  // Create edges (with corresponding set of edge points)
1385  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
1386  {
1387  vector<NodeSharedPtr> edgeNodes;
1388  if (m_conf.m_order > 1) {
1389  for (int j = it->second; j < it->second + n; ++j) {
1390  edgeNodes.push_back(pNodeList[j-1]);
1391  }
1392  }
1393  m_edge.push_back(EdgeSharedPtr(
1394  new Edge(pNodeList[it->first.first-1],
1395  pNodeList[it->first.second-1],
1396  edgeNodes,
1397  m_conf.m_edgeCurveType)));
1398  m_edge.back()->m_id = eid++;
1399  }
1400 
1401  if (m_conf.m_reorient)
1402  {
1403  OrientPrism();
1404  }
1405 
1406  // Create faces
1407  int face_ids[5][4] = {
1408  {0,1,2,3},{0,1,4,-1},{1,2,5,4},{3,2,5,-1},{0,3,5,4}};
1409  int face_edges[5][4];
1410  int faceoffset = 0;
1411  for (int j = 0; j < 5; ++j)
1412  {
1413  vector<NodeSharedPtr> faceVertices;
1414  vector<EdgeSharedPtr> faceEdges;
1415  vector<NodeSharedPtr> faceNodes;
1416  int nEdge = 3 - (j % 2 - 1);
1417 
1418  for (int k = 0; k < nEdge; ++k)
1419  {
1420  faceVertices.push_back(m_vertex[face_ids[j][k]]);
1421  NodeSharedPtr a = m_vertex[face_ids[j][k]];
1422  NodeSharedPtr b = m_vertex[face_ids[j][(k+1) % nEdge]];
1423  for (unsigned int i = 0; i < m_edge.size(); ++i)
1424  {
1425  if ((m_edge[i]->m_n1 == a && m_edge[i]->m_n2 == b) ||
1426  (m_edge[i]->m_n1 == b && m_edge[i]->m_n2 == a))
1427  {
1428  faceEdges.push_back(m_edge[i]);
1429  face_edges[j][k] = i;
1430  break;
1431  }
1432  }
1433  }
1434 
1435  if (m_conf.m_faceNodes)
1436  {
1437  int facenodes = j%2==0 ? n*n : n*(n-1)/2;
1438  faceoffset += facenodes;
1439  int N = 6 + 9*n + faceoffset;
1440  for (int i = 0; i < facenodes; ++i)
1441  {
1442  faceNodes.push_back(pNodeList[N+i]);
1443  }
1444  }
1445  m_face.push_back(FaceSharedPtr(
1446  new Face(faceVertices, faceNodes, faceEdges, m_conf.m_faceCurveType)));
1447  }
1448 
1449  // Re-order edge array to be consistent with Nektar++ ordering.
1450  vector<EdgeSharedPtr> tmp(9);
1451  tmp[0] = m_edge[face_edges[0][0]];
1452  tmp[1] = m_edge[face_edges[0][1]];
1453  tmp[2] = m_edge[face_edges[0][2]];
1454  tmp[3] = m_edge[face_edges[0][3]];
1455  tmp[4] = m_edge[face_edges[1][2]];
1456  tmp[5] = m_edge[face_edges[1][1]];
1457  tmp[6] = m_edge[face_edges[2][1]];
1458  tmp[7] = m_edge[face_edges[3][2]];
1459  tmp[8] = m_edge[face_edges[4][2]];
1460  m_edge = tmp;
1461  }
1462 
1463  /**
1464  * @brief Return the number of nodes defining a prism.
1465  */
1466  unsigned int Prism::GetNumNodes(ElmtConfig pConf)
1467  {
1468  int n = pConf.m_order;
1469  if (pConf.m_faceNodes && pConf.m_volumeNodes)
1470  return (n+1)*(n+1)*(n+2)/2;
1471  else if (pConf.m_faceNodes && !pConf.m_volumeNodes)
1472  return 3*(n+1)*(n+1)+2*(n+1)*(n+2)/2-9*(n+1)+6;
1473  else
1474  return 9*(n+1)-12;
1475  }
1476 
1477  SpatialDomains::GeometrySharedPtr Prism::GetGeom(int coordDim)
1478  {
1479  SpatialDomains::Geometry2DSharedPtr faces[5];
1480  SpatialDomains::PrismGeomSharedPtr ret;
1481 
1482  for (int i = 0; i < 5; ++i)
1483  {
1484  faces[i] = m_face[i]->GetGeom(coordDim);
1485  }
1486 
1487  ret = MemoryManager<SpatialDomains::PrismGeom>::
1488  AllocateSharedPtr(faces);
1489 
1490  return ret;
1491  }
1492 
1493  /**
1494  * @brief .
1495  */
1496  void Prism::Complete(int order)
1497  {
1498  int i, j, pos;
1499 
1500  // Create basis key for a nodal tetrahedron.
1501  LibUtilities::BasisKey B0(
1502  LibUtilities::eOrtho_A, order+1,
1503  LibUtilities::PointsKey(
1504  order+1,LibUtilities::eGaussLobattoLegendre));
1505  LibUtilities::BasisKey B1(
1506  LibUtilities::eOrtho_A, order+1,
1507  LibUtilities::PointsKey(
1508  order+1,LibUtilities::eGaussLobattoLegendre));
1509  LibUtilities::BasisKey B2(
1510  LibUtilities::eOrtho_B, order+1,
1511  LibUtilities::PointsKey(
1512  order+1,LibUtilities::eGaussRadauMAlpha1Beta0));
1513 
1514  // Create a standard nodal prism in order to get the Vandermonde
1515  // matrix to perform interpolation to nodal points.
1516  StdRegions::StdNodalPrismExpSharedPtr nodalPrism =
1517  MemoryManager<StdRegions::StdNodalPrismExp>::AllocateSharedPtr(
1518  B0, B1, B2, LibUtilities::eNodalPrismEvenlySpaced);
1519 
1520  Array<OneD, NekDouble> x, y, z;
1521  nodalPrism->GetNodalPoints(x,y,z);
1522 
1523  SpatialDomains::PrismGeomSharedPtr geom =
1524  boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(
1525  this->GetGeom(3));
1526 
1527  // Create basis key for a prism.
1528  LibUtilities::BasisKey C0(
1529  LibUtilities::eOrtho_A, order+1,
1530  LibUtilities::PointsKey(
1531  order+1,LibUtilities::eGaussLobattoLegendre));
1532  LibUtilities::BasisKey C1(
1533  LibUtilities::eOrtho_A, order+1,
1534  LibUtilities::PointsKey(
1535  order+1,LibUtilities::eGaussLobattoLegendre));
1536  LibUtilities::BasisKey C2(
1537  LibUtilities::eOrtho_B, order+1,
1538  LibUtilities::PointsKey(
1539  order+1,LibUtilities::eGaussRadauMAlpha1Beta0));
1540 
1541  // Create a prism.
1542  LocalRegions::PrismExpSharedPtr prism =
1543  MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(
1544  C0, C1, C2, geom);
1545 
1546  // Get coordinate array for tetrahedron.
1547  int nqtot = prism->GetTotPoints();
1548  Array<OneD, NekDouble> alloc(6*nqtot);
1549  Array<OneD, NekDouble> xi(alloc);
1550  Array<OneD, NekDouble> yi(alloc+ nqtot);
1551  Array<OneD, NekDouble> zi(alloc+2*nqtot);
1552  Array<OneD, NekDouble> xo(alloc+3*nqtot);
1553  Array<OneD, NekDouble> yo(alloc+4*nqtot);
1554  Array<OneD, NekDouble> zo(alloc+5*nqtot);
1555  Array<OneD, NekDouble> tmp;
1556 
1557  prism->GetCoords(xi, yi, zi);
1558 
1559  for (i = 0; i < 3; ++i)
1560  {
1561  Array<OneD, NekDouble> coeffs(nodalPrism->GetNcoeffs());
1562  prism->FwdTrans(alloc+i*nqtot, coeffs);
1563  // Apply Vandermonde matrix to project onto nodal space.
1564  nodalPrism->ModalToNodal(coeffs, tmp=alloc+(i+3)*nqtot);
1565  }
1566 
1567  // Now extract points from the co-ordinate arrays into the
1568  // edge/face/volume nodes. First, extract edge-interior nodes.
1569  for (i = 0; i < 9; ++i)
1570  {
1571  pos = 6 + i*(order-1);
1572  m_edge[i]->m_edgeNodes.clear();
1573  for (j = 0; j < order-1; ++j)
1574  {
1575  m_edge[i]->m_edgeNodes.push_back(
1576  NodeSharedPtr(new Node(0, xo[pos+j], yo[pos+j], zo[pos+j])));
1577  }
1578  }
1579 
1580  // Now extract face-interior nodes.
1581  pos = 6 + 9*(order-1);
1582  for (i = 0; i < 5; ++i)
1583  {
1584  int facesize = i % 2 ? (order-2)*(order-1)/2 : (order-1)*(order-1);
1585  m_face[i]->m_faceNodes.clear();
1586  for (j = 0; j < facesize; ++j)
1587  {
1588  m_face[i]->m_faceNodes.push_back(
1589  NodeSharedPtr(new Node(0, xo[pos+j], yo[pos+j], zo[pos+j])));
1590  }
1591  pos += facesize;
1592  }
1593 
1594  // Finally extract volume nodes.
1595  for (i = pos; i < (order+1)*(order+1)*(order+2)/2; ++i)
1596  {
1597  m_volumeNodes.push_back(
1598  NodeSharedPtr(new Node(0, xo[i], yo[i], zo[i])));
1599  }
1600 
1601  m_conf.m_order = order;
1602  m_conf.m_faceNodes = true;
1603  m_conf.m_volumeNodes = true;
1604  }
1605 
1606  /**
1607  * @brief Orient prism to align degenerate vertices.
1608  *
1609  * Orientation of prismatric elements is required so that the singular
1610  * vertices of triangular faces (which occur as a part of the
1611  * collapsed co-ordinate system) align. The algorithm is based on that
1612  * used in T. Warburton's thesis and in the original Nektar source.
1613  *
1614  * First the points are re-ordered so that the highest global IDs
1615  * represent the two singular points of the prism. Then, if necessary,
1616  * the nodes are rotated either clockwise or counter-clockwise (w.r.t
1617  * to the p-r plane) to correctly align the prism. The #orientation
1618  * variable is set to:
1619  *
1620  * - 0 if the prism is not rotated;
1621  * - 1 if the prism is rotated clockwise;
1622  * - 2 if the prism is rotated counter-clockwise.
1623  *
1624  * This is necessary for some input modules (e.g. #InputNek) which add
1625  * high-order information or bounary conditions to faces.
1626  */
1627  void Prism::OrientPrism()
1628  {
1629  int lid[6], gid[6];
1630 
1631  // Re-order vertices.
1632  for (int i = 0; i < 6; ++i)
1633  {
1634  lid[i] = i;
1635  gid[i] = m_vertex[i]->m_id;
1636  }
1637 
1638  gid[0] = gid[3] = max(gid[0], gid[3]);
1639  gid[1] = gid[2] = max(gid[1], gid[2]);
1640  gid[4] = gid[5] = max(gid[4], gid[5]);
1641 
1642  for (int i = 1; i < 6; ++i)
1643  {
1644  if (gid[0] < gid[i])
1645  {
1646  swap(gid[i], gid[0]);
1647  swap(lid[i], lid[0]);
1648  }
1649  }
1650 
1651  if (lid[0] == 4 || lid[0] == 5)
1652  {
1653  m_orientation = 0;
1654  }
1655  else if (lid[0] == 1 || lid[0] == 2)
1656  {
1657  // Rotate prism clockwise in p-r plane
1658  vector<NodeSharedPtr> vertexmap(6);
1659  vertexmap[0] = m_vertex[4];
1660  vertexmap[1] = m_vertex[0];
1661  vertexmap[2] = m_vertex[3];
1662  vertexmap[3] = m_vertex[5];
1663  vertexmap[4] = m_vertex[1];
1664  vertexmap[5] = m_vertex[2];
1665  m_vertex = vertexmap;
1666  m_orientation = 1;
1667  }
1668  else if (lid[0] == 0 || lid[0] == 3)
1669  {
1670  // Rotate prism counter-clockwise in p-r plane
1671  vector<NodeSharedPtr> vertexmap(6);
1672  vertexmap[0] = m_vertex[1];
1673  vertexmap[1] = m_vertex[4];
1674  vertexmap[2] = m_vertex[5];
1675  vertexmap[3] = m_vertex[2];
1676  vertexmap[4] = m_vertex[0];
1677  vertexmap[5] = m_vertex[3];
1678  m_vertex = vertexmap;
1679  m_orientation = 2;
1680  }
1681  else
1682  {
1683  cerr << "Warning: possible prism orientation problem." << endl;
1684  }
1685  }
1686 
1687 
1688  LibUtilities::ShapeType Hexahedron::m_type = GetElementFactory().
1689  RegisterCreatorFunction(LibUtilities::eHexahedron, Hexahedron::create, "Hexahedron");
1690 
1691  /**
1692  * @brief Create a hexahedral element.
1693  */
1694  Hexahedron::Hexahedron(ElmtConfig pConf,
1695  vector<NodeSharedPtr> pNodeList,
1696  vector<int> pTagList)
1697  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
1698  {
1699  m_tag = "H";
1700  m_dim = 3;
1701  m_taglist = pTagList;
1702  int n = m_conf.m_order-1;
1703 
1704  // Create a map to relate edge nodes to a pair of vertices defining an edge
1705  // This is based on the ordering produced by gmsh.
1706  map<pair<int,int>, int> edgeNodeMap;
1707  map<pair<int,int>, int>::iterator it;
1708  edgeNodeMap[pair<int,int>(1,2)] = 9;
1709  edgeNodeMap[pair<int,int>(1,4)] = 9 + n;
1710  edgeNodeMap[pair<int,int>(1,5)] = 9 + 2*n;
1711  edgeNodeMap[pair<int,int>(2,3)] = 9 + 3*n;
1712  edgeNodeMap[pair<int,int>(2,6)] = 9 + 4*n;
1713  edgeNodeMap[pair<int,int>(3,4)] = 9 + 5*n;
1714  edgeNodeMap[pair<int,int>(3,7)] = 9 + 6*n;
1715  edgeNodeMap[pair<int,int>(4,8)] = 9 + 7*n;
1716  edgeNodeMap[pair<int,int>(5,6)] = 9 + 8*n;
1717  edgeNodeMap[pair<int,int>(5,8)] = 9 + 9*n;
1718  edgeNodeMap[pair<int,int>(6,7)] = 9 + 10*n;
1719  edgeNodeMap[pair<int,int>(7,8)] = 9 + 11*n;
1720 
1721  // Add vertices
1722  for (int i = 0; i < 8; ++i) {
1723  m_vertex.push_back(pNodeList[i]);
1724  }
1725 
1726  // Create edges (with corresponding set of edge points)
1727  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
1728  {
1729  vector<NodeSharedPtr> edgeNodes;
1730  if (m_conf.m_order > 1) {
1731  for (int j = it->second; j < it->second + n; ++j) {
1732  edgeNodes.push_back(pNodeList[j-1]);
1733  }
1734  }
1735  m_edge.push_back(EdgeSharedPtr(new Edge(pNodeList[it->first.first-1],
1736  pNodeList[it->first.second-1],
1737  edgeNodes,
1738  m_conf.m_edgeCurveType)));
1739  }
1740 
1741  // Create faces
1742  int face_ids[6][4] = {
1743  {0,1,2,3},{0,1,5,4},{1,2,6,5},{2,3,7,6},{3,0,4,7},{4,5,6,7}};
1744  for (int j = 0; j < 6; ++j)
1745  {
1746  vector<NodeSharedPtr> faceVertices;
1747  vector<EdgeSharedPtr> faceEdges;
1748  vector<NodeSharedPtr> faceNodes;
1749  for (int k = 0; k < 4; ++k)
1750  {
1751  faceVertices.push_back(m_vertex[face_ids[j][k]]);
1752  NodeSharedPtr a = m_vertex[face_ids[j][k]];
1753  NodeSharedPtr b = m_vertex[face_ids[j][(k+1)%4]];
1754  for (unsigned int i = 0; i < m_edge.size(); ++i)
1755  {
1756  if ( ((*(m_edge[i]->m_n1)==*a) && (*(m_edge[i]->m_n2)==*b))
1757  || ((*(m_edge[i]->m_n1)==*b) && (*(m_edge[i]->m_n2) == *a)) )
1758  {
1759  faceEdges.push_back(m_edge[i]);
1760  break;
1761  }
1762  }
1763  }
1764 
1765  if (m_conf.m_faceNodes)
1766  {
1767  int N = 8 + 12*n + j*n*n;
1768  for (int i = 0; i < n*n; ++i)
1769  {
1770  faceNodes.push_back(pNodeList[N+i]);
1771  }
1772  }
1773  m_face.push_back(FaceSharedPtr(
1774  new Face(faceVertices, faceNodes, faceEdges, m_conf.m_faceCurveType)));
1775  }
1776  }
1777 
1778  SpatialDomains::GeometrySharedPtr Hexahedron::GetGeom(int coordDim)
1779  {
1780  SpatialDomains::QuadGeomSharedPtr faces[6];
1781  SpatialDomains::HexGeomSharedPtr ret;
1782 
1783  for (int i = 0; i < 6; ++i)
1784  {
1785  faces[i] = boost::dynamic_pointer_cast
1786  <SpatialDomains::QuadGeom>(m_face[i]->GetGeom(coordDim));
1787  }
1788 
1789  ret = MemoryManager<SpatialDomains::HexGeom>::
1790  AllocateSharedPtr(faces);
1791 
1792  return ret;
1793  }
1794 
1795  /**
1796  * @brief Return the number of nodes defining a hexahedron.
1797  */
1798  unsigned int Hexahedron::GetNumNodes(ElmtConfig pConf)
1799  {
1800  int n = pConf.m_order;
1801  if (pConf.m_faceNodes && pConf.m_volumeNodes)
1802  return (n+1)*(n+1)*(n+1);
1803  else if (pConf.m_faceNodes && !pConf.m_volumeNodes)
1804  return 6*(n+1)*(n+1)-12*(n+1)+8;
1805  else
1806  return 12*(n+1)-16;
1807  }
1808  }
1809 }
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2