This processing module calculates the Jacobian of elements using SpatialDomains::GeomFactors and the Element::GetGeom method. For now it simply prints a list of elements which have negative Jacobian. More...

#include <ProcessTetSplit.h>

Inheritance diagram for Nektar::Utilities::ProcessTetSplit:

[legend]

Collaboration diagram for Nektar::Utilities::ProcessTetSplit:

[legend]

Public Member Functions
	ProcessTetSplit (MeshSharedPtr m)

virtual	~ProcessTetSplit ()

virtual void	Process ()
	Write mesh to output file. More...

Public Member Functions inherited from Nektar::Utilities::ProcessModule
	ProcessModule ()

	ProcessModule (FieldSharedPtr p_f)

	ProcessModule (MeshSharedPtr p_m)

Public Member Functions inherited from Nektar::Utilities::Module
	Module (FieldSharedPtr p_f)

virtual void	Process (po::variables_map &vm)=0

void	RegisterConfig (string key, string value)
	Register a configuration option with a module. More...

void	PrintConfig ()
	Print out all configuration options for a module. More...

void	SetDefaults ()
	Sets default configuration options for those which have not been set. More...

bool	GetRequireEquiSpaced (void)

void	SetRequireEquiSpaced (bool pVal)

void	EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)

	Module (MeshSharedPtr p_m)

void	RegisterConfig (string key, string value)

void	PrintConfig ()

void	SetDefaults ()

MeshSharedPtr	GetMesh ()

virtual void	ProcessVertices ()
	Extract element vertices. More...

virtual void	ProcessEdges (bool ReprocessEdges=true)
	Extract element edges. More...

virtual void	ProcessFaces (bool ReprocessFaces=true)
	Extract element faces. More...

virtual void	ProcessElements ()
	Generate element IDs. More...

virtual void	ProcessComposites ()
	Generate composites. More...

virtual void	ClearElementLinks ()

Static Public Member Functions
static boost::shared_ptr< Module >	create (MeshSharedPtr m)
	Creates an instance of this class. More...

Static Public Attributes
static ModuleKey	className

Additional Inherited Members
Protected Member Functions inherited from Nektar::Utilities::Module
	Module ()

void	ReorderPrisms (PerMap &perFaces)
	Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...

void	PrismLines (int prism, PerMap &perFaces, set< int > &prismsDone, vector< ElementSharedPtr > &line)

Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr	m_f
	Field object. More...

map< string, ConfigOption >	m_config
	List of configuration values. More...

bool	m_requireEquiSpaced

MeshSharedPtr	m_mesh
	Mesh object. More...

Detailed Description

This processing module calculates the Jacobian of elements using SpatialDomains::GeomFactors and the Element::GetGeom method. For now it simply prints a list of elements which have negative Jacobian.

Definition at line 52 of file ProcessTetSplit.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessTetSplit::ProcessTetSplit ( MeshSharedPtr m )

Definition at line 63 of file ProcessTetSplit.cpp.

References Nektar::Utilities::Module::m_config.

                                                 : ProcessModule(m)
 {
     m_config["nq"] =
         ConfigOption(false, "3", "Number of points in high order elements.");
 }

Nektar::Utilities::ProcessTetSplit::~ProcessTetSplit ( )

virtual

Definition at line 69 of file ProcessTetSplit.cpp.

70 {

71 }

Member Function Documentation

static boost::shared_ptr<Module> Nektar::Utilities::ProcessTetSplit::create ( MeshSharedPtr m )

inlinestatic

Creates an instance of this class.

Definition at line 56 of file ProcessTetSplit.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

     {
         return MemoryManager<ProcessTetSplit>::AllocateSharedPtr(m);
     }

void Nektar::Utilities::ProcessTetSplit::Process ( )

virtual

Write mesh to output file.

Implements Nektar::Utilities::Module.

