Nektar::NekMeshUtils::SurfaceMesh Class Reference

class containing all surface meshing routines methods and classes More...

#include <SurfaceMesh.h>

Collaboration diagram for Nektar::NekMeshUtils::SurfaceMesh:

Public Member Functions
	SurfaceMesh (MeshSharedPtr m, CADSystemSharedPtr cad, OctreeSharedPtr octree, std::vector< unsigned int > sy, NekDouble b)
	Default constructor, requires the cad and octree objects to begin. More...
void	Mesh ()
	Run all linear meshing routines. More...
void	HOSurf ()
	run all high-order surface meshing routines More...
void	Validate ()
	Validate the linear surface mesh. More...
void	Report ()
	Print brief information to screen. More...

Private Attributes
MeshSharedPtr	m_mesh
	mesh object More...
CADSystemSharedPtr	m_cad
	CAD object. More...
OctreeSharedPtr	m_octree
	Octree object. More...
std::map< int, FaceMeshSharedPtr >	m_facemeshes
	map of individual surface meshes from parametric surfaces More...
std::map< int, CurveMeshSharedPtr >	m_curvemeshes
	map of individual curve meshes of the curves in the domain More...
std::vector< unsigned int >	m_symsurfs
	list of sysmetry plane surfaces to build quads onto More...
NekDouble	m_bl
	thickness of the boundary layer if needed More...

Friends
class	MemoryManager< SurfaceMesh >

Detailed Description

class containing all surface meshing routines methods and classes

Definition at line 53 of file SurfaceMesh.h.

Constructor & Destructor Documentation

Nektar::NekMeshUtils::SurfaceMesh::SurfaceMesh ( MeshSharedPtr  m, CADSystemSharedPtr  cad, OctreeSharedPtr  octree, std::vector< unsigned int >  sy, NekDouble  b  )

inline

Default constructor, requires the cad and octree objects to begin.

Definition at line 62 of file SurfaceMesh.h.

        : m_mesh(m), m_cad(cad), m_octree(octree), m_symsurfs(sy), m_bl(b){};

Member Function Documentation

void Nektar::NekMeshUtils::SurfaceMesh::HOSurf

(

)

run all high-order surface meshing routines

Definition at line 241 of file SurfaceMeshHOMesh.cpp.

References ASSERTL0, Nektar::NekMeshUtils::BGFSUpdate(), Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eNodalTriFekete, Nektar::LibUtilities::eQuadrilateral, Nektar::iterator, Nektar::NekMeshUtils::ListOfFaceSpings(), Nektar::LibUtilities::PointsManager(), Nektar::LibUtilities::PrintProgressbar(), Nektar::NekMeshUtils::surf, and Nektar::NekMeshUtils::weights().

