Nektar::NekMeshUtils::FaceMesh Class Reference

class for surface meshes on individual surfaces (paramter plane meshes) More...

#include <FaceMesh.h>

Collaboration diagram for Nektar::NekMeshUtils::FaceMesh:

Public Member Functions
	FaceMesh (const int id, MeshSharedPtr m, CADSurfSharedPtr cad, OctreeSharedPtr oct, const std::map< int, CurveMeshSharedPtr > &cmeshes)
	Default constructor. More...
	FaceMesh (const int id, MeshSharedPtr m, CADSurfSharedPtr cad, OctreeSharedPtr oct, const std::map< int, CurveMeshSharedPtr > &cmeshes, const NekDouble b)
	constructor for building with boundary layer More...
void	Mesh ()
	mesh exectuation command More...

Private Member Functions
void	Stretching ()
	Calculate the paramter plane streching factor. More...
void	Smoothing ()
	performs node smoothing on face More...
void	DiagonalSwap ()
	performs diagonal swapping of edges More...
void	BuildLocalMesh ()
	build a local version of mesh elements More...
void	OptimiseLocalMesh ()
	function which calls the optimisation routines More...
bool	Validate ()
	Validate the surface mesh base on the octree and real dimensions of the edges. More...
void	OrientateCurves ()
	Get the boundries of the surface and extracts the nodes from the curve meshes in the correct order. More...
void	AddNewPoint (Array< OneD, NekDouble > uv)
	adds a new stiener point to the triangulation for meshing More...
void	MakeBL ()
	adds a quad layer around any interior loops More...

Private Attributes
MeshSharedPtr	m_mesh
	mesh pointer More...
CADSurfSharedPtr	m_cadsurf
	CAD surface. More...
OctreeSharedPtr	m_octree
	Octree object. More...
std::map< int, CurveMeshSharedPtr >	m_curvemeshes
	Map of the curve meshes which bound the surfaces. More...
std::vector< EdgeLoop >	m_edgeloops
	data structure containing the edges, their order and oreientation for the surface More...
int	m_id
	id of the surface mesh More...
std::vector< std::vector < NodeSharedPtr > >	orderedLoops
	list of boundary nodes in their order loops More...
std::vector< NodeSharedPtr >	m_stienerpoints
	list of stiener points in the triangulation More...
NekDouble	m_str
	pplane stretching More...
std::vector< std::vector < NodeSharedPtr > >	m_connec
	triangle connectiviities More...
NodeSet	m_localNodes
	local set of nodes More...
EdgeSet	m_localEdges
	local set of edges More...
std::vector< ElementSharedPtr >	m_localElements
	local list of elements More...
NekDouble	m_bl
	boundary layer thickness More...
bool	m_makebl
	should build boundary layer More...
std::vector< std::pair < NodeSharedPtr, NodeSharedPtr > >	blpairs
	list of node links between node on loop and its corresponding interior node in quads More...

Friends
class	MemoryManager< FaceMesh >

Detailed Description

class for surface meshes on individual surfaces (paramter plane meshes)

Definition at line 53 of file FaceMesh.h.

Constructor & Destructor Documentation

Nektar::NekMeshUtils::FaceMesh::FaceMesh ( const int  id, MeshSharedPtr  m, CADSurfSharedPtr  cad, OctreeSharedPtr  oct, const std::map< int, CurveMeshSharedPtr > &  cmeshes  )

inline

Default constructor.

Definition at line 61 of file FaceMesh.h.

References m_cadsurf, m_edgeloops, and m_makebl.

        : m_mesh(m), m_cadsurf(cad), m_octree(oct), m_curvemeshes(cmeshes),
          m_id(id)

    {
        m_edgeloops = m_cadsurf->GetEdges();
        m_makebl    = false;
    };

Nektar::NekMeshUtils::FaceMesh::FaceMesh ( const int  id, MeshSharedPtr  m, CADSurfSharedPtr  cad, OctreeSharedPtr  oct, const std::map< int, CurveMeshSharedPtr > &  cmeshes, const NekDouble  b  )

inline

constructor for building with boundary layer

Definition at line 77 of file FaceMesh.h.

References m_cadsurf, m_edgeloops, and m_makebl.

        : m_mesh(m), m_cadsurf(cad), m_octree(oct), m_curvemeshes(cmeshes),
          m_id(id), m_bl(b)

    {
        m_edgeloops = m_cadsurf->GetEdges();
        m_makebl    = true;
    };

Member Function Documentation

void Nektar::NekMeshUtils::FaceMesh::AddNewPoint ( Array< OneD, NekDouble >  uv )

private

adds a new stiener point to the triangulation for meshing

Definition at line 976 of file FaceMesh.cpp.

{
    // adds a new point but checks that there are no other points nearby first
    Array<OneD, NekDouble> np = m_cadsurf->P(uv);
    NekDouble npDelta         = m_octree->Query(np);

    NodeSharedPtr n = boost::shared_ptr<Node>(
        new Node(m_mesh->m_numNodes++, np[0], np[1], np[2]));

    bool add = true;

    for (int i = 0; i < orderedLoops.size(); i++)
    {
        for (int j = 0; j < orderedLoops[i].size(); j++)
        {
            NekDouble r = orderedLoops[i][j]->Distance(n);

            if (r < npDelta / 2.0)
            {
                add = false;
                break;
            }
        }
    }

    if (add)
    {
        for (int i = 0; i < m_stienerpoints.size(); i++)
        {
            NekDouble r = m_stienerpoints[i]->Distance(n);

            if (r < npDelta / 2.0)
            {
                add = false;
                break;
            }
        }
    }

    if (add)
    {
        n->SetCADSurf(m_id, m_cadsurf, uv);
        m_stienerpoints.push_back(n);
    }
}

void Nektar::NekMeshUtils::FaceMesh::BuildLocalMesh ( )

private

build a local version of mesh elements

Definition at line 745 of file FaceMesh.cpp.

References Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::NekMeshUtils::GetElementFactory(), and Nektar::iterator.

{
    /*************************
    // build a local set of nodes edges and elemenets for optimstaion prior to
    putting them into m_mesh
    */
    int tricomp = m_mesh->m_numcomp++;

    // first build quads is bl surface
    if (m_makebl)
    {
        int quadcomp = m_mesh->m_numcomp++;
        for (int i = 0; i < blpairs.size() - 1;
             i++) // this wont work for more than one sym plane loop
        {
            ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
            vector<NodeSharedPtr> ns;
            ns.push_back(blpairs[i].first);
            ns.push_back(blpairs[i].second);
            ns.push_back(blpairs[i + 1].second);
            ns.push_back(blpairs[i + 1].first);
            vector<int> tags;
            tags.push_back(quadcomp);
            ElementSharedPtr E = GetElementFactory().CreateInstance(
                LibUtilities::eQuadrilateral, conf, ns, tags);
            E->CADSurfId = m_id;
            m_localElements.push_back(E);
        }
        {
            ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
            vector<NodeSharedPtr> ns;
            ns.push_back(blpairs.back().first);
            ns.push_back(blpairs.back().second);
            ns.push_back(blpairs[0].second);
            ns.push_back(blpairs[0].first);
            vector<int> tags;
            tags.push_back(quadcomp);
            ElementSharedPtr E = GetElementFactory().CreateInstance(
                LibUtilities::eQuadrilateral, conf, ns, tags);
            E->CADSurfId = m_id;
            m_localElements.push_back(E);
        }
        for (int i = 0; i < m_localElements.size(); i++)
        {
            vector<NodeSharedPtr> nods = m_localElements[i]->GetVertexList();
            for (int j = 0; j < nods.size(); j++)
            {
                // nodes are already unique some will insert some wont
                m_localNodes.insert(nods[j]);
            }
            vector<EdgeSharedPtr> edgs = m_localElements[i]->GetEdgeList();
            for (int j = 0; j < edgs.size(); j++)
            {
                // look for edge in m_mesh edgeset from curves
                EdgeSet::iterator s = m_mesh->m_edgeSet.find(edgs[j]);
                if (!(s == m_mesh->m_edgeSet.end()))
                {
                    edgs[j] = *s;
                    m_localElements[i]->SetEdge(j, edgs[j]);
                }

                pair<EdgeSet::iterator, bool> test =
                    m_localEdges.insert(edgs[j]);

                if (test.second)
                {
                    (*test.first)
                        ->m_elLink.push_back(
                            pair<ElementSharedPtr, int>(m_localElements[i], j));
                }
                else
                {
                    m_localElements[i]->SetEdge(j, *test.first);
                    (*test.first)
                        ->m_elLink.push_back(
                            pair<ElementSharedPtr, int>(m_localElements[i], j));
                }
            }
            m_localElements[i]->SetId(i);
        }
    }

    for (int i = 0; i < m_connec.size(); i++)
    {
        ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);

        vector<int> tags;
        tags.push_back(tricomp);
        ElementSharedPtr E = GetElementFactory().CreateInstance(
            LibUtilities::eTriangle, conf, m_connec[i], tags);
        E->CADSurfId = m_id;

        vector<NodeSharedPtr> nods = E->GetVertexList();
        for (int j = 0; j < nods.size(); j++)
        {
            // nodes are already unique some will insert some wont
            m_localNodes.insert(nods[j]);
        }
        vector<EdgeSharedPtr> edgs = E->GetEdgeList();
        for (int j = 0; j < edgs.size(); j++)
        {
            // look for edge in m_mesh edgeset from curves
            EdgeSet::iterator s = m_mesh->m_edgeSet.find(edgs[j]);
            if (!(s == m_mesh->m_edgeSet.end()))
            {
                edgs[j] = *s;
                E->SetEdge(j, edgs[j]);
            }

            pair<EdgeSet::iterator, bool> test = m_localEdges.insert(edgs[j]);

            if (test.second)
            {
                (*test.first)
                    ->m_elLink.push_back(pair<ElementSharedPtr, int>(E, j));
            }
            else
            {
                E->SetEdge(j, *test.first);
                (*test.first)
                    ->m_elLink.push_back(pair<ElementSharedPtr, int>(E, j));
            }
        }
        E->SetId(m_localElements.size());
        m_localElements.push_back(E);
    }
}

void Nektar::NekMeshUtils::FaceMesh::DiagonalSwap ( )

private

performs diagonal swapping of edges

TODO fix this bit of code which figures out the node defect, doesnt work on quads or relfex angles

Definition at line 363 of file FaceMesh.cpp.

References ASSERTL0, Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::NekMeshUtils::GetElementFactory(), and Nektar::iterator.