Definition at line 73 of file ProcessTetSplit.cpp.

 {
     int nodeId = m_mesh->m_vertexSet.size();
 
     // Set up map which identifies edges (as pairs of vertex ids)
     // including their vertices to the offset/stride in the 3d array of
     // collapsed co-ordinates. Note that this map also includes the
     // diagonal edges along quadrilateral faces which will be used to
     // add high-order information to the split tetrahedra.
     map<ipair, ipair> edgeMap;
     map<ipair, ipair>::iterator it;
 
     int nq  = m_config["nq"].as<int>();
     int ne  = nq - 2;            // Edge interior
     int nft = (ne - 1) * ne / 2; // Face interior (triangle)
     int nfq = ne * ne;           // Face interior (quad)
     int o   = 6;
 
     // Standard prismatic edges (0->8)
     edgeMap[ipair(0, 1)] = ipair(o, 1);
     edgeMap[ipair(1, 2)] = ipair(o + 1 * ne, 1);
     edgeMap[ipair(3, 2)] = ipair(o + 2 * ne, 1);
     edgeMap[ipair(0, 3)] = ipair(o + 3 * ne, 1);
     edgeMap[ipair(0, 4)] = ipair(o + 4 * ne, 1);
     edgeMap[ipair(1, 4)] = ipair(o + 5 * ne, 1);
     edgeMap[ipair(2, 5)] = ipair(o + 6 * ne, 1);
     edgeMap[ipair(3, 5)] = ipair(o + 7 * ne, 1);
     edgeMap[ipair(4, 5)] = ipair(o + 8 * ne, 1);
 
     // Face 0 diagonals
     o += 9 * ne;
     edgeMap[ipair(0, 2)] = ipair(o, ne + 1);
     edgeMap[ipair(1, 3)] = ipair(o + ne - 1, ne - 1);
 
     // Face 2 diagonals
     o += nfq + nft;
     edgeMap[ipair(1, 5)] = ipair(o, ne + 1);
     edgeMap[ipair(2, 4)] = ipair(o + ne - 1, ne - 1);
 
     // Face 4 diagonals
     o += nfq + nft;
     edgeMap[ipair(0, 5)] = ipair(o, ne + 1);
     edgeMap[ipair(3, 4)] = ipair(o + ne - 1, ne - 1);
 
     // Default PointsType.
     LibUtilities::PointsType pt = LibUtilities::ePolyEvenlySpaced;
 
     /*
      * Split all element types into tetrahedra. This is based on the
      * algorithm found in:
      *
      * "How to Subdivide Pyramids, Prisms and Hexahedra into
      * Tetrahedra", J. Dompierre et al.
      */
 
     // Denotes a set of indices inside m_mesh->m_element[m_mesh->m_expDim-1]
     // which are
     // to be removed. These are precisely the quadrilateral boundary
     // faces which will be replaced by two triangular faces.
     set<int> toRemove;
 
     // Represents table 2 of paper; each row i represents a rotation of
     // the prism nodes such that vertex i is placed at position 0.
     static int indir[6][6] = {{0, 1, 2, 3, 4, 5},
                               {1, 2, 0, 4, 5, 3},
                               {2, 0, 1, 5, 3, 4},
                               {3, 5, 4, 0, 2, 1},
                               {4, 3, 5, 1, 0, 2},
                               {5, 4, 3, 2, 1, 0}};
 
     // Represents table 3 of paper; the first three rows are the three
     // tetrahedra if the first condition is met; the latter three rows
     // are the three tetrahedra if the second condition is met.
     static int prismTet[6][4] = {{0, 1, 2, 5},
                                  {0, 1, 5, 4},
                                  {0, 4, 5, 3},
                                  {0, 1, 2, 4},
                                  {0, 4, 2, 5},
                                  {0, 4, 5, 3}};
 
     // Represents the order of tetrahedral edges (in Nektar++ ordering).
     static int tetEdges[6][2] = {
         {0, 1}, {1, 2}, {0, 2}, {0, 3}, {1, 3}, {2, 3}};
 
     // A tetrahedron nodes -> faces map.
     static int tetFaceNodes[4][3] = {
         {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
 
     static NodeSharedPtr stdPrismNodes[6] = {
         NodeSharedPtr(new Node(0, -1, -1, -1)),
         NodeSharedPtr(new Node(1, 1, -1, -1)),
         NodeSharedPtr(new Node(2, 1, 1, -1)),
         NodeSharedPtr(new Node(3, -1, 1, -1)),
         NodeSharedPtr(new Node(4, -1, -1, 1)),
         NodeSharedPtr(new Node(5, -1, 1, 1))};
 
     // Generate equally spaced nodal points on triangle.
     Array<OneD, NekDouble> rp(nq * (nq + 1) / 2), sp(nq * (nq + 1) / 2);
     LibUtilities::PointsKey elec(nq, LibUtilities::eNodalTriEvenlySpaced);
     LibUtilities::PointsManager()[elec]->GetPoints(rp, sp);
 
     // Make a copy of the element list.
     vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
     m_mesh->m_element[m_mesh->m_expDim].clear();
 
     for (int i = 0; i < el.size(); ++i)
     {
         if (el[i]->GetConf().m_e != LibUtilities::ePrism)
         {
             m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
             continue;
         }
 
         vector<NodeSharedPtr> nodeList(6);
 
         // Map Nektar++ ordering (vertices 0,1,2,3 are base quad) to
         // paper ordering (vertices 0,1,2 are first triangular face).
         int mapPrism[6] = {0, 1, 4, 3, 2, 5};
         for (int j = 0; j < 6; ++j)
         {
             nodeList[j] = el[i]->GetVertex(mapPrism[j]);
         }
 
         // Determine minimum ID of the nodes in this prism.
         int minElId = nodeList[0]->m_id;
         int minId = 0;
         for (int j = 1; j < 6; ++j)
         {
             int curId = nodeList[j]->m_id;
             if (curId < minElId)
             {
                 minElId = curId;
                 minId   = j;
             }
         }
 
         int offset;
 
         // Split prism using paper criterion.
         int id1 = min(nodeList[indir[minId][1]]->m_id,
                       nodeList[indir[minId][5]]->m_id);
         int id2 = min(nodeList[indir[minId][2]]->m_id,
                       nodeList[indir[minId][4]]->m_id);
 
         if (id1 < id2)
         {
             offset = 0;
         }
         else if (id1 > id2)
         {
             offset = 3;
         }
         else
         {
             cerr << "Connectivity issue with prism->tet splitting." << endl;
             abort();
         }
 
         // Create local prismatic region so that co-ordinates of the
         // mapped element can be read from.
         SpatialDomains::PrismGeomSharedPtr geomLayer =
             boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(
                 el[i]->GetGeom(m_mesh->m_spaceDim));
         LibUtilities::BasisKey B0(
             LibUtilities::eOrtho_A,
             nq,
             LibUtilities::PointsKey(nq, LibUtilities::eGaussLobattoLegendre));
         LibUtilities::BasisKey B1(
             LibUtilities::eOrtho_B,
             nq,
             LibUtilities::PointsKey(nq, LibUtilities::eGaussRadauMAlpha1Beta0));
         LocalRegions::PrismExpSharedPtr qs =
             MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(
                 B0, B0, B1, geomLayer);
 
         // Get the coordinates of the high order prismatic element.
         Array<OneD, NekDouble> wsp(3 * nq * nq * nq);
         Array<OneD, Array<OneD, NekDouble> > x(3);
         x[0] = wsp;
         x[1] = wsp + 1 * nq * nq * nq;
         x[2] = wsp + 2 * nq * nq * nq;
         qs->GetCoords(x[0], x[1], x[2]);
 
         LibUtilities::BasisKey NB0(
             LibUtilities::eModified_A, nq, LibUtilities::PointsKey(nq, pt));
         LibUtilities::BasisKey NB1(
             LibUtilities::eModified_B, nq, LibUtilities::PointsKey(nq, pt));
 
         // Process face data. Initially put coordinates into equally
         // spaced nodal distribution.
         StdRegions::StdNodalPrismExpSharedPtr nodalPrism =
             MemoryManager<StdRegions::StdNodalPrismExp>::AllocateSharedPtr(
                 B0, B0, B1, LibUtilities::eNodalPrismEvenlySpaced);
 
         int nCoeffs = nodalPrism->GetNcoeffs();
         Array<OneD, NekDouble> wsp2(3 * nCoeffs);
         Array<OneD, Array<OneD, NekDouble> > xn(3);
         xn[0] = wsp2;
         xn[1] = wsp2 + 1 * nCoeffs;
         xn[2] = wsp2 + 2 * nCoeffs;
 
         for (int j = 0; j < 3; ++j)
         {
             Array<OneD, NekDouble> tmp(nodalPrism->GetNcoeffs());
             qs->FwdTrans(x[j], tmp);
             nodalPrism->ModalToNodal(tmp, xn[j]);
         }
 
         // Store map of nodes.
         map<int, int> prismVerts;
         for (int j = 0; j < 6; ++j)
         {
             prismVerts[el[i]->GetVertex(j)->m_id] = j;
         }
 
         for (int j = 0; j < 3; ++j)
         {
             vector<NodeSharedPtr> tetNodes(4);
 
             // Extract vertices for tetrahedron.
             for (int k = 0; k < 4; ++k)
             {
                 tetNodes[k] = nodeList[indir[minId][prismTet[j + offset][k]]];
             }
 
             // Add high order information to tetrahedral edges.
             for (int k = 0; k < 6; ++k)
             {
                 // Determine prismatic nodes which correspond with this
                 // edge. Apply prism map to this to get Nektar++
                 // ordering (this works since as a permutation,
                 // prismMap^2 = id).
                 int n1 = mapPrism
                     [indir[minId][prismTet[j + offset][tetEdges[k][0]]]];
                 int n2 = mapPrism
                     [indir[minId][prismTet[j + offset][tetEdges[k][1]]]];
 
                 // Find offset/stride
                 it = edgeMap.find(pair<int, int>(n1, n2));
                 if (it == edgeMap.end())
                 {
                     it = edgeMap.find(pair<int, int>(n2, n1));
                     if (it == edgeMap.end())
                     {
                         cerr << "Couldn't find prism edges " << n1 << " " << n2
                              << endl;
                         abort();
                     }
                     // Extract vertices -- reverse order.
                     for (int l = ne - 1; l >= 0; --l)
                     {
                         int pos = it->second.first + l * it->second.second;
                         tetNodes.push_back(NodeSharedPtr(new Node(
                             nodeId++, xn[0][pos], xn[1][pos], xn[2][pos])));
                     }
                 }
                 else
                 {
                     // Extract vertices -- forwards order.
                     for (int l = 0; l < ne; ++l)
                     {
                         int pos = it->second.first + l * it->second.second;
                         tetNodes.push_back(NodeSharedPtr(new Node(
                             nodeId++, xn[0][pos], xn[1][pos], xn[2][pos])));
                     }
                 }
             }
 
             // Create new tetrahedron with edge curvature.
             vector<int> tags = el[i]->GetTagList();
             ElmtConfig conf(LibUtilities::eTetrahedron, nq - 1, false, false);
             ElementSharedPtr elmt = GetElementFactory().CreateInstance(
                 LibUtilities::eTetrahedron, conf, tetNodes, tags);
 
             // Extract interior face data.
             for (int k = 0; k < 4; ++k)
             {
                 // First determine which nodes of prism are being used
                 // for this face.
                 FaceSharedPtr face = elmt->GetFace(k);
                 vector<NodeSharedPtr> triNodes(3);
 
                 for (int l = 0; l < 3; ++l)
                 {
                     NodeSharedPtr v = face->m_vertexList[l];
                     triNodes[l]     = stdPrismNodes[prismVerts[v->m_id]];
                 }
 
                 // Create a triangle with the standard nodes of the
                 // prism.
                 vector<int> tags;
                 ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);
                 ElementSharedPtr elmt = GetElementFactory().CreateInstance(
                     LibUtilities::eTriangle, conf, triNodes, tags);
                 SpatialDomains::GeometrySharedPtr triGeom = elmt->GetGeom(3);
                 triGeom->FillGeom();
                 o                 = 3 + 3 * ne;
                 face->m_curveType = LibUtilities::eNodalTriEvenlySpaced;
 
                 for (int l = 0; l < nft; ++l)
                 {
                     Array<OneD, NekDouble> tmp1(2), tmp2(3);
                     tmp1[0] = rp[o + l];
                     tmp1[1] = sp[o + l];
                     tmp2[0] = triGeom->GetCoord(0, tmp1);
                     tmp2[1] = triGeom->GetCoord(1, tmp1);
                     tmp2[2] = triGeom->GetCoord(2, tmp1);
 
                     // PhysEvaluate on prism geometry.
                     NekDouble xc = geomLayer->GetCoord(0, tmp2);
                     NekDouble yc = geomLayer->GetCoord(1, tmp2);
                     NekDouble zc = geomLayer->GetCoord(2, tmp2);
                     face->m_faceNodes.push_back(
                         NodeSharedPtr(new Node(nodeId++, xc, yc, zc)));
                 }
             }
 
             m_mesh->m_element[m_mesh->m_expDim].push_back(elmt);
         }
 
         // Now check to see if this one of the quadrilateral faces is
         // associated with a boundary condition. If it is, we split the
         // face into two triangles and mark the existing face for
         // removal.
         //
         // Note that this algorithm has significant room for improvement
         // and is likely one of the least optimal approachs - however
         // implementation is simple.
         for (int fid = 0; fid < 5; fid += 2)
         {
             int bl = el[i]->GetBoundaryLink(fid);
 
             if (bl == -1)
             {
                 continue;
             }
 
             vector<NodeSharedPtr> triNodeList(3);
             vector<int> faceNodes(3);
             vector<int> tmp;
             vector<int> tagBE;
             ElmtConfig bconf(LibUtilities::eTriangle, 1, true, true);
             ElementSharedPtr elmt;
 
             // Mark existing boundary face for removal.
             toRemove.insert(bl);
             tagBE = m_mesh->m_element[m_mesh->m_expDim - 1][bl]->GetTagList();
 
             // First loop over tets.
             for (int j = 0; j < 3; ++j)
             {
                 // Now loop over faces.
                 for (int k = 0; k < 4; ++k)
                 {
                     // Finally determine Nektar++ local node numbers for
                     // this face.
                     for (int l = 0; l < 3; ++l)
                     {
                         faceNodes[l] = mapPrism
                             [indir[minId]
                                   [prismTet[j + offset][tetFaceNodes[k][l]]]];
                     }
 
                     tmp = faceNodes;
                     sort(faceNodes.begin(), faceNodes.end());
 
                     // If this face matches a triple denoting a split
                     // quad face, add the face to the expansion list.
                     if ((fid == 0 && ((faceNodes[0] == 0 && faceNodes[1] == 1 &&
                                        faceNodes[2] == 2) ||
                                       (faceNodes[0] == 0 && faceNodes[1] == 2 &&
                                        faceNodes[2] == 3) ||
                                       (faceNodes[0] == 0 && faceNodes[1] == 1 &&
                                        faceNodes[2] == 3) ||
                                       (faceNodes[0] == 1 && faceNodes[1] == 2 &&
                                        faceNodes[2] == 3))) ||
                         (fid == 2 && ((faceNodes[0] == 1 && faceNodes[1] == 2 &&
                                        faceNodes[2] == 5) ||
                                       (faceNodes[0] == 1 && faceNodes[1] == 4 &&
                                        faceNodes[2] == 5) ||
                                       (faceNodes[0] == 1 && faceNodes[1] == 2 &&
                                        faceNodes[2] == 4) ||
                                       (faceNodes[0] == 2 && faceNodes[1] == 4 &&
                                        faceNodes[2] == 5))) ||
                         (fid == 4 && ((faceNodes[0] == 0 && faceNodes[1] == 3 &&
                                        faceNodes[2] == 5) ||
                                       (faceNodes[0] == 0 && faceNodes[1] == 4 &&
                                        faceNodes[2] == 5) ||
                                       (faceNodes[0] == 0 && faceNodes[1] == 3 &&
                                        faceNodes[2] == 4) ||
                                       (faceNodes[0] == 3 && faceNodes[1] == 4 &&
                                        faceNodes[2] == 5))))
                     {
                         triNodeList[0] = nodeList[mapPrism[tmp[0]]];
                         triNodeList[1] = nodeList[mapPrism[tmp[1]]];
                         triNodeList[2] = nodeList[mapPrism[tmp[2]]];
                         elmt = GetElementFactory().CreateInstance(
                             LibUtilities::eTriangle, bconf, triNodeList, tagBE);
                         m_mesh->m_element[m_mesh->m_expDim - 1].push_back(elmt);
                     }
                 }
             }
         }
     }
 
     // Remove 2D elements.
     vector<ElementSharedPtr> tmp;
     for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim - 1].size(); ++i)
     {
         set<int>::iterator it = toRemove.find(i);
         if (it == toRemove.end())
         {
             tmp.push_back(m_mesh->m_element[m_mesh->m_expDim - 1][i]);
         }
     }
 
     m_mesh->m_element[m_mesh->m_expDim - 1] = tmp;
 
     // Re-process mesh to eliminate duplicate vertices and edges.
     ProcessVertices();
     ProcessEdges();
     ProcessFaces();
     ProcessElements();
     ProcessComposites();
 }

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessTetSplit::className

static

Initial value:

=
    GetModuleFactory().RegisterCreatorFunction(
        ModuleKey(eProcessModule, "tetsplit"),
        ProcessTetSplit::create,
        "Split prismatic elements to tetrahedra")

Definition at line 60 of file ProcessTetSplit.h.

Public Member Functions

Static Public Member Functions

Static Public Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation