Nektar::Utilities::ProcessVarOpti Class Reference

#include <ProcessVarOpti.h>

Inheritance diagram for Nektar::Utilities::ProcessVarOpti:

Collaboration diagram for Nektar::Utilities::ProcessVarOpti:

Public Member Functions
	ProcessVarOpti (MeshSharedPtr m)
virtual	~ProcessVarOpti ()
virtual void	Process ()
Public Member Functions inherited from Nektar::NekMeshUtils::ProcessModule
NEKMESHUTILS_EXPORT	ProcessModule (MeshSharedPtr p_m)
Public Member Functions inherited from Nektar::NekMeshUtils::Module
NEKMESHUTILS_EXPORT	Module (MeshSharedPtr p_m)
NEKMESHUTILS_EXPORT void	RegisterConfig (std::string key, std::string value)
	Register a configuration option with a module. More...
NEKMESHUTILS_EXPORT void	PrintConfig ()
	Print out all configuration options for a module. More...
NEKMESHUTILS_EXPORT void	SetDefaults ()
	Sets default configuration options for those which have not been set. More...
NEKMESHUTILS_EXPORT MeshSharedPtr	GetMesh ()
virtual NEKMESHUTILS_EXPORT void	ProcessVertices ()
	Extract element vertices. More...
virtual NEKMESHUTILS_EXPORT void	ProcessEdges (bool ReprocessEdges=true)
	Extract element edges. More...
virtual NEKMESHUTILS_EXPORT void	ProcessFaces (bool ReprocessFaces=true)
	Extract element faces. More...
virtual NEKMESHUTILS_EXPORT void	ProcessElements ()
	Generate element IDs. More...
virtual NEKMESHUTILS_EXPORT void	ProcessComposites ()
	Generate composites. More...
virtual NEKMESHUTILS_EXPORT void	ClearElementLinks ()

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

Static Public Attributes
static ModuleKey	className

Private Types
typedef std::map< int, std::vector< ElUtilSharedPtr > >	NodeElMap

Private Member Functions
void	Analytics ()
std::map < LibUtilities::ShapeType, DerivUtilSharedPtr >	BuildDerivUtil (int o)
void	GetElementMap (int o, std::map< LibUtilities::ShapeType, DerivUtilSharedPtr > derMap)
std::vector< ElUtilSharedPtr >	GetLockedElements (NekDouble thres)
std::vector< std::vector < NodeSharedPtr > >	CreateColoursets (std::vector< NodeSharedPtr > remain)
std::vector< std::vector < NodeSharedPtr > >	GetColouredNodes (std::vector< ElUtilSharedPtr > elLock)
void	RemoveLinearCurvature ()

Private Attributes
NodeElMap	m_nodeElMap
std::vector< ElUtilSharedPtr >	m_dataSet
ResidualSharedPtr	m_res
optiType	m_opti

Additional Inherited Members
Protected Member Functions inherited from Nektar::NekMeshUtils::Module
NEKMESHUTILS_EXPORT void	ReorderPrisms (PerMap &perFaces)
	Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...