{
    /// TODO fix this bit of code which figures out the node defect, doesnt work
    /// on quads or relfex angles
    map<int, int> idealConnec;
    map<int, int> actualConnec;
    map<int, vector<EdgeSharedPtr> > nodetoedge;
    // figure out ideal node count and actual node count
    EdgeSet::iterator eit;
    for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
    {
        nodetoedge[(*eit)->m_n1->m_id].push_back(*eit);
        nodetoedge[(*eit)->m_n2->m_id].push_back(*eit);
    }
    NodeSet::iterator nit;
    for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
    {
        if ((*nit)->GetNumCadCurve() == 0)
        {
            // node is interior
            idealConnec[(*nit)->m_id] = 6;
        }
        else
        {
            // need to identify the two other nodes on the boundary to find
            // interior angle
            vector<NodeSharedPtr> ns;
            vector<EdgeSharedPtr> e = nodetoedge[(*nit)->m_id];
            for (int i = 0; i < e.size(); i++)
            {
                if (!e[i]->onCurve)
                    continue; // the linking nodes are not going to exist on
                              // interior edges

                if (e[i]->m_n1 == (*nit))
                    ns.push_back(e[i]->m_n2);
                else
                    ns.push_back(e[i]->m_n1);
            }
            ASSERTL0(ns.size() == 2,
                     "failed to find 2 nodes in the angle system");

            idealConnec[(*nit)->m_id] =
                ceil((*nit)->Angle(ns[0], ns[1]) / 3.142 * 3) + 1;
        }
    }
    for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
    {
        actualConnec[(*nit)->m_id] = nodetoedge[(*nit)->m_id].size();
    }

    // edgeswapping fun times
    // perfrom edge swap based on node defect and then angle
    for (int q = 0; q < 4; q++)
    {
        int edgesStart = m_localEdges.size();
        EdgeSet edges = m_localEdges;
        m_localEdges.clear();

        int swappedEdges = 0;

        EdgeSet::iterator it;

        for (it = edges.begin(); it != edges.end(); it++)
        {
            EdgeSharedPtr e = *it;

            if (e->m_elLink.size() != 2)
            {
                m_localEdges.insert(e);
                continue;
            }
            if (e->m_elLink[0].first->GetConf().m_e ==
                    LibUtilities::eQuadrilateral ||
                e->m_elLink[1].first->GetConf().m_e ==
                    LibUtilities::eQuadrilateral)
            {
                m_localEdges.insert(e);
                continue;
            }

            ElementSharedPtr tri1 = e->m_elLink[0].first;
            ElementSharedPtr tri2 = e->m_elLink[1].first;

            NodeSharedPtr n1 = e->m_n1;
            NodeSharedPtr n2 = e->m_n2;

            vector<NodeSharedPtr> nt = tri1->GetVertexList();

            // identify node a,b,c,d of the swapping
            NodeSharedPtr A, B, C, D;
            if (nt[0] != n1 && nt[0] != n2)
            {
                C = nt[0];
                B = nt[1];
                A = nt[2];
            }
            else if (nt[1] != n1 && nt[1] != n2)
            {
                C = nt[1];
                B = nt[2];
                A = nt[0];
            }
            else if (nt[2] != n1 && nt[2] != n2)
            {
                C = nt[2];
                B = nt[0];
                A = nt[1];
            }
            else
            {
                ASSERTL0(false, "failed to identify verticies in tri1");
            }

            nt = tri2->GetVertexList();

            if (nt[0] != n1 && nt[0] != n2)
            {
                D = nt[0];
            }
            else if (nt[1] != n1 && nt[1] != n2)
            {
                D = nt[1];
            }
            else if (nt[2] != n1 && nt[2] != n2)
            {
                D = nt[2];
            }
            else
            {
                ASSERTL0(false, "failed to identify verticies in tri2");
            }

            // determine signed area of alternate config
            Array<OneD, NekDouble> ai, bi, ci, di;
            ai = A->GetCADSurfInfo(m_id);
            bi = B->GetCADSurfInfo(m_id);
            ci = C->GetCADSurfInfo(m_id);
            di = D->GetCADSurfInfo(m_id);

            NekDouble CDA, CBD;

            CDA = 0.5 * (-di[0] * ci[1] + ai[0] * ci[1] + ci[0] * di[1] -
                         ai[0] * di[1] - ci[0] * ai[1] + di[0] * ai[1]);

            CBD = 0.5 * (-bi[0] * ci[1] + di[0] * ci[1] + ci[0] * bi[1] -
                         di[0] * bi[1] - ci[0] * di[1] + bi[0] * di[1]);

            // if signed area of the swapping triangles is less than zero
            // that configuration is invalid and swap cannot be performed
            if (!(CDA > 0.001 && CBD > 0.001))
            {
                m_localEdges.insert(e);
                continue;
            }

            bool swap = false; // assume do not swap

            if (q < 2)
            {
                int nodedefectbefore = 0;
                nodedefectbefore +=
                    abs(actualConnec[A->m_id] - idealConnec[A->m_id]);
                nodedefectbefore +=
                    abs(actualConnec[B->m_id] - idealConnec[B->m_id]);
                nodedefectbefore +=
                    abs(actualConnec[C->m_id] - idealConnec[C->m_id]);
                nodedefectbefore +=
                    abs(actualConnec[D->m_id] - idealConnec[D->m_id]);

                int nodedefectafter = 0;
                nodedefectafter +=
                    abs(actualConnec[A->m_id] - 1 - idealConnec[A->m_id]);
                nodedefectafter +=
                    abs(actualConnec[B->m_id] - 1 - idealConnec[B->m_id]);
                nodedefectafter +=
                    abs(actualConnec[C->m_id] + 1 - idealConnec[C->m_id]);
                nodedefectafter +=
                    abs(actualConnec[D->m_id] + 1 - idealConnec[D->m_id]);

                if (nodedefectafter < nodedefectbefore)
                {
                    swap = true;
                }
            }
            else
            {
                NekDouble minanglebefore = C->Angle(A, B);
                minanglebefore           = min(minanglebefore, A->Angle(B, C));
                minanglebefore           = min(minanglebefore, B->Angle(A, C));
                minanglebefore           = min(minanglebefore, B->Angle(A, D));
                minanglebefore           = min(minanglebefore, A->Angle(B, D));
                minanglebefore           = min(minanglebefore, D->Angle(A, B));

                NekDouble minangleafter = C->Angle(B, D);
                minangleafter           = min(minangleafter, D->Angle(B, C));
                minangleafter           = min(minangleafter, B->Angle(C, D));
                minangleafter           = min(minangleafter, C->Angle(A, D));
                minangleafter           = min(minangleafter, A->Angle(C, D));
                minangleafter           = min(minangleafter, D->Angle(A, C));

                if (minangleafter > minanglebefore)
                {
                    swap = true;
                }
            }

            if (swap)
            {
                actualConnec[A->m_id]--;
                actualConnec[B->m_id]--;
                actualConnec[C->m_id]++;
                actualConnec[D->m_id]++;

                // make the 4 other edges
                EdgeSharedPtr CA, AD, DB, BC, CAt, ADt, DBt, BCt;
                CAt = boost::shared_ptr<Edge>(new Edge(C, A));
                ADt = boost::shared_ptr<Edge>(new Edge(A, D));
                DBt = boost::shared_ptr<Edge>(new Edge(D, B));
                BCt = boost::shared_ptr<Edge>(new Edge(B, C));

                vector<EdgeSharedPtr> es = tri1->GetEdgeList();
                for (int i = 0; i < 3; i++)
                {
                    if (es[i] == CAt)
                    {
                        CA = es[i];
                    }
                    if (es[i] == BCt)
                    {
                        BC = es[i];
                    }
                }
                es = tri2->GetEdgeList();
                for (int i = 0; i < 3; i++)
                {
                    if (es[i] == DBt)
                    {
                        DB = es[i];
                    }
                    if (es[i] == ADt)
                    {
                        AD = es[i];
                    }
                }

                // now sort out links for the 4 edges surrounding the patch
                vector<pair<ElementSharedPtr, int> > links;

                links = CA->m_elLink;
                CA->m_elLink.clear();
                for (int i = 0; i < links.size(); i++)
                {
                    if (links[i].first->GetId() == tri1->GetId())
                        continue;
                    CA->m_elLink.push_back(links[i]);
                }

                links = BC->m_elLink;
                BC->m_elLink.clear();
                for (int i = 0; i < links.size(); i++)
                {
                    if (links[i].first->GetId() == tri1->GetId())
                        continue;
                    BC->m_elLink.push_back(links[i]);
                }

                links = AD->m_elLink;
                AD->m_elLink.clear();
                for (int i = 0; i < links.size(); i++)
                {
                    if (links[i].first->GetId() == tri2->GetId())
                        continue;
                    AD->m_elLink.push_back(links[i]);
                }

                links = DB->m_elLink;
                DB->m_elLink.clear();
                for (int i = 0; i < links.size(); i++)
                {
                    if (links[i].first->GetId() == tri2->GetId())
                        continue;
                    DB->m_elLink.push_back(links[i]);
                }

                EdgeSharedPtr newe = boost::shared_ptr<Edge>(new Edge(C, D));

                vector<NodeSharedPtr> t1, t2;
                t1.push_back(B);
                t1.push_back(D);
                t1.push_back(C);
                t2.push_back(A);
                t2.push_back(C);
                t2.push_back(D);

                ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);
                vector<int> tags = tri1->GetTagList();

                int id1 = tri1->GetId();
                int id2 = tri2->GetId();

                ElementSharedPtr ntri1 = GetElementFactory().CreateInstance(
                    LibUtilities::eTriangle, conf, t1, tags);
                tags                   = tri2->GetTagList();
                ElementSharedPtr ntri2 = GetElementFactory().CreateInstance(
                    LibUtilities::eTriangle, conf, t2, tags);

                ntri1->SetId(id1);
                ntri2->SetId(id2);
                ntri1->CADSurfId = m_id;
                ntri2->CADSurfId = m_id;

                vector<EdgeSharedPtr> t1es = ntri1->GetEdgeList();
                for (int i = 0; i < 3; i++)
                {
                    if (t1es[i] == DB)
                    {
                        ntri1->SetEdge(i, DB);
                        DB->m_elLink.push_back(
                            pair<ElementSharedPtr, int>(ntri1, i));
                    }
                    else if (t1es[i] == BC)
                    {
                        ntri1->SetEdge(i, BC);
                        BC->m_elLink.push_back(
                            pair<ElementSharedPtr, int>(ntri1, i));
                    }
                    else if (t1es[i] == newe)
                    {
                        ntri1->SetEdge(i, newe);
                        newe->m_elLink.push_back(
                            pair<ElementSharedPtr, int>(ntri1, i));
                    }
                    else
                    {
                        ASSERTL0(false, "weird edge in new tri 1");
                    }
                }
                vector<EdgeSharedPtr> t2es = ntri2->GetEdgeList();
                for (int i = 0; i < 3; i++)
                {
                    if (t2es[i] == CA)
                    {
                        ntri2->SetEdge(i, CA);
                        CA->m_elLink.push_back(
                            pair<ElementSharedPtr, int>(ntri2, i));
                    }
                    else if (t2es[i] == AD)
                    {
                        ntri2->SetEdge(i, AD);
                        AD->m_elLink.push_back(
                            pair<ElementSharedPtr, int>(ntri2, i));
                    }
                    else if (t2es[i] == newe)
                    {
                        ntri2->SetEdge(i, newe);
                        newe->m_elLink.push_back(
                            pair<ElementSharedPtr, int>(ntri2, i));
                    }
                    else
                    {
                        ASSERTL0(false, "weird edge in new tri 2");
                    }
                }

                m_localEdges.insert(newe);

                m_localElements[id1] = ntri1;
                m_localElements[id2] = ntri2;

                swappedEdges++;
            }
            else
            {
                m_localEdges.insert(e);
            }
        }

        ASSERTL0(m_localEdges.size() == edgesStart, "mismatch edge count");
    }
}