{
    if (m_mesh->m_verbose)
        cout << endl << "\tHigh-Order Surface meshing" << endl;

    // this bit of code sets up information for the standard edge and face.
    // and a mapping for node ordering for spring optimistaion

    LibUtilities::PointsKey ekey(m_mesh->m_nummode,
                                 LibUtilities::eGaussLobattoLegendre);
    Array<OneD, NekDouble> gll;

    LibUtilities::PointsManager()[ekey]->GetPoints(gll);

    LibUtilities::PointsKey pkey(m_mesh->m_nummode,
                                 LibUtilities::eNodalTriFekete);
    Array<OneD, NekDouble> u, v;

    int nq = m_mesh->m_nummode;

    LibUtilities::PointsManager()[pkey]->GetPoints(u, v);

    set<pair<int, int> > springs     = ListOfFaceSpings(nq);
    map<pair<int, int>, NekDouble> z = weights(springs, u, v);

    // because edges are listed twice need a method to not repeat over them
    EdgeSet completedEdges;

    // loop over all the faces in the surface mesh, check all three edges for
    // high order info, if nothing high-order the edge.

    for (int i = 0; i < m_mesh->m_element[2].size(); i++)
    {
        if (m_mesh->m_verbose)
        {
            LibUtilities::PrintProgressbar(
                i, m_mesh->m_element[2].size(), "\t\tSurface elements");
        }

        if (m_mesh->m_element[2][i]->GetConf().m_e ==
            LibUtilities::eQuadrilateral)
        {
            // not setup for high-order quads yet
            continue;
        }

        int surf           = m_mesh->m_element[2][i]->CADSurfId;
        CADSurfSharedPtr s = m_cad->GetSurf(surf);

        FaceSharedPtr f = m_mesh->m_element[2][i]->GetFaceLink();
        vector<EdgeSharedPtr> surfedges =
            m_mesh->m_element[2][i]->GetEdgeList();

        vector<EdgeSharedPtr> edges = f->m_edgeList;
        for (int j = 0; j < edges.size(); j++)
        {
            // test insert the edge into completedEdges
            // if the edge already exists move on
            // if not figure out its high-order information

            EdgeSet::iterator test = completedEdges.find(edges[j]);

            if (test != completedEdges.end())
            {
                continue;
            }

            EdgeSharedPtr e = edges[j];

            // the edges in the element are different to those in the face
            // the cad information is stored in the element edges which are not
            // in the m_mesh->m_edgeSet group.
            // need to link them together and copy the cad information to be
            // able to identify how to make it high-order
            bool foundsurfaceedge = false;
            for (int k = 0; k < surfedges.size(); k++)
            {
                if (surfedges[k] == e)
                {
                    e->onCurve       = surfedges[k]->onCurve;
                    e->CADCurveId    = surfedges[k]->CADCurveId;
                    e->CADCurve      = surfedges[k]->CADCurve;
                    foundsurfaceedge = true;
                }
            }
            ASSERTL0(foundsurfaceedge,
                     "cannot find corresponding surface edge");

            if (e->onCurve)
            {
                int cid             = e->CADCurveId;
                CADCurveSharedPtr c = e->CADCurve;
                NekDouble tb        = e->m_n1->GetCADCurveInfo(cid);
                NekDouble te        = e->m_n2->GetCADCurveInfo(cid);

                // distrobute points along curve as inital guess
                Array<OneD, NekDouble> ti(m_mesh->m_nummode);
                for (int k = 0; k < m_mesh->m_nummode; k++)
                {
                    ti[k] =
                        tb * (1.0 - gll[k]) / 2.0 + te * (1.0 + gll[k]) / 2.0;
                }

                Array<OneD, NekDouble> xi(nq - 2);
                for (int k = 1; k < nq - 1; k++)
                {
                    xi[k - 1] = ti[k];
                }

                OptiEdgeSharedPtr opti =
                    MemoryManager<OptiEdge>::AllocateSharedPtr(ti, gll, c);

                DNekMat B(
                    nq - 2, nq - 2, 0.0); // approximate hessian (I to start)
                for (int k = 0; k < nq - 2; k++)
                {
                    B(k, k) = 1.0;
                }
                DNekMat H(nq - 2,
                          nq - 2,
                          0.0); // approximate inverse hessian (I to start)
                for (int k = 0; k < nq - 2; k++)
                {
                    H(k, k) = 1.0;
                }

                DNekMat J = opti->dF(xi);

                Array<OneD, NekDouble> bnds = c->Bounds();

                bool repeat = true;
                int itct = 0;
                while (repeat)
                {
                    NekDouble Norm = 0;
                    for (int k = 0; k < nq - 2; k++)
                    {
                        Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
                                (bnds[1] - bnds[0]);
                    }
                    Norm = sqrt(Norm);

                    if (Norm < 1E-5)
                    {
                        repeat = false;
                        break;
                    }
                    if (itct > 100)
                    {
                        cout << "failed to optimise on curve" << endl;
                        for (int k = 0; k < nq; k++)
                        {
                            ti[k] = tb * (1.0 - gll[k]) / 2.0 +
                                    te * (1.0 + gll[k]) / 2.0;
                        }
                        exit(-1);
                        break;
                    }
                    itct++;

                    if (!BGFSUpdate(opti, J, B, H))
                    {
                        cout << "BFGS reported no update, curve on "
                             << c->GetId() << endl;
                        break;
                    }
                }
                // need to pull the solution out of opti
                ti = opti->GetSolution();

                vector<CADSurfSharedPtr> s = c->GetAdjSurf();

                ASSERTL0(s.size() == 2, "Number of common surfs should be 2");

                vector<NodeSharedPtr> honodes(m_mesh->m_nummode - 2);
                for (int k = 1; k < m_mesh->m_nummode - 1; k++)
                {
                    Array<OneD, NekDouble> loc = c->P(ti[k]);
                    NodeSharedPtr nn = boost::shared_ptr<Node>(
                        new Node(0, loc[0], loc[1], loc[2]));

                    nn->SetCADCurve(cid, c, ti[k]);
                    Array<OneD, NekDouble> uv = s[0]->locuv(loc);
                    nn->SetCADSurf(s[0]->GetId(), s[0], uv);
                    uv = s[1]->locuv(loc);
                    nn->SetCADSurf(s[1]->GetId(), s[1], uv);
                    honodes[k - 1] = nn;
                }

                e->m_edgeNodes = honodes;
                e->m_curveType = LibUtilities::eGaussLobattoLegendre;
                completedEdges.insert(e);
            }
            else
            {
                // edge is on surface and needs 2d optimisation
                CADSurfSharedPtr s = m_cad->GetSurf(surf);
                Array<OneD, NekDouble> uvb, uve;
                uvb = e->m_n1->GetCADSurfInfo(surf);
                uve = e->m_n2->GetCADSurfInfo(surf);
                Array<OneD, Array<OneD, NekDouble> > uvi(nq);
                for (int k = 0; k < nq; k++)
                {
                    Array<OneD, NekDouble> uv(2);
                    uv[0] = uvb[0] * (1.0 - gll[k]) / 2.0 +
                            uve[0] * (1.0 + gll[k]) / 2.0;
                    uv[1] = uvb[1] * (1.0 - gll[k]) / 2.0 +
                            uve[1] * (1.0 + gll[k]) / 2.0;
                    uvi[k] = uv;
                }

                Array<OneD, NekDouble> bnds = s->GetBounds();
                Array<OneD, NekDouble> all(2 * nq);
                for (int k = 0; k < nq; k++)
                {
                    all[k * 2 + 0] = uvi[k][0];
                    all[k * 2 + 1] = uvi[k][1];
                }

                Array<OneD, NekDouble> xi(2 * (nq - 2));
                for (int k = 1; k < nq - 1; k++)
                {
                    xi[(k - 1) * 2 + 0] = all[k * 2 + 0];
                    xi[(k - 1) * 2 + 1] = all[k * 2 + 1];
                }

                OptiEdgeSharedPtr opti =
                    MemoryManager<OptiEdge>::AllocateSharedPtr(all, gll, s);

                DNekMat B(2 * (nq - 2),
                          2 * (nq - 2),
                          0.0); // approximate hessian (I to start)
                for (int k = 0; k < 2 * (nq - 2); k++)
                {
                    B(k, k) = 1.0;
                }
                DNekMat H(2 * (nq - 2),
                          2 * (nq - 2),
                          0.0); // approximate inverse hessian (I to start)
                for (int k = 0; k < 2 * (nq - 2); k++)
                {
                    H(k, k) = 1.0;
                }

                DNekMat J = opti->dF(xi);

                bool repeat = true;
                int itct    = 0;
                while (repeat)
                {
                    NekDouble Norm = 0;
                    for (int k = 0; k < 2 * (nq - 2); k++)
                    {
                        if (k % 2 == 0)
                        {
                            Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
                                    (bnds[1] - bnds[0]);
                        }
                        else
                        {
                            Norm += J(k, 0) * J(k, 0) / (bnds[3] - bnds[2]) /
                                    (bnds[3] - bnds[2]);
                        }
                    }
                    Norm = sqrt(Norm);

                    if (Norm < 1E-5)
                    {
                        repeat = false;
                        break;
                    }

                    if (itct > 100)
                    {
                        cout << "failed to optimise on edge" << endl;
                        for (int k = 0; k < nq; k++)
                        {
                            Array<OneD, NekDouble> uv(2);
                            uv[0] = uvb[0] * (1.0 - gll[k]) / 2.0 +
                                    uve[0] * (1.0 + gll[k]) / 2.0;
                            uv[1] = uvb[1] * (1.0 - gll[k]) / 2.0 +
                                    uve[1] * (1.0 + gll[k]) / 2.0;
                            uvi[k] = uv;
                        }
                        exit(-1);
                        break;
                    }
                    itct++;

                    if (!BGFSUpdate(opti, J, B, H))
                    {
                        cout << "BFGS reported no update, edge on " << surf
                             << endl;
                        // exit(-1);
                        break;
                    }
                }

                all = opti->GetSolution();

                // need to put all backinto uv
                for (int k = 0; k < nq; k++)
                {
                    uvi[k][0] = all[k * 2 + 0];
                    uvi[k][1] = all[k * 2 + 1];
                }

                vector<NodeSharedPtr> honodes(nq - 2);
                for (int k = 1; k < nq - 1; k++)
                {
                    Array<OneD, NekDouble> loc;
                    loc              = s->P(uvi[k]);
                    NodeSharedPtr nn = boost::shared_ptr<Node>(
                        new Node(0, loc[0], loc[1], loc[2]));
                    nn->SetCADSurf(s->GetId(), s, uvi[k]);
                    honodes[k - 1] = nn;
                }

                e->m_edgeNodes = honodes;
                e->m_curveType = LibUtilities::eGaussLobattoLegendre;
                completedEdges.insert(e);
            }
        }

        ASSERTL0(nq <= 5, "not setup for high-orders yet");

        /*vector<NodeSharedPtr> vertices = f->m_vertexList;

        SpatialDomains::Geometry2DSharedPtr geom = f->GetGeom(3);
        geom->FillGeom();
        StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
        Array<OneD, NekDouble> coeffs0 = geom->GetCoeffs(0);
        Array<OneD, NekDouble> coeffs1 = geom->GetCoeffs(1);
        Array<OneD, NekDouble> coeffs2 = geom->GetCoeffs(2);

        Array<OneD, NekDouble> xc(nq*nq);
        Array<OneD, NekDouble> yc(nq*nq);
        Array<OneD, NekDouble> zc(nq*nq);

        xmap->BwdTrans(coeffs0,xc);
        xmap->BwdTrans(coeffs1,yc);
        xmap->BwdTrans(coeffs2,zc);


        //build an array of all uvs
        Array<OneD, Array<OneD, NekDouble> > uvi(np);
        int ctr = 0;
        for(int j = 0; j < vertices.size(); j++)
        {
            uvi[ctr++] = vertices[j]->GetCADSurfInfo(surf);
        }
        for(int j = 0; j < edges.size(); j++)
        {
            vector<NodeSharedPtr> ns = edges[j]->m_edgeNodes;
            if(!(edges[j]->m_n1 == vertices[j]))
            {
                reverse(ns.begin(),ns.end());
            }
            for(int k = 0; k < ns.size(); k++)
            {
                uvi[ctr++] = ns[k]->GetCADSurfInfo(surf);
            }
        }
        for(int j = np-ni; j < np; j++)
        {
            Array<OneD, NekDouble> xp(2);
            xp[0] = u[j];
            xp[1] = v[j];

            Array<OneD, NekDouble> xyz(3);
            xyz[0] = xmap->PhysEvaluate(xp, xc);
            xyz[1] = xmap->PhysEvaluate(xp, yc);
            xyz[2] = xmap->PhysEvaluate(xp, zc);

            Array<OneD, NekDouble> uv(2);
            s->ProjectTo(xyz,uv);
            uvi[ctr++] = uv;
        }

        OptiFaceSharedPtr opti = MemoryManager<OptiFace>::
                                    AllocateSharedPtr(uvi, z, springs, s);

        DNekMat B(2*ni,2*ni,0.0); //approximate hessian (I to start)
        for(int k = 0; k < 2*ni; k++)
        {
            B(k,k) = 1.0;
        }
        DNekMat H(2*ni,2*ni,0.0); //approximate inverse hessian (I to start)
        for(int k = 0; k < 2*ni; k++)
        {
            H(k,k) = 1.0;
        }

        Array<OneD,NekDouble> xi(ni*2);
        for(int k = np - ni; k < np; k++)
        {
            xi[(k-np+ni)*2+0] = uvi[k][0];
            xi[(k-np+ni)*2+1] = uvi[k][1];
        }

        DNekMat J = opti->dF(xi);

        bool repeat = true;
        int itct = 0;
        while(repeat)
        {
            NekDouble Norm = 0;
            for(int k = 0; k < nq - 2; k++)
            {
                Norm += J(k,0)*J(k,0);
            }
            Norm = sqrt(Norm);

            if(Norm < 1E-8)
            {
                repeat = false;
                break;
            }
            if(itct > 100)
            {
                cout << "failed to optimise on face " << s->GetId() << endl;
                exit(-1);
                break;
            }
            itct++;
            cout << "Norm " << Norm << endl;

            BGFSUpdate(opti, J, B, H); //all will be updated
        }

        uvi = opti->GetSolution();

        vector<NodeSharedPtr> honodes;
        for(int j = np - ni; j < np; j++)
        {
            Array<OneD, NekDouble> loc;
            loc = s->P(uvi[j]);
            NodeSharedPtr nn = boost::shared_ptr<Node>(new
        Node(0,loc[0],loc[1],loc[2]));
            nn->SetCADSurf(surf, s, uvi[j]);
            honodes.push_back(nn);
        }

        f->m_faceNodes = honodes;
        f->m_curveType = LibUtilities::eNodalTriFekete;*/
    }

    if (m_mesh->m_verbose)
        cout << endl;
}