NEKMESHUTILS_EXPORT void	PrismLines (int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
Protected Attributes inherited from Nektar::NekMeshUtils::Module
MeshSharedPtr	m_mesh
	Mesh object. More...
std::map< std::string, ConfigOption >	m_config
	List of configuration values. More...

Detailed Description

Definition at line 83 of file ProcessVarOpti.h.

Member Typedef Documentation

typedef std::map<int, std::vector<ElUtilSharedPtr> > Nektar::Utilities::ProcessVarOpti::NodeElMap

private

Definition at line 101 of file ProcessVarOpti.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessVarOpti::ProcessVarOpti

(

MeshSharedPtr

m

)

Definition at line 67 of file ProcessVarOpti.cpp.

References Nektar::NekMeshUtils::Module::m_config.

                                              : ProcessModule(m)
{
    // clang-format off
    m_config["linearelastic"] =
        ConfigOption(true, "", "Optimise for linear elasticity");
    m_config["winslow"] =
        ConfigOption(true, "", "Optimise for winslow");
    m_config["roca"] =
        ConfigOption(true, "", "Optimise for roca method");
    m_config["hyperelastic"] =
        ConfigOption(true, "", "Optimise for hyper elasticity");
    m_config["numthreads"] =
        ConfigOption(false, "1", "Number of threads");
    m_config["restol"] =
        ConfigOption(false, "1e-6", "Tolerance criterion");
    m_config["maxiter"] =
        ConfigOption(false, "500", "Maximum number of iterations");
    m_config["nq"] =
        ConfigOption(false, "-1", "Order of mesh");
    m_config["region"] =
        ConfigOption(false, "0.0", "create regions based on target");
    m_config["resfile"] =
        ConfigOption(false, "", "writes residual values to file");
    m_config["histfile"] =
        ConfigOption(false, "", "histogram of scaled jac");
    m_config["overint"] =
        ConfigOption(false, "6", "over integration order");
    m_config["analytics"] =
        ConfigOption(false, "", "basic analytics module");
    // clang-format on
}

Nektar::Utilities::ProcessVarOpti::~ProcessVarOpti ( )

virtual

Definition at line 99 of file ProcessVarOpti.cpp.

{
}

Member Function Documentation

void Nektar::Utilities::ProcessVarOpti::Analytics ( )

private

Definition at line 447 of file ProcessVarOpti.cpp.

References BuildDerivUtil(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), GetElementMap(), Nektar::Utilities::GetNodeOptiFactory(), m_dataSet, Nektar::NekMeshUtils::Module::m_mesh, m_nodeElMap, m_opti, and m_res.

Referenced by Process().

{
    // Grab the first element from the list
    ElementSharedPtr elmt = m_mesh->m_element[m_mesh->m_expDim][0];

    // Get curved nodes
    vector<NodeSharedPtr> nodes;
    elmt->GetCurvedNodes(nodes);

    // We're going to investigate only the first node (corner node)
    NodeSharedPtr node = nodes[4];

    // Loop over overintegration orders
    const int nPoints       = 200;
    const int overInt       = 40;
    const NekDouble originX = -1.0;
    const NekDouble originY = -1.0;
    const NekDouble length  = 2.0;
    const NekDouble dx      = length / (nPoints - 1);

    cout << "# overint = " << overInt << endl;
    cout << "# Columns: x, y, over-integration orders (0 -> " << overInt - 1
         << "), "
         << " min(scaledJac)" << endl;

    // Loop over square defined by (originX, originY), length
    for (int k = 0; k < nPoints; ++k)
    {
        node->m_y = originY + k * dx;
        for (int j = 0; j < nPoints; ++j)
        {
            node->m_x = originX + j * dx;
            cout << node->m_x << " " << node->m_y << " ";

            NekDouble minJacNew;

            for (int i = 0; i < overInt; ++i)
            {
                // Clear any existing node to element mapping.
                m_dataSet.clear();
                m_nodeElMap.clear();

                // Build deriv utils and element map.
                map<LibUtilities::ShapeType, DerivUtilSharedPtr> derivUtils =
                    BuildDerivUtil(i);

                // Reconstruct element map
                GetElementMap(i, derivUtils);

                for (int j = 0; j < m_dataSet.size(); j++)
                {
                    m_dataSet[j]->Evaluate();
                    m_dataSet[j]->InitialMinJac();
                }

                // Create NodeOpti object.
                NodeOptiSharedPtr nodeOpti =
                    GetNodeOptiFactory().CreateInstance(
                        m_mesh->m_spaceDim * 11, node,
                        m_nodeElMap.find(node->m_id)->second, m_res, derivUtils,
                        m_opti);

                minJacNew = 0.0;

                // Evaluate functional.
                nodeOpti->CalcMinJac();
                cout << nodeOpti->GetFunctional<2>(minJacNew) << " ";
                // NekDouble eigen;
                // nodeOpti->GetFunctional<2>(minJacNew);
                // nodeOpti->MinEigen<2>(eigen);
                // cout << eigen << " ";
            }

            cout << minJacNew << endl;
        }
    }
}

map< LibUtilities::ShapeType, DerivUtilSharedPtr > Nektar::Utilities::ProcessVarOpti::BuildDerivUtil ( int  o )

private

Definition at line 49 of file PreProcessing.cpp.

References ASSERTL0, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::eNodalPrismElec, Nektar::LibUtilities::eNodalPrismSPI, Nektar::LibUtilities::eNodalQuadElec, Nektar::LibUtilities::eNodalTetElec, Nektar::LibUtilities::eNodalTetSPI, Nektar::LibUtilities::eNodalTriElec, Nektar::LibUtilities::eNodalTriSPI, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTetrahedron, Nektar::LibUtilities::eTriangle, Nektar::LibUtilities::NodalUtil::GetInterpolationMatrix(), Nektar::LibUtilities::NodalUtil::GetVandermonde(), Nektar::LibUtilities::NodalUtil::GetVandermondeForDeriv(), Nektar::iterator, and Nektar::LibUtilities::PointsManager().

Referenced by Analytics(), and Process().

{
    // build Vandermonde information
    map<LibUtilities::ShapeType, DerivUtilSharedPtr> ret;

    // Typedef for points types used in the variational optimser. First entry is
    // the evaluation points; second entry is the distributions used for the
    // full element (includes element boundary)
    typedef std::pair<LibUtilities::PointsType, LibUtilities::PointsType>
        PTypes;

    map<LibUtilities::ShapeType, PTypes> typeMap;
    map<LibUtilities::ShapeType, PTypes>::iterator it;

    if (m_mesh->m_nummode + o <= 11)
    {
        typeMap[LibUtilities::eTriangle] =
            PTypes(LibUtilities::eNodalTriSPI, LibUtilities::eNodalTriElec);
        typeMap[LibUtilities::eTetrahedron] =
            PTypes(LibUtilities::eNodalTetSPI, LibUtilities::eNodalTetElec);
        typeMap[LibUtilities::ePrism] =
            PTypes(LibUtilities::eNodalPrismSPI, LibUtilities::eNodalPrismElec);
    }

    typeMap[LibUtilities::eQuadrilateral] =
        PTypes(LibUtilities::eNodalQuadElec, LibUtilities::eNodalQuadElec);
    // typeMap[LibUtilities::eHexahedron] =
    //    PTypes(LibUtilities::eNodalHexElec, LibUtilities::eNodalHexElec);

    for (it = typeMap.begin(); it != typeMap.end(); it++)
    {
        PTypes pType           = it->second;
        DerivUtilSharedPtr der = boost::shared_ptr<DerivUtil>(new DerivUtil());

        LibUtilities::PointsKey pkey1(m_mesh->m_nummode, pType.second);
        LibUtilities::PointsKey pkey2(m_mesh->m_nummode + o, pType.first);

        const int pDim  = pkey1.GetPointsDim();
        const int order = m_mesh->m_nummode - 1;

        Array<OneD, Array<OneD, NekDouble> > u1(pDim), u2(pDim);

        switch (pDim)
        {
            case 2:
            {
                LibUtilities::PointsManager()[pkey1]->GetPoints(u1[0], u1[1]);
                LibUtilities::PointsManager()[pkey2]->GetPoints(u2[0], u2[1]);
                break;
            }
            case 3:
            {
                LibUtilities::PointsManager()[pkey1]->GetPoints(u1[0], u1[1],
                                                                u1[2]);
                LibUtilities::PointsManager()[pkey2]->GetPoints(u2[0], u2[1],
                                                                u2[2]);
                break;
            }
        }

        der->ptsStd = u1[0].num_elements();
        der->pts    = u2[0].num_elements();

        LibUtilities::NodalUtil *nodalUtil = NULL;

        if (it->first == LibUtilities::eTriangle)
        {
            nodalUtil =
                new LibUtilities::NodalUtilTriangle(order, u1[0], u1[1]);
        }
        else if (it->first == LibUtilities::eQuadrilateral)
        {
            nodalUtil = new LibUtilities::NodalUtilQuad(order, u1[0], u1[1]);
        }
        else if (it->first == LibUtilities::eTetrahedron)
        {
            nodalUtil = new LibUtilities::NodalUtilTetrahedron(order, u1[0],
                                                               u1[1], u1[2]);
        }
        else if (it->first == LibUtilities::ePrism)
        {
            nodalUtil =
                new LibUtilities::NodalUtilPrism(order, u1[0], u1[1], u1[2]);
        }
        else if (it->first == LibUtilities::eHexahedron)
        {
            nodalUtil =
                new LibUtilities::NodalUtilHex(order, u1[0], u1[1], u1[2]);
        }
        else
        {
            ASSERTL0(false,
                     "Unknown element type for derivative utility setup");
        }

        NekMatrix<NekDouble> interp = *nodalUtil->GetInterpolationMatrix(u2);
        NekMatrix<NekDouble> Vandermonde  = *nodalUtil->GetVandermonde();
        NekMatrix<NekDouble> VandermondeI = Vandermonde;
        VandermondeI.Invert();

        for (int i = 0; i < pDim; ++i)
        {
            der->VdmDStd[i] =
                *nodalUtil->GetVandermondeForDeriv(i) * VandermondeI;
            der->VdmD[i] = interp * der->VdmDStd[i];
        }

        Array<OneD, NekDouble> qds =
            LibUtilities::PointsManager()[pkey2]->GetW();
        NekVector<NekDouble> quadWi(qds);
        der->quadW = quadWi;

        ret[it->first] = der;
        //delete nodalUtil;
    }

    return ret;
}

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

inlinestatic

Creates an instance of this class.

Definition at line 87 of file ProcessVarOpti.h.

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

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

vector< vector< NodeSharedPtr > > Nektar::Utilities::ProcessVarOpti::CreateColoursets ( std::vector< NodeSharedPtr >  remain )

private

Definition at line 455 of file PreProcessing.cpp.

References ASSERTL0, Nektar::iterator, and Nektar::LibUtilities::PrintProgressbar().

{
    vector<vector<NodeSharedPtr> > retPart;

    // loop until all free nodes have been sorted
    while (remain.size() > 0)
    {
        vector<NodeSharedPtr> layer;  // one colourset
        set<int> locked;
        set<int> completed;
        for (int i = 0; i < remain.size(); i++)
        {
            NodeElMap::iterator it = m_nodeElMap.find(remain[i]->m_id); //try to find node within all elements
            ASSERTL0(it != m_nodeElMap.end(), "could not find node");

            vector<ElUtilSharedPtr> &elUtils = it->second; // identify the vector of all associated elements of the node

            bool islocked = false; // suppose node is not locked
            for (int j = 0; j < elUtils.size(); j++) // loop over all associated elements of the node
            {
                set<int>::iterator sit = locked.find(elUtils[j]->GetId()); // check all nodes of the element
                if (sit != locked.end())                          //if node is within the set of locked nodes
                {
                    islocked = true; // if yes, flag node as locked
                    break;              // and go to next node
                }
            }
            if (!islocked) // if the node is not locked
            {
                layer.push_back(remain[i]);   // insert node into colourset
                completed.insert(remain[i]->m_id);    // insert sorted node into "completed" list
                for (int j = 0; j < elUtils.size(); j++)   // loop over all other nodes of the same element
                {
                    locked.insert(elUtils[j]->GetId());   // and flag these nodes as locked
                }
            }
        }

        // identify nodes which are not sorted, yet and create new "remain" vector
        vector<NodeSharedPtr> tmp = remain;
        remain.clear();
        for (int i = 0; i < tmp.size(); i++)
        {
            set<int>::iterator sit = completed.find(tmp[i]->m_id);
            if (sit == completed.end())
            {
                remain.push_back(tmp[i]);
            }
        }

        // include layer or colourset into vector of coloursets
        retPart.push_back(layer);

        // print out progress
        if(m_mesh->m_verbose)
        {
            LibUtilities::PrintProgressbar(m_res->n - remain.size(), m_res->n, "Node Coloring");
        }

    }
    return retPart;
}

vector< vector< NodeSharedPtr > > Nektar::Utilities::ProcessVarOpti::GetColouredNodes ( std::vector< ElUtilSharedPtr >  elLock )

private

Definition at line 186 of file PreProcessing.cpp.

References ASSERTL0, and Nektar::iterator.

Referenced by Process().

{

    // create set of nodes to be ignored and hence not included in the coloursets
    NodeSet ignoredNodes;
    for (int i = 0; i < elLock.size(); i++)
    {
        vector<NodeSharedPtr> nodes;
        elLock[i]->GetEl()->GetCurvedNodes(nodes);
        for (int j = 0; j < nodes.size(); j++)
        {
            ignoredNodes.insert(nodes[j]);
        }
    }

    // create set of nodes which are at the boundary and hence not included in the colourset
    NodeSet boundaryNodes;

    switch (m_mesh->m_spaceDim)
    {
        case 2:
        {
            EdgeSet::iterator it;
            for(it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end(); it++)
            {
                if((*it)->m_elLink.size() == 2)
                {
                    continue;
                }

                boundaryNodes.insert((*it)->m_n1);
                boundaryNodes.insert((*it)->m_n2);
                for(int i = 0; i < (*it)->m_edgeNodes.size(); i++)
                {
                    boundaryNodes.insert((*it)->m_edgeNodes[i]);
                }
            }
            break;
        }
        case 3:
        {
            if(!m_mesh->m_cad)
            {
                FaceSet::iterator it;
                for (it = m_mesh->m_faceSet.begin();
                     it != m_mesh->m_faceSet.end(); it++)
                {
                    if ((*it)->m_elLink.size() == 2)
                    {
                        continue;
                    }

                    vector<NodeSharedPtr> vs = (*it)->m_vertexList;
                    for (int j = 0; j < vs.size(); j++)
                    {
                        boundaryNodes.insert(vs[j]);
                    }

                    vector<EdgeSharedPtr> es = (*it)->m_edgeList;
                    for (int j = 0; j < es.size(); j++)
                    {
                        for (int k = 0; k < es[j]->m_edgeNodes.size(); k++)
                        {
                            boundaryNodes.insert(es[j]->m_edgeNodes[k]);
                        }
                    }

                    for (int i = 0; i < (*it)->m_faceNodes.size(); i++)
                    {
                        boundaryNodes.insert((*it)->m_faceNodes[i]);
                    }
                }
            }
            else
            {
                //if we have CAD therefore the only fixed nodes exist on vertices only
                NodeSet::iterator nit;
                for (nit = m_mesh->m_vertexSet.begin(); nit != m_mesh->m_vertexSet.end(); ++nit)
                {
                    if((*nit)->GetNumCadCurve() > 1)
                    {
                        boundaryNodes.insert((*nit));
                    }
                }
            }
            break;
        }
        default:
            ASSERTL0(false,"space dim issue");
    }


    //create vector of free nodes which "remain", hence will be included in the coloursets
    vector<NodeSharedPtr> remainEdgeVertex;
    vector<NodeSharedPtr> remainFace;
    vector<NodeSharedPtr> remainVolume;
    m_res->nDoF = 0;

    // check if vertex nodes are in boundary or ignored nodes, otherwise add to EDGE-VERTEX remain nodes
    NodeSet::iterator nit;
    for (nit = m_mesh->m_vertexSet.begin(); nit != m_mesh->m_vertexSet.end();
         ++nit)
    {
        NodeSet::iterator nit2 = boundaryNodes.find(*nit);
        NodeSet::iterator nit3 = ignoredNodes.find(*nit);
        if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
        {
            remainEdgeVertex.push_back(*nit);
            if ((*nit)->GetNumCadCurve() == 1)
            {
                m_res->nDoF++;
            }
            else if ((*nit)->GetNumCADSurf() == 1)
            {
                m_res->nDoF += 2;
            }
            else
            {
                m_res->nDoF += m_mesh->m_spaceDim;
            }
        }
    }

    // check if edge nodes are in boundary or ignored nodes, otherwise add to EDGE-VERTEX remain nodes
    EdgeSet::iterator eit;
    for (eit = m_mesh->m_edgeSet.begin(); eit != m_mesh->m_edgeSet.end(); eit++)
    {
        vector<NodeSharedPtr> n = (*eit)->m_edgeNodes;
        for (int j = 0; j < n.size(); j++)
        {
            NodeSet::iterator nit2 = boundaryNodes.find(n[j]);
            NodeSet::iterator nit3 = ignoredNodes.find(n[j]);
            if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
            {
                remainEdgeVertex.push_back(n[j]);
                if (n[j]->GetNumCadCurve() == 1)
                {
                    m_res->nDoF++;
                }
                else if (n[j]->GetNumCADSurf() == 1)
                {
                    m_res->nDoF += 2;
                }
                else
                {
                    m_res->nDoF += m_mesh->m_spaceDim;
                }
            }
        }
    }

    // check if face nodes are in boundary or ignored nodes, otherwise add to FACE remain nodes
    FaceSet::iterator fit;
    for (fit = m_mesh->m_faceSet.begin(); fit != m_mesh->m_faceSet.end(); fit++)
    {
        for (int j = 0; j < (*fit)->m_faceNodes.size(); j++)
        {
            NodeSet::iterator nit2 = boundaryNodes.find((*fit)->m_faceNodes[j]);
            NodeSet::iterator nit3 = ignoredNodes.find((*fit)->m_faceNodes[j]);
            if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
            {
                remainFace.push_back((*fit)->m_faceNodes[j]);
                if ((*fit)->m_faceNodes[j]->GetNumCADSurf() == 1)
                {
                    m_res->nDoF += 2;
                }
                else
                {
                    m_res->nDoF += m_mesh->m_spaceDim;
                }
            }
        }
    }

    // check if volume nodes are in boundary or ignored nodes, otherwise add to VOLUME remain nodes
    for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); i++)
    {
        vector<NodeSharedPtr> ns =
            m_mesh->m_element[m_mesh->m_expDim][i]->GetVolumeNodes();
        for (int j = 0; j < ns.size(); j++)
        {
            NodeSet::iterator nit2 = boundaryNodes.find(ns[j]);
            NodeSet::iterator nit3 = ignoredNodes.find(ns[j]);
            if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
            {
                remainVolume.push_back(ns[j]);
                m_res->nDoF += m_mesh->m_spaceDim;
            }
        }
    }

    // size of all free nodes to be included in the coloursets
    m_res->n = remainEdgeVertex.size()
                + remainFace.size() + remainVolume.size();

    // data structure for coloursets, that will ultimately contain all free nodes
    vector<vector<NodeSharedPtr> > ret;
    vector<vector<NodeSharedPtr> > retPart;


    // edge and vertex nodes
    // create vector num_el of number of associated elements of each node
    vector<int> num_el(remainEdgeVertex.size());
    for (int i = 0; i < remainEdgeVertex.size(); i++)
    {
        //try to find node within all elements
        NodeElMap::iterator it = m_nodeElMap.find(remainEdgeVertex[i]->m_id);
        vector<ElUtilSharedPtr> &elUtils = it->second;
        num_el[i] = elUtils.size();
    }
    // finding the permutation according to num_el
    vector<int> permNode(remainEdgeVertex.size());
    for (int i = 0; i < remainEdgeVertex.size(); ++i)
    {
        permNode[i] = i;
    }
    std::sort(permNode.begin(), permNode.end(), NodeComparator(num_el));
    // applying the permutation to remainEdgeVertex
    vector<NodeSharedPtr> remainEdgeVertexSort(remainEdgeVertex.size());
    for (int i = 0; i < remainEdgeVertex.size(); ++i)
    {
        int j = permNode[i];
        remainEdgeVertexSort[i] = remainEdgeVertex[j];
    }

    retPart = CreateColoursets(remainEdgeVertexSort);
    if(m_mesh->m_verbose)
    {
        cout << "Number of Edge/Vertex Coloursets: " << retPart.size() << endl;
    }
    for (int i = 0; i < retPart.size(); i++)
    {
        ret.push_back(retPart[i]);
    }


    // face nodes
    retPart = CreateColoursets(remainFace);
    if(m_mesh->m_verbose)
    {
        cout << "Number of Face Coloursets: " << retPart.size() << endl;
    }
    for (int i = 0; i < retPart.size(); i++)
    {
        ret.push_back(retPart[i]);
    }


    // volume nodes
    retPart = CreateColoursets(remainVolume);
    if(m_mesh->m_verbose)
    {
        cout << "Number of Volume Coloursets: " << retPart.size() << endl;
    }
    for (int i = 0; i < retPart.size(); i++)
    {
        ret.push_back(retPart[i]);
    }


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

    return ret;
}

