36 #include <boost/unordered_map.hpp>
55 typedef std::pair<int, int>
ipair;
60 ProcessTetSplit::create,
61 "Split prismatic elements to tetrahedra");
66 ConfigOption(
false,
"3",
"Number of points in high order elements.");
75 int nodeId =
m_mesh->m_vertexSet.size();
82 map<ipair, ipair> edgeMap;
87 int nft = (ne - 1) * ne / 2;
105 edgeMap[
ipair(1, 3)] =
ipair(o + ne - 1, ne - 1);
110 edgeMap[
ipair(2, 4)] =
ipair(o + ne - 1, ne - 1);
115 edgeMap[
ipair(3, 4)] =
ipair(o + ne - 1, ne - 1);
136 static int indir[6][6] = {{0, 1, 2, 3, 4, 5},
146 static int prismTet[6][4] = {{0, 1, 2, 5},
154 static int tetEdges[6][2] = {
155 {0, 1}, {1, 2}, {0, 2}, {0, 3}, {1, 3}, {2, 3}};
158 static int tetFaceNodes[4][3] = {
159 {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
175 vector<ElementSharedPtr> el =
m_mesh->m_element[
m_mesh->m_expDim];
178 for (
int i = 0; i < el.size(); ++i)
186 vector<NodeSharedPtr> nodeList(6);
190 int mapPrism[6] = {0, 1, 4, 3, 2, 5};
191 for (
int j = 0; j < 6; ++j)
193 nodeList[j] = el[i]->GetVertex(mapPrism[j]);
197 int minElId = nodeList[0]->m_id;
199 for (
int j = 1; j < 6; ++j)
201 int curId = nodeList[j]->m_id;
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);
227 cerr <<
"Connectivity issue with prism->tet splitting." << endl;
235 el[i]->GetGeom(
m_mesh->m_spaceDim));
246 B0, B0, B1, geomLayer);
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]);
267 int nCoeffs = nodalPrism->GetNcoeffs();
271 xn[1] = wsp2 + 1 * nCoeffs;
272 xn[2] = wsp2 + 2 * nCoeffs;
274 for (
int j = 0; j < 3; ++j)
277 qs->FwdTrans(x[j], tmp);
278 nodalPrism->ModalToNodal(tmp, xn[j]);
282 map<int, int> prismVerts;
283 for (
int j = 0; j < 6; ++j)
285 prismVerts[el[i]->GetVertex(j)->m_id] = j;
288 for (
int j = 0; j < 3; ++j)
290 vector<NodeSharedPtr> tetNodes(4);
293 for (
int k = 0; k < 4; ++k)
295 tetNodes[k] = nodeList[indir[minId][prismTet[j + offset][k]]];
299 for (
int k = 0; k < 6; ++k)
306 [indir[minId][prismTet[j + offset][tetEdges[k][0]]]];
308 [indir[minId][prismTet[j + offset][tetEdges[k][1]]]];
311 it = edgeMap.find(pair<int, int>(n1, n2));
312 if (it == edgeMap.end())
314 it = edgeMap.find(pair<int, int>(n2, n1));
315 if (it == edgeMap.end())
317 cerr <<
"Couldn't find prism edges " << n1 <<
" " << n2
322 for (
int l = ne - 1; l >= 0; --l)
324 int pos = it->second.first + l * it->second.second;
326 nodeId++, xn[0][pos], xn[1][pos], xn[2][pos])));
332 for (
int l = 0; l < ne; ++l)
334 int pos = it->second.first + l * it->second.second;
336 nodeId++, xn[0][pos], xn[1][pos], xn[2][pos])));
342 vector<int> tags = el[i]->GetTagList();
348 for (
int k = 0; k < 4; ++k)
353 vector<NodeSharedPtr> triNodes(3);
355 for (
int l = 0; l < 3; ++l)
358 triNodes[l] = stdPrismNodes[prismVerts[v->m_id]];
372 for (
int l = 0; l < nft; ++l)
377 tmp2[0] = triGeom->GetCoord(0, tmp1);
378 tmp2[1] = triGeom->GetCoord(1, tmp1);
379 tmp2[2] = triGeom->GetCoord(2, tmp1);
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(
401 for (
int fid = 0; fid < 5; fid += 2)
403 int bl = el[i]->GetBoundaryLink(fid);
410 vector<NodeSharedPtr> triNodeList(3);
411 vector<int> faceNodes(3);
419 tagBE =
m_mesh->m_element[
m_mesh->m_expDim - 1][bl]->GetTagList();
422 for (
int j = 0; j < 3; ++j)
425 for (
int k = 0; k < 4; ++k)
429 for (
int l = 0; l < 3; ++l)
431 faceNodes[l] = mapPrism
433 [prismTet[j + offset][tetFaceNodes[k][l]]]];
437 sort(faceNodes.begin(), faceNodes.end());
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))))
466 triNodeList[0] = nodeList[mapPrism[tmp[0]]];
467 triNodeList[1] = nodeList[mapPrism[tmp[1]]];
468 triNodeList[2] = nodeList[mapPrism[tmp[2]]];
471 m_mesh->m_element[
m_mesh->m_expDim - 1].push_back(elmt);
479 vector<ElementSharedPtr> tmp;
480 for (
int i = 0; i <
m_mesh->m_element[
m_mesh->m_expDim - 1].size(); ++i)
483 if (it == toRemove.end())
485 tmp.push_back(
m_mesh->m_element[
m_mesh->m_expDim - 1][i]);
Basic information about an element.
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.
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 .
MeshSharedPtr m_mesh
Mesh object.
ElementFactory & GetElementFactory()
Represents a point in the domain.
Gauss Radau pinned at x=-1, .
virtual void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
Principle Orthogonal Functions .
1D Evenly-spaced points using Lagrange polynomial
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 .
boost::shared_ptr< Node > NodeSharedPtr
PointsManagerT & PointsManager(void)
Principle Orthogonal Functions .
Defines a specification for a set of points.
virtual ~ProcessTetSplit()
boost::shared_ptr< PrismExp > PrismExpSharedPtr
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.
virtual void Process()
Write mesh to output file.
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
virtual void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
3D Evenly-spaced points on a Prism
2D Evenly-spaced points on a Triangle
boost::shared_ptr< Element > ElementSharedPtr
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
boost::shared_ptr< Geometry > GeometrySharedPtr
Describes the specification for a Basis.
1D Gauss-Lobatto-Legendre quadrature points
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.