void Nektar::NekMeshUtils::SurfaceMesh::Mesh

(

)

Run all linear meshing routines.

Definition at line 49 of file SurfaceMesh.cpp.

References Nektar::StdRegions::find(), Nektar::iterator, and Nektar::LibUtilities::PrintProgressbar().

{
    if (m_mesh->m_verbose)
        cout << endl << "Surface meshing" << endl;

    if (m_mesh->m_verbose)
        cout << endl << "\tCurve meshing:" << endl << endl;

    m_mesh->m_numNodes = m_cad->GetNumVerts();

    // linear mesh all curves
    for (int i = 1; i <= m_cad->GetNumCurve(); i++)
    {
        if (m_mesh->m_verbose)
        {
            LibUtilities::PrintProgressbar(
                i, m_cad->GetNumCurve(), "Curve progress");
        }

        m_curvemeshes[i] = MemoryManager<CurveMesh>::AllocateSharedPtr(
            i, m_mesh, m_cad->GetCurve(i), m_octree);

        m_curvemeshes[i]->Mesh();
    }

    if (m_mesh->m_verbose)
        cout << endl << "\tFace meshing:" << endl << endl;

    // linear mesh all surfaces
    for (int i = 1; i <= m_cad->GetNumSurf(); i++)
    {
        if (m_mesh->m_verbose)
        {
            LibUtilities::PrintProgressbar(
                i, m_cad->GetNumSurf(), "Face progress");
        }

        vector<unsigned int>::iterator it;

        it = find(m_symsurfs.begin(), m_symsurfs.end(), i);
        if (it == m_symsurfs.end())
        {
            m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
                i, m_mesh, m_cad->GetSurf(i), m_octree, m_curvemeshes);
        }
        else
        {
            m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
                i, m_mesh, m_cad->GetSurf(i), m_octree, m_curvemeshes, m_bl);
        }

        m_facemeshes[i]->Mesh();
    }
}

