Nektar::NekMeshUtils::FaceMesh Class Reference

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

#include <FaceMesh.h>

Public Member Functions
	FaceMesh (const int id, MeshSharedPtr m, const std::map< int, CurveMeshSharedPtr > &cmeshes, const int comp)
	Default constructor. More...
void	Mesh ()
	mesh exectuation command More...
bool	ValidateCurves ()
	validate the curve meshes More...
void	ValidateLoops ()
	validate the curve meshes considering the loops More...

Private Member Functions
void	OrientateCurves ()
	Get the boundries of the surface and extracts the nodes from the curve meshes in the correct order. More...
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	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...
std::map< int, CurveMeshSharedPtr >	m_curvemeshes
	Map of the curve meshes which bound the surfaces. More...
std::vector< EdgeLoopSharedPtr >	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...
NodeSet	m_inBoundary
	set of nodes which are in the boundary (easier to identify conflicts with) More...
int	m_compId
	identity to put into element tags More...

Friends
class	MemoryManager< FaceMesh >

Detailed Description

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

Definition at line 51 of file FaceMesh.h.

Constructor & Destructor Documentation

FaceMesh()

Nektar::NekMeshUtils::FaceMesh::FaceMesh ( const int  id, MeshSharedPtr  m, const std::map< int, CurveMeshSharedPtr > &  cmeshes, const int  comp  )

inline

Default constructor.

Definition at line 59 of file FaceMesh.h.

        : m_mesh(m), m_curvemeshes(cmeshes), m_id(id), m_compId(comp)
 
    {
        m_cadsurf = m_mesh->m_cad->GetSurf(m_id);
        m_edgeloops = m_cadsurf->GetEdges();
    };

References m_cadsurf, m_edgeloops, m_id, and m_mesh.

Member Function Documentation

AddNewPoint()

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

private

adds a new stiener point to the triangulation for meshing

Definition at line 1129 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_mesh->m_octree->Query(np);
 
    NodeSharedPtr n = std::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_cadsurf, uv);
        m_stienerpoints.push_back(n);
    }
}

BuildLocalMesh()

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

private

build a local version of mesh elements

Definition at line 902 of file FaceMesh.cpp.

{
    /*************************
    // build a local set of nodes edges and elemenets for optimstaion prior to
    putting them into m_mesh
    */
 
    for (int i = 0; i < m_connec.size(); i++)
    {
        ElmtConfig conf(LibUtilities::eTriangle, 1, false, false,
                        m_mesh->m_spaceDim != 3);
 
        vector<int> tags;
        tags.push_back(m_compId);
        ElementSharedPtr E = GetElementFactory().CreateInstance(
            LibUtilities::eTriangle, conf, m_connec[i], tags);
        E->m_parentCAD = m_cadsurf;
 
        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]);
        }
 
        E->SetId(m_localElements.size());
        m_localElements.push_back(E);
    }
 
    for (int i = 0; i < m_localElements.size(); ++i)
    {
        for (int j = 0; j < m_localElements[i]->GetEdgeCount(); ++j)
        {
            pair<EdgeSet::iterator, bool> testIns;
            EdgeSharedPtr ed = m_localElements[i]->GetEdge(j);
            // look for edge in m_mesh edgeset from curves
            EdgeSet::iterator s = m_mesh->m_edgeSet.find(ed);
            if (!(s == m_mesh->m_edgeSet.end()))
            {
                ed = *s;
                m_localElements[i]->SetEdge(j, *s);
            }
 
            testIns = m_localEdges.insert(ed);
 
            if (testIns.second)
            {
                EdgeSharedPtr ed2 = *testIns.first;
                ed2->m_elLink.push_back(
                    pair<ElementSharedPtr, int>(m_localElements[i], j));
            }
            else
            {
                EdgeSharedPtr e2 = *(testIns.first);
                m_localElements[i]->SetEdge(j, e2);
                e2->m_elLink.push_back(
                    pair<ElementSharedPtr, int>(m_localElements[i], j));
            }
        }
    }
}

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

DiagonalSwap()

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

