Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessTetSplit.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessTetSplit.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: Split prisms -> tets
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <boost/unordered_map.hpp>
37 
41 #include <LocalRegions/PrismExp.h>
42 
44 
45 #include "ProcessTetSplit.h"
46 
47 using namespace std;
48 using namespace Nektar::NekMeshUtils;
49 
50 namespace Nektar
51 {
52 namespace Utilities
53 {
54 
55 typedef std::pair<int, int> ipair;
56 
57 ModuleKey ProcessTetSplit::className =
59  ModuleKey(eProcessModule, "tetsplit"),
60  ProcessTetSplit::create,
61  "Split prismatic elements to tetrahedra");
62 
63 ProcessTetSplit::ProcessTetSplit(MeshSharedPtr m) : ProcessModule(m)
64 {
65  m_config["nq"] =
66  ConfigOption(false, "3", "Number of points in high order elements.");
67 }
68 
70 {
71 }
72 
74 {
75  int nodeId = m_mesh->m_vertexSet.size();
76 
77  // Set up map which identifies edges (as pairs of vertex ids)
78  // including their vertices to the offset/stride in the 3d array of
79  // collapsed co-ordinates. Note that this map also includes the
80  // diagonal edges along quadrilateral faces which will be used to
81  // add high-order information to the split tetrahedra.
82  map<ipair, ipair> edgeMap;
84 
85  int nq = m_config["nq"].as<int>();
86  int ne = nq - 2; // Edge interior
87  int nft = (ne - 1) * ne / 2; // Face interior (triangle)
88  int nfq = ne * ne; // Face interior (quad)
89  int o = 6;
90 
91  // Standard prismatic edges (0->8)
92  edgeMap[ipair(0, 1)] = ipair(o, 1);
93  edgeMap[ipair(1, 2)] = ipair(o + 1 * ne, 1);
94  edgeMap[ipair(3, 2)] = ipair(o + 2 * ne, 1);
95  edgeMap[ipair(0, 3)] = ipair(o + 3 * ne, 1);
96  edgeMap[ipair(0, 4)] = ipair(o + 4 * ne, 1);
97  edgeMap[ipair(1, 4)] = ipair(o + 5 * ne, 1);
98  edgeMap[ipair(2, 5)] = ipair(o + 6 * ne, 1);
99  edgeMap[ipair(3, 5)] = ipair(o + 7 * ne, 1);
100  edgeMap[ipair(4, 5)] = ipair(o + 8 * ne, 1);
101 
102  // Face 0 diagonals
103  o += 9 * ne;
104  edgeMap[ipair(0, 2)] = ipair(o, ne + 1);
105  edgeMap[ipair(1, 3)] = ipair(o + ne - 1, ne - 1);
106 
107  // Face 2 diagonals
108  o += nfq + nft;
109  edgeMap[ipair(1, 5)] = ipair(o, ne + 1);
110  edgeMap[ipair(2, 4)] = ipair(o + ne - 1, ne - 1);
111 
112  // Face 4 diagonals
113  o += nfq + nft;
114  edgeMap[ipair(0, 5)] = ipair(o, ne + 1);
115  edgeMap[ipair(3, 4)] = ipair(o + ne - 1, ne - 1);
116 
117  // Default PointsType.
119 
120  /*
121  * Split all element types into tetrahedra. This is based on the
122  * algorithm found in:
123  *
124  * "How to Subdivide Pyramids, Prisms and Hexahedra into
125  * Tetrahedra", J. Dompierre et al.
126  */
127 
128  // Denotes a set of indices inside m_mesh->m_element[m_mesh->m_expDim-1]
129  // which are
130  // to be removed. These are precisely the quadrilateral boundary
131  // faces which will be replaced by two triangular faces.
132  set<int> toRemove;
133 
134  // Represents table 2 of paper; each row i represents a rotation of
135  // the prism nodes such that vertex i is placed at position 0.
136  static int indir[6][6] = {{0, 1, 2, 3, 4, 5},
137  {1, 2, 0, 4, 5, 3},
138  {2, 0, 1, 5, 3, 4},
139  {3, 5, 4, 0, 2, 1},
140  {4, 3, 5, 1, 0, 2},
141  {5, 4, 3, 2, 1, 0}};
142 
143  // Represents table 3 of paper; the first three rows are the three
144  // tetrahedra if the first condition is met; the latter three rows
145  // are the three tetrahedra if the second condition is met.
146  static int prismTet[6][4] = {{0, 1, 2, 5},
147  {0, 1, 5, 4},
148  {0, 4, 5, 3},
149  {0, 1, 2, 4},
150  {0, 4, 2, 5},
151  {0, 4, 5, 3}};
152 
153  // Represents the order of tetrahedral edges (in Nektar++ ordering).
154  static int tetEdges[6][2] = {
155  {0, 1}, {1, 2}, {0, 2}, {0, 3}, {1, 3}, {2, 3}};
156 
157  // A tetrahedron nodes -> faces map.
158  static int tetFaceNodes[4][3] = {
159  {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
160 
161  static NodeSharedPtr stdPrismNodes[6] = {
162  NodeSharedPtr(new Node(0, -1, -1, -1)),
163  NodeSharedPtr(new Node(1, 1, -1, -1)),
164  NodeSharedPtr(new Node(2, 1, 1, -1)),
165  NodeSharedPtr(new Node(3, -1, 1, -1)),
166  NodeSharedPtr(new Node(4, -1, -1, 1)),
167  NodeSharedPtr(new Node(5, -1, 1, 1))};
168 
169  // Generate equally spaced nodal points on triangle.
170  Array<OneD, NekDouble> rp(nq * (nq + 1) / 2), sp(nq * (nq + 1) / 2);
172  LibUtilities::PointsManager()[elec]->GetPoints(rp, sp);
173 
174  // Make a copy of the element list.
175  vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
176  m_mesh->m_element[m_mesh->m_expDim].clear();
177 
178  for (int i = 0; i < el.size(); ++i)
179  {
180  if (el[i]->GetConf().m_e != LibUtilities::ePrism)
181  {
182  m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
183  continue;
184  }
185 
186  vector<NodeSharedPtr> nodeList(6);
187 
188  // Map Nektar++ ordering (vertices 0,1,2,3 are base quad) to
189  // paper ordering (vertices 0,1,2 are first triangular face).
190  int mapPrism[6] = {0, 1, 4, 3, 2, 5};
191  for (int j = 0; j < 6; ++j)
192  {
193  nodeList[j] = el[i]->GetVertex(mapPrism[j]);
194  }
195 
196  // Determine minimum ID of the nodes in this prism.
197  int minElId = nodeList[0]->m_id;
198  int minId = 0;
199  for (int j = 1; j < 6; ++j)
200  {
201  int curId = nodeList[j]->m_id;
202  if (curId < minElId)
203  {
204  minElId = curId;
205  minId = j;
206  }
207  }
208 
209  int offset;
210 
211  // Split prism using paper criterion.
212  int id1 = min(nodeList[indir[minId][1]]->m_id,
213  nodeList[indir[minId][5]]->m_id);
214  int id2 = min(nodeList[indir[minId][2]]->m_id,
215  nodeList[indir[minId][4]]->m_id);
216 
217  if (id1 < id2)
218  {
219  offset = 0;
220  }
221  else if (id1 > id2)
222  {
223  offset = 3;
224  }
225  else
226  {
227  cerr << "Connectivity issue with prism->tet splitting." << endl;
228  abort();
229  }
230 
231  // Create local prismatic region so that co-ordinates of the
232  // mapped element can be read from.
234  boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(
235  el[i]->GetGeom(m_mesh->m_spaceDim));
238  nq,
242  nq,
246  B0, B0, B1, geomLayer);
247 
248  // Get the coordinates of the high order prismatic element.
249  Array<OneD, NekDouble> wsp(3 * nq * nq * nq);
251  x[0] = wsp;
252  x[1] = wsp + 1 * nq * nq * nq;
253  x[2] = wsp + 2 * nq * nq * nq;
254  qs->GetCoords(x[0], x[1], x[2]);
255 
260 
261  // Process face data. Initially put coordinates into equally
262  // spaced nodal distribution.
266 
267  int nCoeffs = nodalPrism->GetNcoeffs();
268  Array<OneD, NekDouble> wsp2(3 * nCoeffs);
270  xn[0] = wsp2;
271  xn[1] = wsp2 + 1 * nCoeffs;
272  xn[2] = wsp2 + 2 * nCoeffs;
273 
274  for (int j = 0; j < 3; ++j)
275  {
276  Array<OneD, NekDouble> tmp(nodalPrism->GetNcoeffs());
277  qs->FwdTrans(x[j], tmp);
278  nodalPrism->ModalToNodal(tmp, xn[j]);
279  }
280 
281  // Store map of nodes.
282  map<int, int> prismVerts;
283  for (int j = 0; j < 6; ++j)
284  {
285  prismVerts[el[i]->GetVertex(j)->m_id] = j;
286  }
287 
288  for (int j = 0; j < 3; ++j)
289  {
290  vector<NodeSharedPtr> tetNodes(4);
291 
292  // Extract vertices for tetrahedron.
293  for (int k = 0; k < 4; ++k)
294  {
295  tetNodes[k] = nodeList[indir[minId][prismTet[j + offset][k]]];
296  }
297 
298  // Add high order information to tetrahedral edges.
299  for (int k = 0; k < 6; ++k)
300  {
301  // Determine prismatic nodes which correspond with this
302  // edge. Apply prism map to this to get Nektar++
303  // ordering (this works since as a permutation,
304  // prismMap^2 = id).
305  int n1 = mapPrism
306  [indir[minId][prismTet[j + offset][tetEdges[k][0]]]];
307  int n2 = mapPrism
308  [indir[minId][prismTet[j + offset][tetEdges[k][1]]]];
309 
310  // Find offset/stride
311  it = edgeMap.find(pair<int, int>(n1, n2));
312  if (it == edgeMap.end())
313  {
314  it = edgeMap.find(pair<int, int>(n2, n1));
315  if (it == edgeMap.end())
316  {
317  cerr << "Couldn't find prism edges " << n1 << " " << n2
318  << endl;
319  abort();
320  }
321  // Extract vertices -- reverse order.
322  for (int l = ne - 1; l >= 0; --l)
323  {
324  int pos = it->second.first + l * it->second.second;
325  tetNodes.push_back(NodeSharedPtr(new Node(
326  nodeId++, xn[0][pos], xn[1][pos], xn[2][pos])));
327  }
328  }
329  else
330  {
331  // Extract vertices -- forwards order.
332  for (int l = 0; l < ne; ++l)
333  {
334  int pos = it->second.first + l * it->second.second;
335  tetNodes.push_back(NodeSharedPtr(new Node(
336  nodeId++, xn[0][pos], xn[1][pos], xn[2][pos])));
337  }
338  }
339  }
340 
341  // Create new tetrahedron with edge curvature.
342  vector<int> tags = el[i]->GetTagList();
343  ElmtConfig conf(LibUtilities::eTetrahedron, nq - 1, false, false);
345  LibUtilities::eTetrahedron, conf, tetNodes, tags);
346 
347  // Extract interior face data.
348  for (int k = 0; k < 4; ++k)
349  {
350  // First determine which nodes of prism are being used
351  // for this face.
352  FaceSharedPtr face = elmt->GetFace(k);
353  vector<NodeSharedPtr> triNodes(3);
354 
355  for (int l = 0; l < 3; ++l)
356  {
357  NodeSharedPtr v = face->m_vertexList[l];
358  triNodes[l] = stdPrismNodes[prismVerts[v->m_id]];
359  }
360 
361  // Create a triangle with the standard nodes of the
362  // prism.
363  vector<int> tags;
364  ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);
366  LibUtilities::eTriangle, conf, triNodes, tags);
367  SpatialDomains::GeometrySharedPtr triGeom = elmt->GetGeom(3);
368  triGeom->FillGeom();
369  o = 3 + 3 * ne;
370  face->m_curveType = LibUtilities::eNodalTriEvenlySpaced;
371 
372  for (int l = 0; l < nft; ++l)
373  {
374  Array<OneD, NekDouble> tmp1(2), tmp2(3);
375  tmp1[0] = rp[o + l];
376  tmp1[1] = sp[o + l];
377  tmp2[0] = triGeom->GetCoord(0, tmp1);
378  tmp2[1] = triGeom->GetCoord(1, tmp1);
379  tmp2[2] = triGeom->GetCoord(2, tmp1);
380 
381  // PhysEvaluate on prism geometry.
382  NekDouble xc = geomLayer->GetCoord(0, tmp2);
383  NekDouble yc = geomLayer->GetCoord(1, tmp2);
384  NekDouble zc = geomLayer->GetCoord(2, tmp2);
385  face->m_faceNodes.push_back(
386  NodeSharedPtr(new Node(nodeId++, xc, yc, zc)));
387  }
388  }
389 
390  m_mesh->m_element[m_mesh->m_expDim].push_back(elmt);
391  }
392 
393  // Now check to see if this one of the quadrilateral faces is
394  // associated with a boundary condition. If it is, we split the
395  // face into two triangles and mark the existing face for
396  // removal.
397  //
398  // Note that this algorithm has significant room for improvement
399  // and is likely one of the least optimal approachs - however
400  // implementation is simple.
401  for (int fid = 0; fid < 5; fid += 2)
402  {
403  int bl = el[i]->GetBoundaryLink(fid);
404 
405  if (bl == -1)
406  {
407  continue;
408  }
409 
410  vector<NodeSharedPtr> triNodeList(3);
411  vector<int> faceNodes(3);
412  vector<int> tmp;
413  vector<int> tagBE;
414  ElmtConfig bconf(LibUtilities::eTriangle, 1, true, true);
415  ElementSharedPtr elmt;
416 
417  // Mark existing boundary face for removal.
418  toRemove.insert(bl);
419  tagBE = m_mesh->m_element[m_mesh->m_expDim - 1][bl]->GetTagList();
420 
421  // First loop over tets.
422  for (int j = 0; j < 3; ++j)
423  {
424  // Now loop over faces.
425  for (int k = 0; k < 4; ++k)
426  {
427  // Finally determine Nektar++ local node numbers for
428  // this face.
429  for (int l = 0; l < 3; ++l)
430  {
431  faceNodes[l] = mapPrism
432  [indir[minId]
433  [prismTet[j + offset][tetFaceNodes[k][l]]]];
434  }
435 
436  tmp = faceNodes;
437  sort(faceNodes.begin(), faceNodes.end());
438 
439  // If this face matches a triple denoting a split
440  // quad face, add the face to the expansion list.
441  if ((fid == 0 && ((faceNodes[0] == 0 && faceNodes[1] == 1 &&
442  faceNodes[2] == 2) ||
443  (faceNodes[0] == 0 && faceNodes[1] == 2 &&
444  faceNodes[2] == 3) ||
445  (faceNodes[0] == 0 && faceNodes[1] == 1 &&
446  faceNodes[2] == 3) ||
447  (faceNodes[0] == 1 && faceNodes[1] == 2 &&
448  faceNodes[2] == 3))) ||
449  (fid == 2 && ((faceNodes[0] == 1 && faceNodes[1] == 2 &&
450  faceNodes[2] == 5) ||
451  (faceNodes[0] == 1 && faceNodes[1] == 4 &&
452  faceNodes[2] == 5) ||
453  (faceNodes[0] == 1 && faceNodes[1] == 2 &&
454  faceNodes[2] == 4) ||
455  (faceNodes[0] == 2 && faceNodes[1] == 4 &&
456  faceNodes[2] == 5))) ||
457  (fid == 4 && ((faceNodes[0] == 0 && faceNodes[1] == 3 &&
458  faceNodes[2] == 5) ||
459  (faceNodes[0] == 0 && faceNodes[1] == 4 &&
460  faceNodes[2] == 5) ||
461  (faceNodes[0] == 0 && faceNodes[1] == 3 &&
462  faceNodes[2] == 4) ||
463  (faceNodes[0] == 3 && faceNodes[1] == 4 &&
464  faceNodes[2] == 5))))
465  {
466  triNodeList[0] = nodeList[mapPrism[tmp[0]]];
467  triNodeList[1] = nodeList[mapPrism[tmp[1]]];
468  triNodeList[2] = nodeList[mapPrism[tmp[2]]];
470  LibUtilities::eTriangle, bconf, triNodeList, tagBE);
471  m_mesh->m_element[m_mesh->m_expDim - 1].push_back(elmt);
472  }
473  }
474  }
475  }
476  }
477 
478  // Remove 2D elements.
479  vector<ElementSharedPtr> tmp;
480  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim - 1].size(); ++i)
481  {
482  set<int>::iterator it = toRemove.find(i);
483  if (it == toRemove.end())
484  {
485  tmp.push_back(m_mesh->m_element[m_mesh->m_expDim - 1][i]);
486  }
487  }
488 
489  m_mesh->m_element[m_mesh->m_expDim - 1] = tmp;
490 
491  // Re-process mesh to eliminate duplicate vertices and edges.
492  ProcessVertices();
493  ProcessEdges();
494  ProcessFaces();
495  ProcessElements();
497 }
498 }
499 }
Basic information about an element.
Definition: Element.h:58
pair< ModuleType, string > ModuleKey
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
map< string, ConfigOption > m_config
List of configuration values.
Principle Modified Functions .
Definition: BasisType.h:49
STL namespace.
MeshSharedPtr m_mesh
Mesh object.
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
Represents a point in the domain.
Definition: Node.h:60
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
virtual void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
Principle Orthogonal Functions .
Definition: BasisType.h:47
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:63
virtual void ProcessVertices()
Extract element vertices.
boost::shared_ptr< StdNodalPrismExp > StdNodalPrismExpSharedPtr
std::pair< int, int > ipair
virtual void ProcessElements()
Generate element IDs.
Principle Modified Functions .
Definition: BasisType.h:50
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
PointsManagerT & PointsManager(void)
Principle Orthogonal Functions .
Definition: BasisType.h:46
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
boost::shared_ptr< PrismExp > PrismExpSharedPtr
Definition: PrismExp.h:217
virtual void ProcessComposites()
Generate composites.
Represents a command-line configuration option.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:137
virtual void Process()
Write mesh to output file.
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
virtual void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
3D Evenly-spaced points on a Prism
Definition: PointsType.h:73
2D Evenly-spaced points on a Triangle
Definition: PointsType.h:70
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
Definition: Face.h:378
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
Describes the specification for a Basis.
Definition: Basis.h:50
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215