void Nektar::NekMeshUtils::FaceMesh::MakeBL ( )

private

adds a quad layer around any interior loops

Definition at line 147 of file FaceMesh.cpp.

References ASSERTL0.

{
    for (int i = 1; i < orderedLoops.size(); i++) // dont do this to first loop
    {
        for (int j = 0; j < orderedLoops[i].size(); j++)
        {
            // for each of the nodes make a new node which exists off the
            // surface
            vector<pair<int, CADSurfSharedPtr> > surfs =
                orderedLoops[i][j]->GetCADSurfs();
            ASSERTL0(surfs.size() > 1,
                     "point does not have enough surfs to make quad");

            Array<OneD, NekDouble> AN(3, 0.0);
            // make a averaged normal ignoring surface mid
            for (int s = 0; s < surfs.size(); s++)
            {
                if (surfs[s].first == m_id)
                    continue; // does not contribute to norm

                Array<OneD, NekDouble> uv =
                    orderedLoops[i][j]->GetCADSurfInfo(surfs[s].first);
                Array<OneD, NekDouble> N = surfs[s].second->N(uv);
                for (int k = 0; k < 3; k++)
                    AN[k] += N[k];
            }

            // renormalise normal
            NekDouble mag = sqrt(AN[0] * AN[0] + AN[1] * AN[1] + AN[2] * AN[2]);
            for (int k = 0; k < 3; k++)
                AN[k] /= mag;

            Array<OneD, NekDouble> loc = orderedLoops[i][j]->GetLoc();
            Array<OneD, NekDouble> tp(3);
            for (int k = 0; k < 3; k++)
                tp[k] = m_bl * AN[k] + loc[k];
            // project tp onto to surface to get new point
            Array<OneD, NekDouble> uv(2);
            m_cadsurf->ProjectTo(tp, uv);
            NodeSharedPtr nn = boost::shared_ptr<Node>(
                new Node(m_mesh->m_numNodes++, tp[0], tp[1], tp[2]));
            nn->SetCADSurf(m_id, m_cadsurf, uv);

            blpairs.push_back(
                pair<NodeSharedPtr, NodeSharedPtr>(orderedLoops[i][j], nn));
            // place the new node into ordered loops for the surface
            // triangulation to work With
            orderedLoops[i][j] = nn;
        }
    }
}