private

performs diagonal swapping of edges

Definition at line 506 of file FaceMesh.cpp.

{
    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++)
    {
        // this routine is broken and needs looking at
        NodeSet::iterator f = m_inBoundary.find((*nit));
        if (f == m_inBoundary.end()) // node isnt on curve so skip
        {
            // 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;*/
            idealConnec[(*nit)->m_id] = 4;
        }
    }
    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;
 
            NodeSet::iterator f1 = m_inBoundary.find((*it)->m_n1);
            NodeSet::iterator f2 = m_inBoundary.find((*it)->m_n2);
            if (f1 != m_inBoundary.end() && f2 != m_inBoundary.end())
            {
                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
            // //cout<< A->GetNumCADSurf() << " " << B->GetNumCADSurf() << " "
            // //    << C->GetNumCADSurf() << " " << D->GetNumCADSurf() << endl;
            // ofstream file;
            // file.open("pts.3D");
            // file << "x y z value" << endl;
            // file << A->m_x << " " << A->m_y << " " << A->m_z << endl;
            // file << B->m_x << " " << B->m_y << " " << B->m_z << endl;
            // file << C->m_x << " " << C->m_y << " " << C->m_z << endl;
            // file << D->m_x << " " << D->m_y << " " << D->m_z << endl;
            // file.close();
            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 = std::shared_ptr<Edge>(new Edge(C, A));
                ADt = std::shared_ptr<Edge>(new Edge(A, D));
                DBt = std::shared_ptr<Edge>(new Edge(D, B));
                BCt = std::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 = std::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,
                                m_mesh->m_spaceDim != 3);
                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->m_parentCAD = m_cadsurf;
                ntri2->m_parentCAD = m_cadsurf;
 
                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");
    }
}

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

MakeBL()

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

private

adds a quad layer around any interior loops

Mesh()

void Nektar::NekMeshUtils::FaceMesh::Mesh

(

)

mesh exectuation command

Definition at line 145 of file FaceMesh.cpp.

{
    Stretching();
    OrientateCurves();
 
    int numPoints = 0;
    for (int i = 0; i < orderedLoops.size(); i++)
    {
        numPoints += orderedLoops[i].size();
        for (int j = 0; j < orderedLoops[i].size(); j++)
        {
            m_inBoundary.insert(orderedLoops[i][j]);
        }
    }
 
    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, "number of verts in face is less than 3");
 
    // 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;
 
    // continuously remesh until all triangles conform to the spacing in the
    // octree
    while (repeat)
    {
        repeat = Validate();
        if (!repeat)
        {
            break;
        }
        m_connec.clear();
        pplanemesh->AssignStiener(m_stienerpoints);
        pplanemesh->Mesh();
        pplanemesh->Extract(m_connec);
        meshcounter++;
        //    break;
    }
 
    // build a local version of the mesh (one set of triangles).  this is done
    // so edge connectivity infomration can be used for optimisation
    BuildLocalMesh();
 
    OptimiseLocalMesh();
 
    // make new elements and add to list from list of nodes and connectivity
    // from triangle removing unnesercary infomration from the elements
    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[2].push_back(m_localElements[i]);
    }
 
    if (m_mesh->m_verbose)
    {
        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;
    }
}

References ASSERTL0.

OptimiseLocalMesh()

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

private

function which calls the optimisation routines

Definition at line 243 of file FaceMesh.cpp.

{
    // each optimisation algorithm is based on the work in chapter 19
    DiagonalSwap();
 
    Smoothing();
 
    DiagonalSwap();
 
    Smoothing();
}

OrientateCurves()

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 1175 of file FaceMesh.cpp.

{
    // this could be a second run on orentate so clear some info
    orderedLoops.clear();
 
    // 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] == CADOrientation::eForwards)
            {
                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);
    }
}

References Nektar::StdRegions::eForwards.

Smoothing()

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

private

performs node smoothing on face

Definition at line 255 of file FaceMesh.cpp.