void Nektar::Utilities::ProcessVarOpti::GetElementMap ( int  o, std::map< LibUtilities::ShapeType, DerivUtilSharedPtr >  derMap  )

private

Definition at line 520 of file PreProcessing.cpp.

References ASSERTL0.

Referenced by Analytics(), and Process().

{
    for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); i++)
    {
        ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
        vector<NodeSharedPtr> ns;
        el->GetCurvedNodes(ns);
        ElUtilSharedPtr d = boost::shared_ptr<ElUtil>(new ElUtil(
            el, derMap[el->GetShapeType()], m_res, m_mesh->m_nummode, o));
        m_dataSet.push_back(d);
    }

    for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); i++)
    {
        ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
        vector<NodeSharedPtr> ns;
        el->GetCurvedNodes(ns);

        for (int j = 0; j < ns.size(); j++)
        {
            m_nodeElMap[ns[j]->m_id].push_back(m_dataSet[i]);
        }

        ASSERTL0(derMap[el->GetShapeType()]->ptsStd == ns.size(),
                 "mismatch node count");
    }
}

vector< ElUtilSharedPtr > Nektar::Utilities::ProcessVarOpti::GetLockedElements ( NekDouble  thres )

private

Definition at line 549 of file PreProcessing.cpp.

References Nektar::iterator.

