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