{
    EdgeSet::iterator eit;
    NodeSet::iterator nit;
 
    Array<OneD, NekDouble> bounds = m_cadsurf->GetBounds();
 
    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++)
        {
            NodeSet::iterator f = m_inBoundary.find((*nit));
 
            // node is on curve so skip
            if (f != m_inBoundary.end())
            {
                continue;
            }
 
            // this can be real nodes or dummy nodes depending on the system
            vector<NodeSharedPtr> connodes;
 
            vector<EdgeSharedPtr> edges  = connectingedges[(*nit)->m_id];
            vector<ElementSharedPtr> els = connectingelements[(*nit)->m_id];
 
            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();
 
                    // elememt is adjacent to J therefore no intersection on IJ
                    if (v[0] == J || v[1] == J || v[2] == J)
                    {
                        continue;
                    }
 
                    // 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 = std::shared_ptr<Node>(
                        new Node(0, locd[0], locd[1], locd[2]));
                    dn->SetCADSurf(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> u0(2);
            u0[0] = 0.0;
            u0[1] = 0.0;
 
            for (int i = 0; i < nodesystem.size(); i++)
            {
                Array<OneD, NekDouble> uj = nodesystem[i]->GetCADSurfInfo(m_id);
                u0[0] += uj[0] / nodesystem.size();
                u0[1] += uj[1] / nodesystem.size();
            }
 
            /*Array<OneD, NekDouble> pu0 = m_cadsurf->P(u0);
            NekDouble di = m_mesh->m_octree->Query(pu0);
            Array<OneD, NekDouble> F(2, 0.0), dF(4, 0.0);
            for (int i = 0; i < nodesystem.size(); i++)
            {
                Array<OneD, NekDouble> rj = nodesystem[i]->GetLoc();
                Array<OneD, NekDouble> uj = nodesystem[i]->GetCADSurfInfo(m_id);
                NekDouble dj  = m_mesh->m_octree->Query(rj);
                NekDouble d   = (di + dj) / 2.0;
                NekDouble wij = sqrt((rj[0] - pu0[0]) * (rj[0] - pu0[0]) +
                                     (rj[1] - pu0[1]) * (rj[1] - pu0[1]) +
                                     (rj[2] - pu0[2]) * (rj[2] - pu0[2])) -
                                d;
 
                NekDouble umag = sqrt((uj[0] - u0[0]) * (uj[0] - u0[0]) +
                                      (uj[1] - u0[1]) * (uj[1] - u0[1]));
 
                F[0] += wij * (uj[0] - u0[0]) / umag;
                F[1] += wij * (uj[1] - u0[1]) / umag;
 
                Array<OneD, NekDouble> d1 = m_cadsurf->D1(u0);
 
                Array<OneD, NekDouble> dw(2, 0.0);
                dw[0] = -2.0 *
                        ((rj[0] - pu0[0]) * d1[3] + (rj[1] - pu0[1]) * d1[4] +
                         (rj[2] - pu0[2]) * d1[5]);
                dw[1] = -2.0 *
                        ((rj[0] - pu0[0]) * d1[6] + (rj[1] - pu0[1]) * d1[7] +
                         (rj[2] - pu0[2]) * d1[8]);
 
                Array<OneD, NekDouble> drhs(4, 0.0);
                drhs[0] = 2.0 * ((uj[0] - u0[0]) * (uj[0] - u0[0]) - umag) /
                          umag / umag;
                drhs[1] =
                    2.0 * ((uj[0] - u0[0]) * (uj[1] - u0[1])) / umag / umag;
                drhs[2] = drhs[1];
                drhs[3] = 2.0 * ((uj[1] - u0[1]) * (uj[1] - u0[1]) - umag) /
                          umag / umag;
 
                dF[0] += dw[0] * (uj[0] - u0[0]) / umag + wij * drhs[0];
 
                dF[1] += dw[0] * (uj[1] - u0[1]) / umag + wij * drhs[1];
 
                dF[2] += dw[1] * (uj[0] - u0[0]) / umag + wij * drhs[2];
 
                dF[3] += dw[1] * (uj[1] - u0[1]) / umag + wij * drhs[3];
            }
 
            NekDouble det = dF[0] * dF[3] - dF[1] * dF[2];
            NekDouble tmp = dF[0];
            dF[0]         = dF[3] / det;
            dF[3]         = tmp / det;
            dF[1] *= -1.0 / det;
            dF[2] *= -1.0 / det;
 
            u0[0] -= (dF[0] * F[0] + dF[2] * F[1]);
            u0[1] -= (dF[1] * F[0] + dF[3] * F[1]);*/
 
            bool inbounds = true;
            if (u0[0] < bounds[0])
            {
                inbounds = false;
            }
            else if (u0[0] > bounds[1])
            {
                inbounds = false;
            }
            else if (u0[1] < bounds[2])
            {
                inbounds = false;
            }
            else if (u0[1] > bounds[3])
            {
                inbounds = false;
            }
 
            if (!inbounds)
            {
                continue;
            }
 
            /*Array<OneD, NekDouble> FN(2, 0.0);
            pu0 = m_cadsurf->P(u0);
            di  = m_mesh->m_octree->Query(pu0);
            for (int i = 0; i < nodesystem.size(); i++)
            {
                Array<OneD, NekDouble> rj = nodesystem[i]->GetLoc();
                Array<OneD, NekDouble> uj = nodesystem[i]->GetCADSurfInfo(m_id);
                NekDouble dj  = m_mesh->m_octree->Query(rj);
                NekDouble d   = (di + dj) / 2.0;
                NekDouble wij = sqrt((rj[0] - pu0[0]) * (rj[0] - pu0[0]) +
                                     (rj[1] - pu0[1]) * (rj[1] - pu0[1]) +
                                     (rj[2] - pu0[2]) * (rj[2] - pu0[2])) - d;
 
                NekDouble umag = sqrt((uj[0] - u0[0]) * (uj[0] - u0[0]) +
                                      (uj[1] - u0[1]) * (uj[1] - u0[1]));
 
                FN[0] += wij * (uj[0] - u0[0]) / umag;
                FN[1] += wij * (uj[1] - u0[1]) / umag;
            }
 
            if (F[0] * F[0] + F[1] * F[1] < FN[0] * FN[0] + FN[1] * FN[1])
            {
                continue;
            }*/
 
            Array<OneD, NekDouble> l2 = m_cadsurf->P(u0);
            (*nit)->Move(l2, m_id, u0);
        }
    }
}