void Nektar::NekMeshUtils::FaceMesh::Mesh

(

)

mesh exectuation command

Definition at line 49 of file FaceMesh.cpp.

References ASSERTL0, and Nektar::iterator.

{
    Stretching();

    OrientateCurves();

    if (m_makebl)
    {
        MakeBL();
    }

    int numPoints = 0;
    for (int i = 0; i < orderedLoops.size(); i++)
    {
        numPoints += orderedLoops[i].size();
    }

    stringstream ss;
    ss << "3 points required for triangulation, " << numPoints << " in loop"
       << endl;
    ss << "curves: ";
    for (int i = 0; i < m_edgeloops.size(); i++)
    {
        for (int j = 0; j < m_edgeloops[i].edges.size(); j++)
        {
            ss << m_edgeloops[i].edges[j]->GetId() << " ";
        }
    }

    ASSERTL0(numPoints > 2, ss.str());

    // create interface to triangle thirdparty library
    TriangleInterfaceSharedPtr pplanemesh =
        MemoryManager<TriangleInterface>::AllocateSharedPtr();

    vector<Array<OneD, NekDouble> > centers;
    for (int i = 0; i < m_edgeloops.size(); i++)
    {
        centers.push_back(m_edgeloops[i].center);
    }

    pplanemesh->Assign(orderedLoops, centers, m_id, m_str);

    pplanemesh->Mesh();

    pplanemesh->Extract(m_connec);

    bool repeat     = true;
    int meshcounter = 1;

    while (repeat)
    {
        repeat = Validate();
        if (!repeat)
        {
            break;
        }
        m_connec.clear();
        pplanemesh->AssignStiener(m_stienerpoints);
        pplanemesh->Mesh();
        pplanemesh->Extract(m_connec);
        meshcounter++;
    }

    BuildLocalMesh();

    OptimiseLocalMesh();

    // clear local element links
    EdgeSet::iterator eit;
    for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
    {
        (*eit)->m_elLink.clear();
    }

    // make new elements and add to list from list of nodes and connectivity
    // from triangle
    for (int i = 0; i < m_localElements.size(); i++)
    {
        vector<EdgeSharedPtr> e = m_localElements[i]->GetEdgeList();
        for (int j = 0; j < e.size(); j++)
        {
            e[j]->m_elLink.clear();
        }
        m_mesh->m_element[m_mesh->m_expDim].push_back(m_localElements[i]);
    }

    cout << "\r                                                                "
            "                             ";
    cout << scientific << "\r\t\tFace " << m_id << endl
         << "\t\t\tNodes: " << m_localNodes.size() << endl
         << "\t\t\tEdges: " << m_localEdges.size() << endl
         << "\t\t\tTriangles: " << m_localElements.size() << endl
         << "\t\t\tLoops: " << m_edgeloops.size() << endl
         << "\t\t\tSTR: " << m_str << endl
         << endl;
}

void Nektar::NekMeshUtils::FaceMesh::OptimiseLocalMesh ( )

private

function which calls the optimisation routines

Definition at line 199 of file FaceMesh.cpp.

{
    DiagonalSwap();

    Smoothing();

    DiagonalSwap();

    Smoothing();
}

void Nektar::NekMeshUtils::FaceMesh::OrientateCurves ( )

private

Get the boundries of the surface and extracts the nodes from the curve meshes in the correct order.

Definition at line 1022 of file FaceMesh.cpp.

References ASSERTL0.