void Nektar::NekMeshUtils::SurfaceMesh::Report

(

)

Print brief information to screen.

Definition at line 126 of file SurfaceMesh.cpp.

{
    if (m_mesh->m_verbose)
    {
        int ns = m_mesh->m_vertexSet.size();
        int es = m_mesh->m_edgeSet.size();
        int ts = m_mesh->m_element[2].size();
        int ep = ns - es + ts;
        cout << endl << "\tSurface mesh statistics" << endl;
        cout << "\t\tNodes: " << ns << endl;
        cout << "\t\tEdges: " << es << endl;
        cout << "\t\tTriangles " << ts << endl;
        cout << "\t\tEuler-Poincaré characteristic: " << ep << endl;
    }
}

void Nektar::NekMeshUtils::SurfaceMesh::Validate

(

)

Validate the linear surface mesh.

Definition at line 105 of file SurfaceMesh.cpp.

References ASSERTL0, and Nektar::iterator.

{
    if (m_mesh->m_verbose)
        cout << endl << "\tVerifying surface mesh" << endl;

    EdgeSet::iterator it;

    for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end(); it++)
    {
        if ((*it)->m_elLink.size() != 2)
        {
            if (m_mesh->m_verbose)
                cout << "\t\tFailed" << endl;
            ASSERTL0(false, "edge not listed twice");
        }
    }

    if (m_mesh->m_verbose)
        cout << "\t\tPassed" << endl;
}