References ASSERTL0.

Stretching()

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

private

Calculate the paramter plane streching factor.

Definition at line 964 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;
}

Validate()

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

private

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

Definition at line 1018 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> r(3), a(3);
 
        vector<Array<OneD, NekDouble> > info;
 
        for (int j = 0; j < 3; j++)
        {
            info.push_back(m_connec[i][j]->GetCADSurfInfo(m_id));
        }
 
        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]);
 
        a[0] = m_connec[i][0]->Angle(m_connec[i][1]->GetLoc(),
                                     m_connec[i][2]->GetLoc(),
                                     m_cadsurf->N(info[0]));
        a[1] = m_connec[i][1]->Angle(m_connec[i][2]->GetLoc(),
                                     m_connec[i][0]->GetLoc(),
                                     m_cadsurf->N(info[1]));
        a[2] = m_connec[i][2]->Angle(m_connec[i][0]->GetLoc(),
                                     m_connec[i][1]->GetLoc(),
                                     m_cadsurf->N(info[2]));
 
        NekDouble d1 = m_mesh->m_octree->Query(m_connec[i][0]->GetLoc());
        NekDouble d2 = m_mesh->m_octree->Query(m_connec[i][1]->GetLoc());
        NekDouble d3 = m_mesh->m_octree->Query(m_connec[i][2]->GetLoc());
 
        Array<OneD, NekDouble> uvc(2);
        uvc[0] = (info[0][0] + info[1][0] + info[2][0]) / 3.0;
        uvc[1] = (info[0][1] + info[1][1] + info[2][1]) / 3.0;
 
        Array<OneD, NekDouble> locc = m_cadsurf->P(uvc);
        NekDouble d4 = m_mesh->m_octree->Query(locc);
 
        NekDouble d = (d1 + d2 + d3 + d4) / 4.0;
 
        vector<bool> valid(3);
        valid[0] = r[0] < d * 1.5;
        valid[1] = r[1] < d * 1.5;
        valid[2] = r[2] < d * 1.5;
 
        vector<bool> angValid(3);
        angValid[0] = a[0] / M_PI * 180.0 > 20.0 && a[0] / M_PI * 180.0 < 120.0;
        angValid[1] = a[1] / M_PI * 180.0 > 20.0 && a[1] / M_PI * 180.0 < 120.0;
        angValid[2] = a[2] / M_PI * 180.0 > 20.0 && a[2] / M_PI * 180.0 < 120.0;
 
        int numValid    = 0;
        int numAngValid = 0;
        for (int j = 0; j < 3; j++)
        {
            if (valid[j])
            {
                numValid++;
            }
            if (angValid[j])
            {
                numAngValid++;
            }
        }
 
        // if numvalid is zero no work to be done
        /*if (numValid != 3)
        {
            AddNewPoint(uvc);
        }*/
 
        if (numValid != 3 || numAngValid != 3)
        {
            // break the bad edge
            /*int a=0, b=0;
            if(!valid[0])
            {
                a = 0;
                b = 1;
            }
            else if(!valid[1])
            {
                a = 1;
                b = 2;
            }
            else if(!valid[2])
            {
                a = 2;
                b = 0;
            }
 
            Array<OneD, NekDouble> uvn(2);
            uvn[0] = (info[a][0] + info[b][0]) / 2.0;
            uvn[1] = (info[a][1] + info[b][1]) / 2.0;*/
            AddNewPoint(uvc);
        }
    }
 
    if (m_stienerpoints.size() == pointBefore)
    {
        return false;
    }
    else
    {
        return true;
    }
}