Referenced by Process().

{
    vector<ElUtilSharedPtr> elBelowThres;
    for (int i = 0; i < m_dataSet.size(); ++i)
    {
        if (m_dataSet[i]->GetScaledJac() < thres)
        {
            elBelowThres.push_back(m_dataSet[i]);
        }
    }

    boost::unordered_set<int> inmesh;
    pair<boost::unordered_set<int>::iterator, bool> t;
    vector<ElUtilSharedPtr> totest;

    for (int i = 0; i < elBelowThres.size(); i++)
    {
        t = inmesh.insert(elBelowThres[i]->GetId());

        vector<FaceSharedPtr> f = elBelowThres[i]->GetEl()->GetFaceList();
        for (int j = 0; j < f.size(); j++)
        {
            for (int k = 0; k < f[j]->m_elLink.size(); k++)
            {
                if (f[j]->m_elLink[k].first->GetId() ==
                    elBelowThres[i]->GetId())
                    continue;

                t = inmesh.insert(f[j]->m_elLink[k].first->GetId());
                if (t.second)
                {
                    totest.push_back(
                        m_dataSet[f[j]->m_elLink[k].first->GetId()]);
                }
            }
        }
    }

    for (int i = 0; i < 6; i++)
    {
        vector<ElUtilSharedPtr> tmp = totest;
        totest.clear();
        for (int j = 0; j < tmp.size(); j++)
        {
            vector<FaceSharedPtr> f = tmp[j]->GetEl()->GetFaceList();
            for (int k = 0; k < f.size(); k++)
            {
                for (int l = 0; l < f[k]->m_elLink.size(); l++)
                {
                    if (f[k]->m_elLink[l].first->GetId() == tmp[j]->GetId())
                        continue;

                    t = inmesh.insert(f[k]->m_elLink[l].first->GetId());
                    if (t.second)
                    {
                        totest.push_back(
                            m_dataSet[f[k]->m_elLink[l].first->GetId()]);
                    }
                }
            }
        }
    }

    // now need to invert the list
    vector<ElUtilSharedPtr> ret;
    for (int i = 0; i < m_dataSet.size(); ++i)
    {
        boost::unordered_set<int>::iterator s =
            inmesh.find(m_dataSet[i]->GetId());
        if (s == inmesh.end())
        {
            ret.push_back(m_dataSet[i]);
        }
    }

    return ret;
}