Friends And Related Function Documentation

friend class MemoryManager< SurfaceMesh >

friend

Definition at line 56 of file SurfaceMesh.h.

Member Data Documentation

NekDouble Nektar::NekMeshUtils::SurfaceMesh::m_bl

private

thickness of the boundary layer if needed

Definition at line 103 of file SurfaceMesh.h.

CADSystemSharedPtr Nektar::NekMeshUtils::SurfaceMesh::m_cad

private

CAD object.

Definition at line 93 of file SurfaceMesh.h.

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

private

map of individual curve meshes of the curves in the domain

Definition at line 99 of file SurfaceMesh.h.

std::map<int, FaceMeshSharedPtr> Nektar::NekMeshUtils::SurfaceMesh::m_facemeshes

private

map of individual surface meshes from parametric surfaces

Definition at line 97 of file SurfaceMesh.h.

MeshSharedPtr Nektar::NekMeshUtils::SurfaceMesh::m_mesh

private

mesh object

Definition at line 91 of file SurfaceMesh.h.

OctreeSharedPtr Nektar::NekMeshUtils::SurfaceMesh::m_octree

private

Octree object.

Definition at line 95 of file SurfaceMesh.h.

std::vector<unsigned int> Nektar::NekMeshUtils::SurfaceMesh::m_symsurfs

private

list of sysmetry plane surfaces to build quads onto

Definition at line 101 of file SurfaceMesh.h.