ValidateCurves()

bool Nektar::NekMeshUtils::FaceMesh::ValidateCurves

(

)

validate the curve meshes

Definition at line 46 of file FaceMesh.cpp.

{
    vector<int> curvesInSurface;
    for (int i = 0; i < m_edgeloops.size(); i++)
    {
        for (int j = 0; j < m_edgeloops[i]->edges.size(); j++)
        {
            curvesInSurface.push_back(m_edgeloops[i]->edges[j]->GetId());
        }
    }
 
    bool error = false;
 
    for (int i = 0; i < curvesInSurface.size(); i++)
    {
        vector<EdgeSharedPtr> es =
            m_curvemeshes[curvesInSurface[i]]->GetMeshEdges();
 
        for (int j = i; j < curvesInSurface.size(); j++)
        {
            if (i == j)
            {
                continue;
            }
 
            vector<EdgeSharedPtr> es2 =
                m_curvemeshes[curvesInSurface[j]]->GetMeshEdges();
 
            for (int l = 0; l < es.size(); l++)
            {
                Array<OneD, NekDouble> P1 = es[l]->m_n1->GetCADSurfInfo(m_id);
                Array<OneD, NekDouble> P2 = es[l]->m_n2->GetCADSurfInfo(m_id);
                for (int k = 0; k < es2.size(); k++)
                {
                    if (es[l]->m_n1 == es2[k]->m_n1 ||
                        es[l]->m_n1 == es2[k]->m_n2 ||
                        es[l]->m_n2 == es2[k]->m_n1 ||
                        es[l]->m_n2 == es2[k]->m_n2)
                    {
                        continue;
                    }
 
                    Array<OneD, NekDouble> P3 =
                        es2[k]->m_n1->GetCADSurfInfo(m_id);
                    Array<OneD, NekDouble> P4 =
                        es2[k]->m_n2->GetCADSurfInfo(m_id);
 
                    NekDouble den = (P4[0] - P3[0]) * (P2[1] - P1[1]) -
                                    (P2[0] - P1[0]) * (P4[1] - P3[1]);
 
                    if (fabs(den) < 1e-8)
                    {
                        continue;
                    }
 
                    NekDouble t = ((P1[0] - P3[0]) * (P4[1] - P3[1]) -
                                   (P4[0] - P3[0]) * (P1[1] - P3[1])) /
                                  den;
                    NekDouble u =
                        (P1[0] - P3[0] + t * (P2[0] - P1[0])) / (P4[0] - P3[0]);
 
                    if (t < 1.0 && t > 0.0 && u < 1.0 && u > 0.0)
                    {
                        Array<OneD, NekDouble> uv(2);
                        uv[0] = P1[0] + t * (P2[0] - P1[0]);
                        uv[1] = P1[1] + t * (P2[1] - P1[1]);
                        Array<OneD, NekDouble> loc = m_cadsurf->P(uv);
                        cout << endl
                             << "Curve mesh error at " << loc[0] << " "
                             << loc[1] << " " << loc[2] << " on face " << m_id
                             << endl;
                        error = true;
                    }
                }
            }
        }
    }
    return error;
}