void Nektar::Utilities::ProcessVarOpti::Process ( )

virtual

Implements Nektar::NekMeshUtils::Module.

Definition at line 103 of file ProcessVarOpti.cpp.

References Analytics(), ASSERTL0, BuildDerivUtil(), Nektar::Thread::ThreadMaster::CreateInstance(), Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::Utilities::eHypEl, Nektar::Utilities::eLinEl, Nektar::Utilities::eRoca, Nektar::Utilities::eWins, GetColouredNodes(), GetElementMap(), GetLockedElements(), Nektar::Utilities::GetNodeOptiFactory(), Nektar::iterator, Nektar::NekMeshUtils::Module::m_config, m_dataSet, Nektar::NekMeshUtils::Module::m_mesh, m_nodeElMap, m_opti, m_res, CellMLToNektar.cellml_metadata::p, Nektar::Thread::ThreadMaster::SessionJob, Nektar::Thread::ThreadMaster::SetThreadingType(), Nektar::Timer::Start(), Nektar::Timer::Stop(), and Nektar::Timer::TimePerTest().

{
    if (m_mesh->m_verbose)
    {
        cout << "ProcessVarOpti: Optimising... " << endl;
    }

    if (m_config["linearelastic"].beenSet)
    {
        m_opti = eLinEl;
    }
    else if (m_config["winslow"].beenSet)
    {
        m_opti = eWins;
    }
    else if (m_config["roca"].beenSet)
    {
        m_opti = eRoca;
    }
    else if (m_config["hyperelastic"].beenSet)
    {
        m_opti = eHypEl;
    }
    else
    {
        ASSERTL0(false, "not opti type set");
    }

    const int maxIter      = m_config["maxiter"].as<int>();
    const NekDouble restol = m_config["restol"].as<NekDouble>();

    // m_mesh->m_nummode = m_config["nq"].as<int>();

    EdgeSet::iterator eit;
    bool fd = false;

    if (m_config["nq"].beenSet)
    {
        m_mesh->m_nummode = m_config["nq"].as<int>();
        fd                = true;
    }

    if (!fd)
    {
        for (eit = m_mesh->m_edgeSet.begin(); eit != m_mesh->m_edgeSet.end();
             eit++)
        {
            if ((*eit)->m_edgeNodes.size() > 0)
            {
                m_mesh->m_nummode = (*eit)->m_edgeNodes.size() + 2;
                fd                = true;
                break;
            }
        }
    }
    ASSERTL0(fd, "failed to find order of mesh");

    // Safety feature: limit over-integration order for high-order triangles
    // over order 5.
    int intOrder = m_config["overint"].as<int>();
    intOrder = m_mesh->m_nummode + intOrder <= 11 ?
        intOrder : 11 - m_mesh->m_nummode;

    if (m_mesh->m_verbose)
    {
        cout << "Identified mesh order as: " << m_mesh->m_nummode - 1 << endl;
    }

    if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3)
    {
        ASSERTL0(false, "cannot deal with manifolds");
    }

    m_res      = boost::shared_ptr<Residual>(new Residual);
    m_res->val = 1.0;

    
    m_mesh->MakeOrder(m_mesh->m_nummode - 1,
                      LibUtilities::eGaussLobattoLegendre);

    if (m_config["analytics"].beenSet)
    {
        Analytics();
        return;
    }

    map<LibUtilities::ShapeType, DerivUtilSharedPtr> derivUtils =
        BuildDerivUtil(intOrder);

    GetElementMap(intOrder, derivUtils);

    m_res->startInv = 0;
    m_res->worstJac = numeric_limits<double>::max();
    for (int i = 0; i < m_dataSet.size(); i++)
    {
        m_dataSet[i]->Evaluate();
        m_dataSet[i]->InitialMinJac();
    }

    vector<ElUtilSharedPtr> elLock;

    if (m_config["region"].beenSet)
    {
        elLock = GetLockedElements(m_config["region"].as<NekDouble>());
    }

    
    vector<vector<NodeSharedPtr> > freenodes = GetColouredNodes(elLock);
    vector<vector<NodeOptiSharedPtr> > optiNodes;
    
    // turn the free nodes into optimisable objects with all required data
    set<int> check;
    for (int i = 0; i < freenodes.size(); i++)
    {
        vector<NodeOptiSharedPtr> ns;
        for (int j = 0; j < freenodes[i].size(); j++)
        {
            NodeElMap::iterator it = m_nodeElMap.find(freenodes[i][j]->m_id);
            ASSERTL0(it != m_nodeElMap.end(), "could not find");

            int optiKind = m_mesh->m_spaceDim;

            if (freenodes[i][j]->GetNumCadCurve() == 1)
            {
                optiKind += 10;
            }
            else if (freenodes[i][j]->GetNumCADSurf() == 1)
            {
                optiKind += 20;
            }
            else
            {
                optiKind += 10 * m_mesh->m_expDim;
            }

            set<int>::iterator c = check.find(freenodes[i][j]->m_id);
            ASSERTL0(c == check.end(), "duplicate node");
            check.insert(freenodes[i][j]->m_id);

            ns.push_back(GetNodeOptiFactory().CreateInstance(
                optiKind, freenodes[i][j], it->second, m_res, derivUtils,
                m_opti));
        }
        optiNodes.push_back(ns);
    }

    int nset = optiNodes.size();
    int p    = 0;
    int mn   = numeric_limits<int>::max();
    int mx   = 0;
    for (int i = 0; i < nset; i++)
    {
        p += optiNodes[i].size();
        mn = min(mn, int(optiNodes[i].size()));
        mx = max(mx, int(optiNodes[i].size()));
    }

    if (m_config["histfile"].beenSet)
    {
        ofstream histFile;
        string name = m_config["histfile"].as<string>() + "_start.txt";
        histFile.open(name.c_str());

        for (int i = 0; i < m_dataSet.size(); i++)
        {
            histFile << m_dataSet[i]->GetScaledJac() << endl;
        }
        histFile.close();
    }

    if(m_mesh->m_verbose)
    {
        cout << scientific << endl;
        cout << "N elements:\t\t" << m_mesh->m_element[m_mesh->m_expDim].size() - elLock.size() << endl
             << "N elements invalid:\t" << m_res->startInv << endl
             << "Worst jacobian:\t\t" << m_res->worstJac << endl
             << "N free nodes:\t\t" << m_res->n << endl
             << "N Dof:\t\t\t" << m_res->nDoF << endl
             << "N color sets:\t\t" << nset << endl
             << "Avg set colors:\t\t" << p/nset << endl
             << "Min set:\t\t" << mn << endl
             << "Max set:\t\t" << mx << endl
             << "Residual tolerance:\t" << restol << endl;
    }

    int nThreads = m_config["numthreads"].as<int>();

    int ctr = 0;
    Thread::ThreadMaster tms;
    tms.SetThreadingType("ThreadManagerBoost");
    Thread::ThreadManagerSharedPtr tm =
        tms.CreateInstance(Thread::ThreadMaster::SessionJob, nThreads);

    Timer t;
    t.Start();

    ofstream resFile;
    if (m_config["resfile"].beenSet)
    {
        resFile.open(m_config["resfile"].as<string>().c_str());
    }

    for (int i = 0; i < optiNodes.size(); i++)
    {
        vector<Thread::ThreadJob *> jobs(optiNodes[i].size());
        for (int j = 0; j < optiNodes[i].size(); j++)
        {
            optiNodes[i][j]->CalcMinJac();
        }
    }

    while (m_res->val > restol)
    {
        ctr++;
        m_res->val       = 0.0;
        m_res->func      = 0.0;
        m_res->nReset[0] = 0;
        m_res->nReset[1] = 0;
        m_res->nReset[2] = 0;
        m_res->alphaI    = 0;
        for (int i = 0; i < optiNodes.size(); i++)
        {
            vector<Thread::ThreadJob *> jobs(optiNodes[i].size());
            for (int j = 0; j < optiNodes[i].size(); j++)
            {
                jobs[j] = optiNodes[i][j]->GetJob();
            }

            tm->SetNumWorkers(0);
            tm->QueueJobs(jobs);
            tm->SetNumWorkers(nThreads);
            tm->Wait();
        }

        m_res->startInv = 0;
        m_res->worstJac = numeric_limits<double>::max();

        vector<Thread::ThreadJob *> elJobs(m_dataSet.size());
        for (int i = 0; i < m_dataSet.size(); i++)
        {
            elJobs[i] = m_dataSet[i]->GetJob();
        }

        tm->SetNumWorkers(0);
        tm->QueueJobs(elJobs);
        tm->SetNumWorkers(nThreads);
        tm->Wait();

        if (m_config["resfile"].beenSet)
        {
            resFile << m_res->val << " " << m_res->worstJac << " "
                    << m_res->func << endl;
        }

        if(m_mesh->m_verbose)
        {
            cout << ctr
                 << "\tResidual: " << m_res->val
                 << "\tMin Jac: " << m_res->worstJac
                 << "\tInvalid: " << m_res->startInv
                 << "\tReset nodes: " << m_res->nReset[0] << "/" << m_res->nReset[1]
                 << "/" << m_res->nReset[2] << "\tFunctional: " << m_res->func
                 << endl;
        }

        if(ctr >= maxIter)
        {
            break;
        }
    }

    if (m_config["histfile"].beenSet)
    {
        ofstream histFile;
        string name = m_config["histfile"].as<string>() + "_end.txt";
        histFile.open(name.c_str());

        for (int i = 0; i < m_dataSet.size(); i++)
        {
            histFile << m_dataSet[i]->GetScaledJac() << endl;
        }
        histFile.close();
    }
    if (m_config["resfile"].beenSet)
    {
        resFile.close();
    }

    t.Stop();

    //RemoveLinearCurvature();

    if(m_mesh->m_verbose)
    {
        cout << "Time to compute: " << t.TimePerTest(1) << endl;
        cout << "Invalid at end: " << m_res->startInv << endl;
        cout << "Worst at end: " << m_res->worstJac << endl;
    }
}