{
    // create list of bounding loop nodes
    for (int i = 0; i < m_edgeloops.size(); i++)
    {
        vector<NodeSharedPtr> cE;
        for (int j = 0; j < m_edgeloops[i].edges.size(); j++)
        {
            int cid = m_edgeloops[i].edges[j]->GetId();
            vector<NodeSharedPtr> edgePoints =
                m_curvemeshes[cid]->GetMeshPoints();

            int numPoints = m_curvemeshes[cid]->GetNumPoints();

            if (m_edgeloops[i].edgeo[j] == 0)
            {
                for (int k = 0; k < numPoints - 1; k++)
                {
                    cE.push_back(edgePoints[k]);
                }
            }
            else
            {
                for (int k = numPoints - 1; k > 0; k--)
                {
                    cE.push_back(edgePoints[k]);
                }
            }
        }
        orderedLoops.push_back(cE);
    }

    // loops made need to orientate on which is biggest and define holes
    for (int i = 0; i < orderedLoops.size(); i++)
    {
        NekDouble area = 0.0;
        for (int j = 0; j < orderedLoops[i].size() - 1; j++)
        {
            Array<OneD, NekDouble> n1info, n2info;
            n1info = orderedLoops[i][j]->GetCADSurfInfo(m_id);
            n2info = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);

            area += -n2info[1] * (n2info[0] - n1info[0]) +
                    n1info[0] * (n2info[1] - n1info[1]);
        }
        area *= 0.5;
        m_edgeloops[i].area = area;
    }

    int ct = 0;

    do
    {
        ct = 0;
        for (int i = 0; i < m_edgeloops.size() - 1; i++)
        {
            if (fabs(m_edgeloops[i].area) < fabs(m_edgeloops[i + 1].area))
            {
                // swap
                vector<NodeSharedPtr> orderedlooptmp = orderedLoops[i];
                EdgeLoop edgeLoopstmp                = m_edgeloops[i];

                orderedLoops[i] = orderedLoops[i + 1];
                m_edgeloops[i]  = m_edgeloops[i + 1];

                orderedLoops[i + 1] = orderedlooptmp;
                m_edgeloops[i + 1]  = edgeLoopstmp;

                ct += 1;
            }
        }

    } while (ct > 0);

    for (int i = 0; i < orderedLoops.size(); i++)
    {
        NodeSharedPtr n1, n2;

        n1 = orderedLoops[i][0];
        n2 = orderedLoops[i][1];

        Array<OneD, NekDouble> n1info, n2info;
        n1info = n1->GetCADSurfInfo(m_id);
        n2info = n2->GetCADSurfInfo(m_id);

        Array<OneD, NekDouble> N(2);
        NekDouble mag = sqrt((n1info[0] - n2info[0]) * (n1info[0] - n2info[0]) +
                             (n1info[1] - n2info[1]) * (n1info[1] - n2info[1]));
        ASSERTL0(mag > 1e-30, "infinity");
        N[0] = -1.0 * (n2info[1] - n1info[1]) / mag;
        N[1] = (n2info[0] - n1info[0]) / mag;

        Array<OneD, NekDouble> P(2);
        P[0] = (n1info[0] + n2info[0]) / 2.0 + 1e-6 * N[0];
        P[1] = (n1info[1] + n2info[1]) / 2.0 + 1e-6 * N[1];

        // now test to see if p is inside or outside the shape
        // vector to the right
        int intercepts = 0;
        for (int j = 0; j < orderedLoops[i].size() - 1; j++)
        {
            Array<OneD, NekDouble> nt1, nt2;
            nt1 = orderedLoops[i][j]->GetCADSurfInfo(m_id);
            nt2 = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);

            if (fabs(nt2[1] - nt1[1]) < 1e-30)
                continue;

            NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
            NekDouble S   = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;

            if (!(lam < 0) && !(lam > 1) && S > 0)
            {
                intercepts++;
            }
        }
        {
            Array<OneD, NekDouble> nt1, nt2;
            nt1 = orderedLoops[i].back()->GetCADSurfInfo(m_id);
            nt2 = orderedLoops[i][0]->GetCADSurfInfo(m_id);

            if (fabs(nt2[1] - nt1[1]) < 1e-30)
                continue;

            NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
            NekDouble S   = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;

            if (!(lam < 0) && !(lam > 1) && S > 0)
            {
                intercepts++;
            }
        }
        if (intercepts % 2 == 0)
        {
            P[0]       = (n1info[0] + n2info[0]) / 2.0 - 1e-6 * N[0];
            P[1]       = (n1info[1] + n2info[1]) / 2.0 - 1e-6 * N[1];
            intercepts = 0;
            for (int j = 0; j < orderedLoops[i].size() - 1; j++)
            {
                Array<OneD, NekDouble> nt1, nt2;
                nt1 = orderedLoops[i][j]->GetCADSurfInfo(m_id);
                nt2 = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);

                if (fabs(nt2[1] - nt1[1]) < 1e-30)
                    continue;

                NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
                NekDouble S   = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;

                if (!(lam < 0) && !(lam > 1) && S > 0)
                {
                    intercepts++;
                }
            }
            {
                Array<OneD, NekDouble> nt1, nt2;
                nt1 = orderedLoops[i].back()->GetCADSurfInfo(m_id);
                nt2 = orderedLoops[i][0]->GetCADSurfInfo(m_id);

                if (fabs(nt2[1] - nt1[1]) < 1e-30)
                    continue;

                NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
                NekDouble S   = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;

                if (!(lam < 0) && !(lam > 1) && S > 0)
                {
                    intercepts++;
                }
            }
            if (intercepts % 2 == 0)
            {
                cerr << "still failed to find point inside loop" << endl;
            }
        }

        m_edgeloops[i].center = P;
    }

    if (m_edgeloops[0].area < 0) // reverse the first uvLoop
    {
        vector<NodeSharedPtr> tmp = orderedLoops[0];
        reverse(tmp.begin(), tmp.end());
        orderedLoops[0] = tmp;
    }

    for (int i = 1; i < orderedLoops.size(); i++)
    {
        if (m_edgeloops[i].area > 0) // reverse the loop
        {
            vector<NodeSharedPtr> tmp = orderedLoops[i];
            reverse(tmp.begin(), tmp.end());
            orderedLoops[i] = tmp;
        }
    }
}

void Nektar::NekMeshUtils::FaceMesh::Smoothing ( )

private

performs node smoothing on face

Definition at line 210 of file FaceMesh.cpp.

References ASSERTL0, Nektar::LibUtilities::eQuadrilateral, and Nektar::iterator.

