Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Tetrahedron.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 
37 #include <LocalRegions/TetExp.h>
38 #include <SpatialDomains/TetGeom.h>
39 
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47 namespace NekMeshUtils
48 {
49 
50 LibUtilities::ShapeType Tetrahedron::m_type =
52  LibUtilities::eTetrahedron, Tetrahedron::create, "Tetrahedron");
53 
54 /**
55  * @brief Create a tetrahedron element.
56  */
57 Tetrahedron::Tetrahedron(ElmtConfig pConf,
58  vector<NodeSharedPtr> pNodeList,
59  vector<int> pTagList)
60  : Element(pConf, GetNumNodes(pConf), pNodeList.size())
61 {
62  m_tag = "A";
63  m_dim = 3;
64  m_taglist = pTagList;
65  int n = m_conf.m_order - 1;
66 
67  // Create a map to relate edge nodes to a pair of vertices
68  // defining an edge.
69  map<pair<int, int>, int> edgeNodeMap;
70  map<pair<int, int>, int>::iterator it;
71  edgeNodeMap[pair<int, int>(1, 2)] = 5;
72  edgeNodeMap[pair<int, int>(2, 3)] = 5 + n;
73  edgeNodeMap[pair<int, int>(1, 3)] = 5 + 2 * n;
74  edgeNodeMap[pair<int, int>(1, 4)] = 5 + 3 * n;
75  edgeNodeMap[pair<int, int>(2, 4)] = 5 + 4 * n;
76  edgeNodeMap[pair<int, int>(3, 4)] = 5 + 5 * n;
77 
78  // Add vertices
79  for (int i = 0; i < 4; ++i)
80  {
81  m_vertex.push_back(pNodeList[i]);
82  }
83 
84  // Create edges (with corresponding set of edge points)
85  int eid = 0;
86  for (it = edgeNodeMap.begin(); it != edgeNodeMap.end(); ++it)
87  {
88  vector<NodeSharedPtr> edgeNodes;
89  if (m_conf.m_order > 1)
90  {
91  for (int j = it->second; j < it->second + n; ++j)
92  {
93  edgeNodes.push_back(pNodeList[j - 1]);
94  }
95  }
96  m_edge.push_back(EdgeSharedPtr(new Edge(pNodeList[it->first.first - 1],
97  pNodeList[it->first.second - 1],
98  edgeNodes,
100  m_edge.back()->m_id = eid++;
101  }
102 
103  // Reorient the tet to ensure collapsed coordinates align between
104  // adjacent elements.
105  if (m_conf.m_reorient)
106  {
107  OrientTet();
108  }
109  else
110  {
111  // If we didn't need to orient the tet then set up the
112  // orientation map as the identity mapping.
113  for (int i = 0; i < 4; ++i)
114  {
115  m_orientationMap[i] = i;
116  }
117  }
118 
119  // Face-vertex IDs
120  int face_ids[4][3] = {{0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
121 
122  // Face-edge IDs
123  int face_edges[4][3];
124 
125  // Create faces
126  for (int j = 0; j < 4; ++j)
127  {
128  vector<NodeSharedPtr> faceVertices;
129  vector<EdgeSharedPtr> faceEdges;
130  vector<NodeSharedPtr> faceNodes;
131 
132  // Extract the edges for this face. We need to do this because
133  // of the reorientation which might have been applied (see the
134  // additional note below).
135  for (int k = 0; k < 3; ++k)
136  {
137  faceVertices.push_back(m_vertex[face_ids[j][k]]);
138  NodeSharedPtr a = m_vertex[face_ids[j][k]];
139  NodeSharedPtr b = m_vertex[face_ids[j][(k + 1) % 3]];
140  for (unsigned int i = 0; i < m_edge.size(); ++i)
141  {
142  if (((*(m_edge[i]->m_n1) == *a) &&
143  (*(m_edge[i]->m_n2) == *b)) ||
144  ((*(m_edge[i]->m_n1) == *b) && (*(m_edge[i]->m_n2) == *a)))
145  {
146  face_edges[j][k] = i;
147  faceEdges.push_back(m_edge[i]);
148  break;
149  }
150  }
151  }
152 
153  // When face curvature is supplied, it may have been the case
154  // that our tetrahedron was reoriented. In this case, faces have
155  // different vertex IDs and so we have to rotate the face
156  // curvature so that the two align appropriately.
157  if (m_conf.m_faceNodes)
158  {
159  const int nFaceNodes = n * (n - 1) / 2;
160 
161  // Get the vertex IDs of whatever face we are processing.
162  vector<int> faceIds(3);
163  faceIds[0] = faceVertices[0]->m_id;
164  faceIds[1] = faceVertices[1]->m_id;
165  faceIds[2] = faceVertices[2]->m_id;
166 
167  // Find out the original face number as we were given it in
168  // the constructor using the orientation map.
169  int origFace = -1;
170  for (int i = 0; i < 4; ++i)
171  {
172  if (m_orientationMap[i] == j)
173  {
174  origFace = i;
175  break;
176  }
177  }
178 
179  ASSERTL0(origFace >= 0, "Couldn't find face");
180 
181  // Now get the face nodes for the original face.
182  int N = 4 + 6 * n + origFace * nFaceNodes;
183  for (int i = 0; i < nFaceNodes; ++i)
184  {
185  faceNodes.push_back(pNodeList[N + i]);
186  }
187 
188  // Find the original face vertex IDs.
189  vector<int> origFaceIds(3);
190  origFaceIds[0] = pNodeList[face_ids[origFace][0]]->m_id;
191  origFaceIds[1] = pNodeList[face_ids[origFace][1]]->m_id;
192  origFaceIds[2] = pNodeList[face_ids[origFace][2]]->m_id;
193 
194  // Construct a HOTriangle object which performs the
195  // orientation magically for us.
196  HOTriangle<NodeSharedPtr> hoTri(origFaceIds, faceNodes);
197  hoTri.Align(faceIds);
198 
199  // Copy the face nodes back again.
200  faceNodes = hoTri.surfVerts;
201  }
202 
203  m_face.push_back(FaceSharedPtr(new Face(
204  faceVertices, faceNodes, faceEdges, m_conf.m_faceCurveType)));
205  }
206 
207  vector<EdgeSharedPtr> tmp(6);
208  tmp[0] = m_edge[face_edges[0][0]];
209  tmp[1] = m_edge[face_edges[0][1]];
210  tmp[2] = m_edge[face_edges[0][2]];
211  tmp[3] = m_edge[face_edges[1][2]];
212  tmp[4] = m_edge[face_edges[1][1]];
213  tmp[5] = m_edge[face_edges[2][1]];
214  m_edge = tmp;
215 }
216 
218 {
221 
222  for (int i = 0; i < 4; ++i)
223  {
224  tfaces[i] = boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
225  m_face[i]->GetGeom(coordDim));
226  }
227 
229 
230  return ret;
231 }
232 
233 /**
234  * @brief Return the number of nodes defining a tetrahedron.
235  */
237 {
238  int n = pConf.m_order;
239  if (pConf.m_volumeNodes && pConf.m_faceNodes)
240  return (n + 1) * (n + 2) * (n + 3) / 6;
241  else if (!pConf.m_volumeNodes && pConf.m_faceNodes)
242  return 4 * (n + 1) * (n + 2) / 2 - 6 * (n + 1) + 4;
243  else
244  return 6 * (n + 1) - 8;
245 }
246 
247 /**
248  * @brief .
249  */
250 void Tetrahedron::Complete(int order)
251 {
252  int i, j;
253 
254  // Create basis key for a nodal tetrahedron.
257  order + 1,
258  LibUtilities::PointsKey(order + 1,
262  order + 1,
263  LibUtilities::PointsKey(order + 1,
267  order + 1,
268  LibUtilities::PointsKey(order + 1,
270 
271  // Create a standard nodal tetrahedron in order to get the
272  // Vandermonde matrix to perform interpolation to nodal points.
276 
277  Array<OneD, NekDouble> x, y, z;
278 
279  nodalTet->GetNodalPoints(x, y, z);
280 
282  boost::dynamic_pointer_cast<SpatialDomains::TetGeom>(this->GetGeom(3));
283 
284  // Create basis key for a tetrahedron.
287  order + 1,
288  LibUtilities::PointsKey(order + 1,
292  order + 1,
293  LibUtilities::PointsKey(order + 1,
297  order + 1,
298  LibUtilities::PointsKey(order + 1,
300 
301  // Create a tet.
304  C0, C1, C2, geom);
305 
306  // Get coordinate array for tetrahedron.
307  int nqtot = tet->GetTotPoints();
308  Array<OneD, NekDouble> alloc(6 * nqtot);
309  Array<OneD, NekDouble> xi(alloc);
310  Array<OneD, NekDouble> yi(alloc + nqtot);
311  Array<OneD, NekDouble> zi(alloc + 2 * nqtot);
312  Array<OneD, NekDouble> xo(alloc + 3 * nqtot);
313  Array<OneD, NekDouble> yo(alloc + 4 * nqtot);
314  Array<OneD, NekDouble> zo(alloc + 5 * nqtot);
316 
317  tet->GetCoords(xi, yi, zi);
318 
319  for (i = 0; i < 3; ++i)
320  {
321  Array<OneD, NekDouble> coeffs(nodalTet->GetNcoeffs());
322  tet->FwdTrans(alloc + i * nqtot, coeffs);
323  // Apply Vandermonde matrix to project onto nodal space.
324  nodalTet->ModalToNodal(coeffs, tmp = alloc + (i + 3) * nqtot);
325  }
326 
327  // Now extract points from the co-ordinate arrays into the
328  // edge/face/volume nodes. First, extract edge-interior nodes.
329  for (i = 0; i < 6; ++i)
330  {
331  int pos = 4 + i * (order - 1);
332  m_edge[i]->m_edgeNodes.clear();
333  for (j = 0; j < order - 1; ++j)
334  {
335  m_edge[i]->m_edgeNodes.push_back(NodeSharedPtr(
336  new Node(0, xo[pos + j], yo[pos + j], zo[pos + j])));
337  }
338  }
339 
340  // Now extract face-interior nodes.
341  for (i = 0; i < 4; ++i)
342  {
343  int pos = 4 + 6 * (order - 1) + i * (order - 2) * (order - 1) / 2;
344  m_face[i]->m_faceNodes.clear();
345  for (j = 0; j < (order - 2) * (order - 1) / 2; ++j)
346  {
347  m_face[i]->m_faceNodes.push_back(NodeSharedPtr(
348  new Node(0, xo[pos + j], yo[pos + j], zo[pos + j])));
349  }
350  }
351 
352  // Finally extract volume nodes.
353  int pos = 4 + 6 * (order - 1) + 4 * (order - 2) * (order - 1) / 2;
354  for (i = pos; i < (order + 1) * (order + 2) * (order + 3) / 6; ++i)
355  {
356  m_volumeNodes.push_back(
357  NodeSharedPtr(new Node(0, xo[i], yo[i], zo[i])));
358  }
359 
360  m_conf.m_order = order;
361  m_conf.m_faceNodes = true;
362  m_conf.m_volumeNodes = true;
363 }
364 
365 struct TetOrient
366 {
367  TetOrient(vector<int> nid, int fid) : nid(nid), fid(fid)
368  {
369  }
370  vector<int> nid;
371  int fid;
372 };
373 
374 struct TetOrientHash : std::unary_function<struct TetOrient, std::size_t>
375 {
376  std::size_t operator()(struct TetOrient const &p) const
377  {
378  return boost::hash_range(p.nid.begin(), p.nid.end());
379  }
380 };
381 typedef boost::unordered_set<struct TetOrient, TetOrientHash> TetOrientSet;
382 
383 bool operator==(const struct TetOrient &a, const struct TetOrient &b)
384 {
385  if (a.nid.size() != b.nid.size())
386  {
387  return false;
388  }
389 
390  for (int i = 0; i < a.nid.size(); ++i)
391  {
392  if (a.nid[i] != b.nid[i])
393  {
394  return false;
395  }
396  }
397 
398  return true;
399 }
400 
401 /**
402  * @brief Orient tetrahedron to align degenerate vertices.
403  *
404  * Orientation of tetrahedral elements is required so that the
405  * singular vertices of triangular faces (which occur as a part of the
406  * collapsed co-ordinate system) align. The algorithm is based on that
407  * used in T. Warburton's thesis and in the original Nektar source.
408  *
409  * First the vertices are ordered with the highest global vertex at
410  * the top degenerate point, and the base degenerate point has second
411  * lowest ID. These vertices are swapped if the element is incorrectly
412  * oriented.
413  */
415 {
416  static int face_ids[4][3] = {{0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
417 
418  TetOrientSet faces;
419 
420  // Create a copy of the original vertex ordering. This is used to
421  // construct a mapping, #orientationMap, which maps the original
422  // face ordering to the new face ordering.
423  for (int i = 0; i < 4; ++i)
424  {
425  vector<int> nodes(3);
426 
427  nodes[0] = m_vertex[face_ids[i][0]]->m_id;
428  nodes[1] = m_vertex[face_ids[i][1]]->m_id;
429  nodes[2] = m_vertex[face_ids[i][2]]->m_id;
430 
431  sort(nodes.begin(), nodes.end());
432  struct TetOrient faceNodes(nodes, i);
433  faces.insert(faceNodes);
434  }
435 
436  // Store a copy of the original vertex ordering so we can create a
437  // permutation map later.
438  vector<NodeSharedPtr> origVert = m_vertex;
439 
440  // Order vertices with highest global vertex at top degenerate
441  // point. Place second highest global vertex at base degenerate
442  // point.
443  sort(m_vertex.begin(), m_vertex.end());
444 
445  // Calculate a.(b x c) if negative, reverse order of
446  // non-degenerate points to correctly orientate the tet.
447 
448  NekDouble ax = m_vertex[1]->m_x - m_vertex[0]->m_x;
449  NekDouble ay = m_vertex[1]->m_y - m_vertex[0]->m_y;
450  NekDouble az = m_vertex[1]->m_z - m_vertex[0]->m_z;
451  NekDouble bx = m_vertex[2]->m_x - m_vertex[0]->m_x;
452  NekDouble by = m_vertex[2]->m_y - m_vertex[0]->m_y;
453  NekDouble bz = m_vertex[2]->m_z - m_vertex[0]->m_z;
454  NekDouble cx = m_vertex[3]->m_x - m_vertex[0]->m_x;
455  NekDouble cy = m_vertex[3]->m_y - m_vertex[0]->m_y;
456  NekDouble cz = m_vertex[3]->m_z - m_vertex[0]->m_z;
457 
458  NekDouble nx = (ay * bz - az * by);
459  NekDouble ny = (az * bx - ax * bz);
460  NekDouble nz = (ax * by - ay * bx);
461  NekDouble nmag = sqrt(nx * nx + ny * ny + nz * nz);
462  nx /= nmag;
463  ny /= nmag;
464  nz /= nmag;
465 
466  // distance of top vertex from base
467  NekDouble dist = cx * nx + cy * ny + cz * nz;
468 
469  if (fabs(dist) <= 1e-4)
470  {
471  cerr << "Warning: degenerate tetrahedron, 3rd vertex is = " << dist
472  << " from face" << endl;
473  }
474 
475  if (dist < 0)
476  {
477  swap(m_vertex[0], m_vertex[1]);
478  }
479 
480  nx = (ay * cz - az * cy);
481  ny = (az * cx - ax * cz);
482  nz = (ax * cy - ay * cx);
483  nmag = sqrt(nx * nx + ny * ny + nz * nz);
484  nx /= nmag;
485  ny /= nmag;
486  nz /= nmag;
487 
488  // distance of top vertex from base
489  dist = bx * nx + by * ny + bz * nz;
490 
491  if (fabs(dist) <= 1e-4)
492  {
493  cerr << "Warning: degenerate tetrahedron, 2nd vertex is = " << dist
494  << " from face" << endl;
495  }
496 
497  nx = (by * cz - bz * cy);
498  ny = (bz * cx - bx * cz);
499  nz = (bx * cy - by * cx);
500  nmag = sqrt(nx * nx + ny * ny + nz * nz);
501  nx /= nmag;
502  ny /= nmag;
503  nz /= nmag;
504 
505  // distance of top vertex from base
506  dist = ax * nx + ay * ny + az * nz;
507 
508  if (fabs(dist) <= 1e-4)
509  {
510  cerr << "Warning: degenerate tetrahedron, 1st vertex is = " << dist
511  << " from face" << endl;
512  }
513 
515 
516  // Search for the face in the original set of face nodes. Then use
517  // this to construct the #orientationMap.
518  for (int i = 0; i < 4; ++i)
519  {
520  vector<int> nodes(3);
521 
522  nodes[0] = m_vertex[face_ids[i][0]]->m_id;
523  nodes[1] = m_vertex[face_ids[i][1]]->m_id;
524  nodes[2] = m_vertex[face_ids[i][2]]->m_id;
525  sort(nodes.begin(), nodes.end());
526 
527  struct TetOrient faceNodes(nodes, 0);
528 
529  it = faces.find(faceNodes);
530  m_orientationMap[it->fid] = i;
531 
532  for (int j = 0; j < 4; ++j)
533  {
534  if (m_vertex[i]->m_id == origVert[j]->m_id)
535  {
536  m_origVertMap[j] = i;
537  break;
538  }
539  }
540  }
541 }
542 }
543 }
bool m_faceNodes
Denotes whether the element contains face nodes. For 2D elements, if this is true then the element co...
Definition: Element.h:89
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Basic information about an element.
Definition: Element.h:58
LibUtilities::PointsType m_faceCurveType
Distribution of points in faces.
Definition: Element.h:103
Represents an edge which joins two points.
Definition: Edge.h:61
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Align(std::vector< int > vertId)
Align this surface to a given vertex ID.
Definition: Triangle.h:138
Represents a face comprised of three or more edges.
Definition: Face.h:60
virtual NEKMESHUTILS_EXPORT SpatialDomains::GeometrySharedPtr GetGeom(int coordDim)
Generate a Nektar++ geometry object for this element.
STL namespace.
std::vector< T > surfVerts
The triangle surface vertices – templated so that this can either be nodes or IDs.
Definition: Triangle.h:68
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
ElmtConfig m_conf
Contains configuration of the element.
Definition: Element.h:496
boost::shared_ptr< TetExp > TetExpSharedPtr
Definition: TetExp.h:223
Represents a point in the domain.
Definition: Node.h:60
std::vector< int > m_taglist
List of integers specifying properties of the element.
Definition: Element.h:500
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
bool operator==(ElmtConfig const &c1, ElmtConfig const &c2)
std::size_t operator()(struct TetOrient const &p) const
LibUtilities::PointsType m_edgeCurveType
Distribution of points in edges.
Definition: Element.h:101
unsigned int m_order
Order of the element.
Definition: Element.h:96
Principle Orthogonal Functions .
Definition: BasisType.h:47
boost::unordered_set< struct TetOrient, TetOrientHash > TetOrientSet
std::vector< NodeSharedPtr > m_vertex
List of element vertex nodes.
Definition: Element.h:502
unsigned int m_dim
Dimension of the element.
Definition: Element.h:494
bool m_volumeNodes
Denotes whether the element contains volume (i.e. interior) nodes. These are not supported by either ...
Definition: Element.h:94
std::vector< EdgeSharedPtr > m_edge
List of element edges.
Definition: Element.h:504
Principle Orthogonal Functions .
Definition: BasisType.h:48
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
Principle Orthogonal Functions .
Definition: BasisType.h:46
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
std::vector< NodeSharedPtr > m_volumeNodes
List of element volume nodes.
Definition: Element.h:508
boost::shared_ptr< StdNodalTetExp > StdNodalTetExpSharedPtr
void OrientTet()
Orient tetrahedron to align degenerate vertices.
3D Evenly-spaced points on a Tetrahedron
Definition: PointsType.h:71
A lightweight struct for dealing with high-order triangle alignment.
Definition: Triangle.h:53
std::string m_tag
Tag character describing the element.
Definition: Element.h:498
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:196
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
static NEKMESHUTILS_EXPORT unsigned int GetNumNodes(ElmtConfig pConf)
Return the number of nodes defining a tetrahedron.
unsigned int m_id
ID of the element.
Definition: Element.h:492
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
boost::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:106
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
std::vector< FaceSharedPtr > m_face
List of element faces.
Definition: Element.h:506
bool m_reorient
Denotes whether the element needs to be re-orientated for a spectral element framework.
Definition: Element.h:99
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
Definition: Face.h:378
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
virtual NEKMESHUTILS_EXPORT void Complete(int order)
Describes the specification for a Basis.
Definition: Basis.h:50
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
TetOrient(vector< int > nid, int fid)
Base class for element definitions.
Definition: Element.h:115
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215