ProcessTetSplit.cpp

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////////
//
//  File: ProcessTetSplit.cpp
//
//  For more information, please see: http://www.nektar.info/
//
//  The MIT License
//
//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
//  Department of Aeronautics, Imperial College London (UK), and Scientific
//  Computing and Imaging Institute, University of Utah (USA).
//
//  License for the specific language governing rights and limitations under
//  Permission is hereby granted, free of charge, to any person obtaining a
//  copy of this software and associated documentation files (the "Software"),
//  to deal in the Software without restriction, including without limitation
//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
//  and/or sell copies of the Software, and to permit persons to whom the
//  Software is furnished to do so, subject to the following conditions:
//
//  The above copyright notice and this permission notice shall be included
//  in all copies or substantial portions of the Software.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
//  DEALINGS IN THE SOFTWARE.
//
//  Description: Split prisms -> tets
//
////////////////////////////////////////////////////////////////////////////////

#include <string>
using namespace std;

#include "../MeshElements.h"
#include "../InputModules/InputNek.h"
#include "ProcessTetSplit.h"

#include <StdRegions/StdNodalPrismExp.h>
#include <LocalRegions/PrismExp.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Foundations/ManagerAccess.h>

#include <boost/unordered_map.hpp>

namespace Nektar
{
    namespace Utilities
    {
        typedef std::pair<int,int> ipair;

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

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

        ProcessTetSplit::~ProcessTetSplit()
        {

        }

        void ProcessTetSplit::Process()
        {
            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();
        }
    }
}