{
    EdgeSet::iterator eit;
    NodeSet::iterator nit;

    map<int, vector<EdgeSharedPtr> > connectingedges;

    map<int, vector<ElementSharedPtr> > connectingelements;

    for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
    {
        connectingedges[(*eit)->m_n1->m_id].push_back(*eit);
        connectingedges[(*eit)->m_n2->m_id].push_back(*eit);
    }

    for (int i = 0; i < m_localElements.size(); i++)
    {
        vector<NodeSharedPtr> v = m_localElements[i]->GetVertexList();
        for (int j = 0; j < 3; j++)
        {
            connectingelements[v[j]->m_id].push_back(m_localElements[i]);
        }
    }

    // perform 4 runs of elastic relaxation based on the octree
    for (int q = 0; q < 4; q++)
    {
        for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
        {
            if ((*nit)->GetNumCadCurve() > 0) // node is on curve so skip
                continue;

            vector<NodeSharedPtr> connodes; // this can be real nodes or dummy
                                            // nodes depending on the system

            vector<EdgeSharedPtr> edges  = connectingedges[(*nit)->m_id];
            vector<ElementSharedPtr> els = connectingelements[(*nit)->m_id];

            bool perfrom = true;
            for (int i = 0; i < els.size(); i++)
            {
                if (els[i]->GetConf().m_e == LibUtilities::eQuadrilateral)
                {
                    perfrom = false;
                    break;
                }
            }

            if (!perfrom)
                continue;

            vector<NodeSharedPtr> nodesystem;
            vector<NekDouble> lamp;

            for (int i = 0; i < edges.size(); i++)
            {
                vector<NekDouble> lambda;

                NodeSharedPtr J;
                if (*nit == edges[i]->m_n1)
                    J = edges[i]->m_n2;
                else if (*nit == edges[i]->m_n2)
                    J = edges[i]->m_n1;
                else
                    ASSERTL0(false, "could not find node");

                Array<OneD, NekDouble> ui = (*nit)->GetCADSurfInfo(m_id);
                Array<OneD, NekDouble> uj = J->GetCADSurfInfo(m_id);

                for (int j = 0; j < els.size(); j++)
                {
                    vector<NodeSharedPtr> v = els[j]->GetVertexList();
                    if (v[0] == J || v[1] == J || v[2] == J)
                        continue; // elememt is adjacent to J therefore no
                                  // intersection on IJ

                    // need to find other edge
                    EdgeSharedPtr AtoB;
                    bool found               = false;
                    vector<EdgeSharedPtr> es = els[j]->GetEdgeList();
                    for (int k = 0; k < es.size(); k++)
                    {
                        if (!(es[k]->m_n1 == *nit || es[k]->m_n2 == *nit))
                        {
                            found = true;
                            AtoB  = es[k];
                            break;
                        }
                    }
                    ASSERTL0(found, "failed to find edge to test");

                    Array<OneD, NekDouble> A = AtoB->m_n1->GetCADSurfInfo(m_id);
                    Array<OneD, NekDouble> B = AtoB->m_n2->GetCADSurfInfo(m_id);

                    NekDouble lam = ((A[0] - uj[0]) * (B[1] - A[1]) -
                                     (A[1] - uj[1]) * (B[0] - A[0])) /
                                    ((ui[0] - uj[0]) * (B[1] - A[1]) -
                                     (ui[1] - uj[1]) * (B[0] - A[0]));

                    if (!(lam < 0) && !(lam > 1))
                        lambda.push_back(lam);
                }

                if (lambda.size() > 0)
                {
                    sort(lambda.begin(), lambda.end());
                    // make a new dummy node based on the system
                    Array<OneD, NekDouble> ud(2);
                    ud[0]                       = uj[0] + lambda[0] * (ui[0] - uj[0]);
                    ud[1]                       = uj[1] + lambda[0] * (ui[1] - uj[1]);
                    Array<OneD, NekDouble> locd = m_cadsurf->P(ud);
                    NodeSharedPtr dn = boost::shared_ptr<Node>(
                        new Node(0, locd[0], locd[1], locd[2]));
                    dn->SetCADSurf(m_id, m_cadsurf, ud);

                    nodesystem.push_back(dn);
                    lamp.push_back(lambda[0]);
                }
                else
                {
                    nodesystem.push_back(J);
                    lamp.push_back(1.0);
                }
            }

            Array<OneD, NekDouble> ui(2);
            ui[0] = 0.0;
            ui[1] = 0.0;

            for (int i = 0; i < nodesystem.size(); i++)
            {
                Array<OneD, NekDouble> uj = nodesystem[i]->GetCADSurfInfo(m_id);
                ui[0] += uj[0] / nodesystem.size();
                ui[1] += uj[1] / nodesystem.size();
            }

            Array<OneD, NekDouble> bounds = m_cadsurf->GetBounds();

            Array<OneD, NekDouble> uvn(2);

            uvn[0] = ui[0];
            uvn[1] = ui[1];

            if (!(uvn[0] < bounds[0] || uvn[0] > bounds[1] ||
                  uvn[1] < bounds[2] || uvn[1] > bounds[3]))
            {
                Array<OneD, NekDouble> l2 = m_cadsurf->P(uvn);
                (*nit)->Move(l2, m_id, uvn);
            }
        }
    }
}

void Nektar::NekMeshUtils::FaceMesh::Stretching ( )

private

Calculate the paramter plane streching factor.

Definition at line 873 of file FaceMesh.cpp.

{
    // define a sampling and calculate the aspect ratio of the paramter plane
    m_str                       = 0.0;
    Array<OneD, NekDouble> bnds = m_cadsurf->GetBounds();

    NekDouble dxu = int(bnds[1] - bnds[0] < bnds[3] - bnds[2]
                            ? 40
                            : (bnds[1] - bnds[0]) / (bnds[3] - bnds[2]) * 40);
    NekDouble dxv = int(bnds[3] - bnds[2] < bnds[1] - bnds[0]
                            ? 40
                            : (bnds[3] - bnds[2]) / (bnds[1] - bnds[0]) * 40);

    NekDouble du = (bnds[1] - bnds[0]) / dxu;
    NekDouble dv = (bnds[3] - bnds[2]) / dxv;

    int ct = 0;

    for (int i = 0; i < dxu; i++)
    {
        for (int j = 0; j < dxv; j++)
        {
            Array<OneD, NekDouble> uv(2);
            uv[0] = bnds[0] + i * du;
            uv[1] = bnds[2] + j * dv;
            if (i == dxu - 1)
                uv[0] = bnds[1];
            if (j == dxv - 1)
                uv[1]                = bnds[3];
            Array<OneD, NekDouble> r = m_cadsurf->D1(uv);

            NekDouble ru = sqrt(r[3] * r[3] + r[4] * r[4] + r[5] * r[5]);
            NekDouble rv = sqrt(r[6] * r[6] + r[7] * r[7] + r[8] * r[8]);

            ru *= du;
            rv *= dv;

            if (rv < 1E-8)
                continue;

            m_str += ru / rv;
            ct++;
        }
    }

    m_str /= ct;
}