void Nektar::Utilities::ProcessVarOpti::RemoveLinearCurvature ( )

private

Definition at line 627 of file PreProcessing.cpp.

References ASSERTL0, and Nektar::iterator.

{
    for(int i = 0; i < m_dataSet.size(); i++)
    {
        if(m_dataSet[i]->GetScaledJac() > 0.999)
        {
            ElementSharedPtr el = m_dataSet[i]->GetEl();
            vector<NodeSharedPtr> ns;
            el->SetVolumeNodes(ns);
        }
    }

    map<int, vector<FaceSharedPtr> > edgeToFace;

    FaceSet::iterator fit;
    for(fit = m_mesh->m_faceSet.begin(); fit != m_mesh->m_faceSet.end(); fit++)
    {
        bool rm = true;
        for(int i = 0; i < (*fit)->m_elLink.size(); i++)
        {
            int id = (*fit)->m_elLink[i].first->GetId();
            if(m_dataSet[id]->GetScaledJac() <= 0.999)
            {
                rm = false;
                break;
            }
        }
        if(rm)
        {
            (*fit)->m_faceNodes.clear();
        }

        vector<EdgeSharedPtr> es = (*fit)->m_edgeList;
        for(int i = 0; i < es.size(); i++)
        {
            edgeToFace[es[i]->m_id].push_back(*fit);
        }
    }

    EdgeSet::iterator eit;
    for(eit = m_mesh->m_edgeSet.begin(); eit != m_mesh->m_edgeSet.end(); eit++)
    {
        map<int, vector<FaceSharedPtr> >::iterator it;
        it = edgeToFace.find((*eit)->m_id);
        ASSERTL0(it != edgeToFace.end(),"not found");
        bool rm = true;
        for(int i = 0; i < it->second.size(); i++)
        {
            if(it->second[i]->m_faceNodes.size() > 0)
            {
                rm = false;
                break;
            }
        }
        if(rm)
        {
            (*eit)->m_edgeNodes.clear();
        }
    }
}

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessVarOpti::className

static

Initial value:

=
    GetModuleFactory().RegisterCreatorFunction(
        ModuleKey(eProcessModule, "varopti"), ProcessVarOpti::create,
        "Optimise mesh locations.")

Definition at line 91 of file ProcessVarOpti.h.

std::vector<ElUtilSharedPtr> Nektar::Utilities::ProcessVarOpti::m_dataSet

private

Definition at line 115 of file ProcessVarOpti.h.

Referenced by Analytics(), and Process().

NodeElMap Nektar::Utilities::ProcessVarOpti::m_nodeElMap

private

Definition at line 114 of file ProcessVarOpti.h.

Referenced by Analytics(), and Process().

optiType Nektar::Utilities::ProcessVarOpti::m_opti

private

Definition at line 118 of file ProcessVarOpti.h.

Referenced by Analytics(), and Process().

ResidualSharedPtr Nektar::Utilities::ProcessVarOpti::m_res

private

Definition at line 117 of file ProcessVarOpti.h.

Referenced by Analytics(), and Process().