References CG_Iterations::loc.

ValidateLoops()

void Nektar::NekMeshUtils::FaceMesh::ValidateLoops

(

)

validate the curve meshes considering the loops

Definition at line 126 of file FaceMesh.cpp.

{
    OrientateCurves();
 
    for (int i = 0; i < orderedLoops.size(); i++)
    {
        int numPoints = orderedLoops[i].size();
        if(numPoints == 2)
        {
            //force a remesh of the curves
            for(int j = 0; j < m_edgeloops[i]->edges.size(); j++)
            {
                int cid = m_edgeloops[i]->edges[j]->GetId();
                m_curvemeshes[cid]->ReMesh();
            }
        }
    }
}

Friends And Related Function Documentation

MemoryManager< FaceMesh >

friend class MemoryManager< FaceMesh >

friend

Definition at line 165 of file FaceMesh.h.

Member Data Documentation

m_cadsurf

CADSurfSharedPtr Nektar::NekMeshUtils::FaceMesh::m_cadsurf

private

CAD surface.

Definition at line 137 of file FaceMesh.h.

Referenced by FaceMesh().

m_compId

int Nektar::NekMeshUtils::FaceMesh::m_compId

private

identity to put into element tags

Definition at line 162 of file FaceMesh.h.

m_connec

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

private

triangle connectiviities

Definition at line 152 of file FaceMesh.h.

m_curvemeshes

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

private

Map of the curve meshes which bound the surfaces.

Definition at line 139 of file FaceMesh.h.

m_edgeloops

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

private

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

Definition at line 142 of file FaceMesh.h.

Referenced by FaceMesh().

m_id

int Nektar::NekMeshUtils::FaceMesh::m_id

private

id of the surface mesh

Definition at line 144 of file FaceMesh.h.

Referenced by FaceMesh().

m_inBoundary

NodeSet Nektar::NekMeshUtils::FaceMesh::m_inBoundary

private

set of nodes which are in the boundary (easier to identify conflicts with)

Definition at line 160 of file FaceMesh.h.

m_localEdges

EdgeSet Nektar::NekMeshUtils::FaceMesh::m_localEdges

private

local set of edges

Definition at line 156 of file FaceMesh.h.

m_localElements

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

private

local list of elements

Definition at line 158 of file FaceMesh.h.

m_localNodes

NodeSet Nektar::NekMeshUtils::FaceMesh::m_localNodes

private

local set of nodes

Definition at line 154 of file FaceMesh.h.

m_mesh

MeshSharedPtr Nektar::NekMeshUtils::FaceMesh::m_mesh

private

mesh pointer

Definition at line 135 of file FaceMesh.h.

Referenced by FaceMesh().

m_stienerpoints

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

private

list of stiener points in the triangulation

Definition at line 148 of file FaceMesh.h.

m_str

NekDouble Nektar::NekMeshUtils::FaceMesh::m_str

private

pplane stretching

Definition at line 150 of file FaceMesh.h.

orderedLoops

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

private

list of boundary nodes in their order loops

Definition at line 146 of file FaceMesh.h.