bool Nektar::NekMeshUtils::FaceMesh::Validate ( )

private

Validate the surface mesh base on the octree and real dimensions of the edges.

Definition at line 921 of file FaceMesh.cpp.

{
    // check all edges in the current mesh for length against the octree
    // if the octree is not conformed to add a new point inside the triangle
    // if no new points are added meshing can stop
    int pointBefore = m_stienerpoints.size();
    for (int i = 0; i < m_connec.size(); i++)
    {
        Array<OneD, NekDouble> triDelta(3);

        Array<OneD, NekDouble> r(3);

        r[0] = m_connec[i][0]->Distance(m_connec[i][1]);
        r[1] = m_connec[i][1]->Distance(m_connec[i][2]);
        r[2] = m_connec[i][2]->Distance(m_connec[i][0]);

        triDelta[0] = m_octree->Query(m_connec[i][0]->GetLoc());
        triDelta[1] = m_octree->Query(m_connec[i][1]->GetLoc());
        triDelta[2] = m_octree->Query(m_connec[i][2]->GetLoc());

        int numValid = 0;

        if (r[0] < triDelta[0] && r[2] < triDelta[0])
            numValid++;

        if (r[1] < triDelta[1] && r[0] < triDelta[1])
            numValid++;

        if (r[2] < triDelta[2] && r[1] < triDelta[2])
            numValid++;

        if (numValid != 3)
        {
            Array<OneD, NekDouble> ainfo, binfo, cinfo;
            ainfo = m_connec[i][0]->GetCADSurfInfo(m_id);
            binfo = m_connec[i][1]->GetCADSurfInfo(m_id);
            cinfo = m_connec[i][2]->GetCADSurfInfo(m_id);

            Array<OneD, NekDouble> uvc(2);
            uvc[0] = (ainfo[0] + binfo[0] + cinfo[0]) / 3.0;
            uvc[1] = (ainfo[1] + binfo[1] + cinfo[1]) / 3.0;
            AddNewPoint(uvc);
        }
    }

    if (m_stienerpoints.size() == pointBefore)
    {
        return false;
    }
    else
    {
        return true;
    }
}

Friends And Related Function Documentation

friend class MemoryManager< FaceMesh >

friend

Definition at line 56 of file FaceMesh.h.

Member Data Documentation

std::vector<std::pair<NodeSharedPtr, NodeSharedPtr> > Nektar::NekMeshUtils::FaceMesh::blpairs

private

list of node links between node on loop and its corresponding interior node in quads

Definition at line 177 of file FaceMesh.h.

NekDouble Nektar::NekMeshUtils::FaceMesh::m_bl

private

boundary layer thickness

Definition at line 172 of file FaceMesh.h.

CADSurfSharedPtr Nektar::NekMeshUtils::FaceMesh::m_cadsurf

private

CAD surface.

Definition at line 147 of file FaceMesh.h.

Referenced by FaceMesh().

std::vector<std::vector<NodeSharedPtr> > Nektar::NekMeshUtils::FaceMesh::m_connec

private

triangle connectiviities

Definition at line 164 of file FaceMesh.h.

std::map<int, CurveMeshSharedPtr> Nektar::NekMeshUtils::FaceMesh::m_curvemeshes

private

Map of the curve meshes which bound the surfaces.

Definition at line 151 of file FaceMesh.h.

std::vector<EdgeLoop> Nektar::NekMeshUtils::FaceMesh::m_edgeloops

private

data structure containing the edges, their order and oreientation for the surface

Definition at line 154 of file FaceMesh.h.

Referenced by FaceMesh().

int Nektar::NekMeshUtils::FaceMesh::m_id

private

id of the surface mesh

Definition at line 156 of file FaceMesh.h.

EdgeSet Nektar::NekMeshUtils::FaceMesh::m_localEdges

private

local set of edges

Definition at line 168 of file FaceMesh.h.

std::vector<ElementSharedPtr> Nektar::NekMeshUtils::FaceMesh::m_localElements

private

local list of elements

Definition at line 170 of file FaceMesh.h.

NodeSet Nektar::NekMeshUtils::FaceMesh::m_localNodes

private

local set of nodes

Definition at line 166 of file FaceMesh.h.

bool Nektar::NekMeshUtils::FaceMesh::m_makebl

private

should build boundary layer

Definition at line 174 of file FaceMesh.h.

Referenced by FaceMesh().

MeshSharedPtr Nektar::NekMeshUtils::FaceMesh::m_mesh

private

mesh pointer

Definition at line 145 of file FaceMesh.h.

OctreeSharedPtr Nektar::NekMeshUtils::FaceMesh::m_octree

private

Octree object.

Definition at line 149 of file FaceMesh.h.

std::vector<NodeSharedPtr> Nektar::NekMeshUtils::FaceMesh::m_stienerpoints

private

list of stiener points in the triangulation

Definition at line 160 of file FaceMesh.h.

NekDouble Nektar::NekMeshUtils::FaceMesh::m_str

private

pplane stretching

Definition at line 162 of file FaceMesh.h.

std::vector<std::vector<NodeSharedPtr> > Nektar::NekMeshUtils::FaceMesh::orderedLoops

private

list of boundary nodes in their order loops

Definition at line 158 of file FaceMesh.h.