Nektar::FieldUtils::Iso Class Reference

#include <ProcessIsoContour.h>

Public Member Functions
void	Condense (void)
void	GlobalCondense (std::vector< std::shared_ptr< Iso > > &iso, bool verbose)
void	SeparateRegions (std::vector< std::shared_ptr< Iso > > &iso, int minsize, bool verbose)
void	Smooth (int n_iter, NekDouble lambda, NekDouble mu)
int	GetNVert (void)
void	SetNVert (int n)
int	GetNTris (void)
void	SetNTris (int n)
void	SetFields (const int loc, const Array< OneD, Array< OneD, NekDouble > > &intfields, const int j)
NekDouble	GetFields (const int i, const int j)
void	SetX (int loc, NekDouble val)
void	SetY (int loc, NekDouble val)
void	SetZ (int loc, NekDouble val)
NekDouble	GetX (int loc)
NekDouble	GetY (int loc)
NekDouble	GetZ (int loc)
int	GetVId (int i)
void	ResizeVId (int nconn)
void	SetVId (int i, int j)
void	ResizeFields (int size)
	Iso (int nfields)
	~Iso (void)

Private Attributes
bool	m_condensed
int	m_nvert
int	m_ntris
std::vector< NekDouble >	m_x
std::vector< NekDouble >	m_y
std::vector< NekDouble >	m_z
std::vector< std::vector< NekDouble > >	m_fields
Array< OneD, int >	m_vid

Detailed Description

Definition at line 46 of file ProcessIsoContour.h.

Constructor & Destructor Documentation

Iso()

Nektar::FieldUtils::Iso::Iso ( int  nfields )

inline

Definition at line 155 of file ProcessIsoContour.h.

References m_condensed, m_fields, m_nvert, m_x, m_y, and m_z.

        {
            m_condensed = false;
            m_nvert     = 0;
            m_fields.resize(nfields);
            // set up initial vectors to be 1000 long
            m_x.resize(1000);
            m_y.resize(1000);
            m_z.resize(1000);
            for(int i = 0; i < m_fields.size(); ++i)
            {
                m_fields[i].resize(1000);
            }
        };

~Iso()

Nektar::FieldUtils::Iso::~Iso ( void  )

inline

Definition at line 170 of file ProcessIsoContour.h.

        {
        }

Member Function Documentation

Condense()

void Nektar::FieldUtils::Iso::Condense

(

void

)

Definition at line 672 of file ProcessIsoContour.cpp.

References Nektar::StdRegions::find(), Nektar::FieldUtils::IsoVertex::m_fields, Nektar::FieldUtils::IsoVertex::m_id, Nektar::FieldUtils::IsoVertex::m_x, Nektar::FieldUtils::IsoVertex::m_y, and Nektar::FieldUtils::IsoVertex::m_z.

{
    int i,j,cnt;
    IsoVertex v;
    vector<IsoVertex> vert;

    if(!m_ntris) return;

    if(m_condensed) return;
    m_condensed = true;

    vert.reserve(m_ntris/6);

    m_vid = Array<OneD, int>(3*m_ntris);

    // fill first 3 points and initialise fields
    v.m_fields.resize(m_fields.size());
    for(cnt =0, i = 0; i < 3; ++i)
    {
        v.m_x = m_x[i];
        v.m_y = m_y[i];
        v.m_z = m_z[i];
        for(int f = 0; f < m_fields.size(); ++f)
        {
            v.m_fields[f] = m_fields[f][i];
        }
        v.m_id = cnt;
        vert.push_back(v);
        m_vid[i] = v.m_id;
        ++cnt;
    }

    for(i = 1; i < m_ntris; ++i)
    {
        for(j = 0; j < 3; ++j)
        {
            v.m_x = m_x[3*i+j];
            v.m_y = m_y[3*i+j];
            v.m_z = m_z[3*i+j];

            auto pt = find(vert.begin(),vert.end(),v);
            if(pt != vert.end())
            {
                m_vid[3*i+j] = pt[0].m_id;
            }
            else
            {
                v.m_id = cnt;

                for(int f = 0; f < m_fields.size(); ++f)
                {
                    v.m_fields[f] = m_fields[f][3*i+j];
                }

                vert.push_back(v);

                m_vid[3*i+j] = v.m_id;
                ++cnt;
           }
        }
    }

    // remove elements with multiple vertices
    for(i = 0,cnt=0; i < m_ntris;)
    {
        if((m_vid[3*i]  ==m_vid[3*i+1])||
           (m_vid[3*i]  ==m_vid[3*i+2])||
           (m_vid[3*i+1]==m_vid[3*i+2]))
        {
            cnt++;
            for(j = 3*i; j < 3*(m_ntris-1); ++j)
            {
                m_vid[j] = m_vid[j+3];
            }
            m_ntris--;
        }
        else
        {
            ++i;
        }
    }

    m_nvert  = vert.size();

    m_x.resize(m_nvert);
    m_y.resize(m_nvert);
    m_z.resize(m_nvert);

    for(int f = 0; f < m_fields.size(); ++f)
    {
        m_fields[f].resize(m_nvert);
    }

    for(i = 0; i < m_nvert; ++i)
    {
        m_x[i] = vert[i].m_x;
        m_y[i] = vert[i].m_y;
        m_z[i] = vert[i].m_z;
        for(int f = 0; f < m_fields.size(); ++f)
        {
            m_fields[f][i] = vert[i].m_fields[f];
        }
    }
}

GetFields()

NekDouble Nektar::FieldUtils::Iso::GetFields ( const int  i, const int  j  )

inline

Definition at line 89 of file ProcessIsoContour.h.

References m_fields.

        {
            return m_fields[i][j];
        }

GetNTris()

int Nektar::FieldUtils::Iso::GetNTris ( void  )

inline

Definition at line 65 of file ProcessIsoContour.h.

References m_ntris.

        {
            return m_ntris;
        }

GetNVert()

int Nektar::FieldUtils::Iso::GetNVert ( void  )

inline

Definition at line 55 of file ProcessIsoContour.h.

References m_nvert.

        {
            return m_nvert;
        }

GetVId()

int Nektar::FieldUtils::Iso::GetVId ( int  i )

inline

Definition at line 124 of file ProcessIsoContour.h.

References m_vid.

        {
            return m_vid[i];
        }

GetX()

NekDouble Nektar::FieldUtils::Iso::GetX ( int  loc )

inline

Definition at line 109 of file ProcessIsoContour.h.

References CG_Iterations::loc, and m_x.

        {
            return m_x[loc];
        }

GetY()

NekDouble Nektar::FieldUtils::Iso::GetY ( int  loc )

inline

Definition at line 114 of file ProcessIsoContour.h.

References CG_Iterations::loc, and m_y.

        {
            return m_y[loc];
        }

GetZ()

NekDouble Nektar::FieldUtils::Iso::GetZ ( int  loc )

inline

Definition at line 119 of file ProcessIsoContour.h.

References CG_Iterations::loc, and m_z.

        {
            return m_z[loc];
        }

GlobalCondense()

void Nektar::FieldUtils::Iso::GlobalCondense	(	std::vector< std::shared_ptr< Iso > > &	iso,
		bool	verbose
	)

Definition at line 778 of file ProcessIsoContour.cpp.

References Nektar::NekConstants::kNekZeroTol, and Nektar::LibUtilities::PrintProgressbar().

{
    typedef bg::model::point<NekDouble, 3, bg::cs::cartesian> BPoint;
    typedef std::pair<BPoint, unsigned int> PointPair;

    int    i,j,n;
    int    nelmt;
    int    niso=iso.size();
    int    id1,id2;
    Array<OneD, Array<OneD, int> > vidmap;

    if(m_condensed) return;
    m_condensed = true;

    vidmap = Array<OneD, Array<OneD, int> > (niso);

    m_ntris = 0;
    for(i = 0; i < niso; ++i)
    {
        if(iso[i]->m_ntris)
        {
            m_ntris += iso[i]->m_ntris;
        }
    }

    m_vid = Array<OneD, int>(3*m_ntris);

    m_nvert = 0;
    for(i = 0; i < niso; ++i)
    {
        if(iso[i]->m_ntris)
        {
            m_nvert += iso[i]->m_nvert;
        }
    }

    vector< vector<int> > isocon;
    isocon.resize(niso);
    vector<int> global_to_unique_map;
    global_to_unique_map.resize(m_nvert);

    vector< pair<int,int> > global_to_iso_map;
    global_to_iso_map.resize(m_nvert);

    //Fill vertex array into rtree format
    id2 = 0;
    std::vector<PointPair> inPoints;
    for(i = 0; i < niso; ++i)
    {
        for(id1 = 0; id1 < iso[i]->m_nvert; ++id1)
        {
            inPoints.push_back(PointPair(BPoint( iso[i]->m_x[id1],
                                                 iso[i]->m_y[id1],
                                                 iso[i]->m_z[id1]), id2));
            global_to_unique_map[id2]=id2;
            global_to_iso_map[id2] = make_pair(i,id1);
            id2++;
        }
    }


    if(verbose)
    {
        cout << "\t Process building tree ..." << endl;
    }

    //Build tree
    bgi::rtree<PointPair, bgi::rstar<16> > rtree;
    rtree.insert(inPoints.begin(), inPoints.end());

    //Find neighbours
    int      unique_index = 0;
    int      prog=0;
    for(i = 0; i < m_nvert; ++i)
    {
        if(verbose)
        {
            prog = LibUtilities::PrintProgressbar(i,m_nvert,
                                                  "Nearest verts",prog);
        }

        BPoint queryPoint = inPoints[i].first;


        // check to see if point has been already reset to lower than
        // unique value
        if(global_to_unique_map[i] < unique_index) // do nothing
        {
        }
        else
        {
            // find nearest 10 points within the distance box
            std::vector<PointPair> result;
            rtree.query(bgi::nearest(queryPoint, 10), std::back_inserter(result));

            //see if any values have unique value  already
            set<int> samept;
            int new_index = -1;
            for(id1 = 0; id1 < result.size(); ++id1)
            {
                NekDouble dist = bg::distance(queryPoint, result[id1].first);
                if(dist*dist < NekConstants::kNekZeroTol) // same point
                {
                    id2 = result[id1].second;
                    samept.insert(id2);

                    if(global_to_unique_map[id2] <unique_index)
                    {
                        new_index = global_to_unique_map[id2];
                    }
                }
            }
            if(new_index == -1)
            {
                new_index = unique_index;
                unique_index++;
            }

            // reset all same values to new_index
            global_to_unique_map[i] = new_index;
            for(auto &it : samept)
            {
                global_to_unique_map[it] = new_index;
            }
        }
    }

    if(verbose)
    {
        cout << endl;
    }

    m_nvert = unique_index;

    nelmt = 0;
    // reset m_vid;
    int cnt = 0;
    for(n = 0; n < niso; ++n)
    {
        for(i = 0; i < iso[n]->m_ntris; ++i,nelmt++)
        {
            for(j=0; j < 3;++j)
            {
                m_vid[3*nelmt+j] = global_to_unique_map[iso[n]->m_vid[3*i+j]+cnt];
            }
        }
        cnt += iso[n]->m_nvert;
    }

    m_ntris = nelmt;

    m_x.resize(m_nvert);
    m_y.resize(m_nvert);
    m_z.resize(m_nvert);

    m_fields.resize(iso[0]->m_fields.size());
    for(i = 0; i < iso[0]->m_fields.size(); ++i)
    {
        m_fields[i].resize(m_nvert);
    }

    for(n = 0; n < global_to_unique_map.size(); ++n)
    {

        m_x[global_to_unique_map[n]] = bg::get<0>(inPoints[n].first);
        m_y[global_to_unique_map[n]] = bg::get<1>(inPoints[n].first);
        m_z[global_to_unique_map[n]] = bg::get<2>(inPoints[n].first);

        int isoid = global_to_iso_map[n].first;
        int ptid  = global_to_iso_map[n].second;
        for(j = 0; j < m_fields.size(); ++j)
        {
            m_fields[j][global_to_unique_map[n]] = iso[isoid]->
                m_fields[j][ptid];
        }
    }
}

ResizeFields()

void Nektar::FieldUtils::Iso::ResizeFields ( int  size )

inline

Definition at line 139 of file ProcessIsoContour.h.

References m_fields, m_nvert, m_x, m_y, and m_z.

        {
            if(size > m_x.size()) // add 1000 element to vectors
            {
                m_x.resize(size+100);
                m_y.resize(size+100);
                m_z.resize(size+100);;
                for(int i = 0; i < m_fields.size(); ++i)
                {
                    m_fields[i].resize(size+1000);
                }

            }
            m_nvert = size;
        }

ResizeVId()

void Nektar::FieldUtils::Iso::ResizeVId ( int  nconn )

inline

Definition at line 129 of file ProcessIsoContour.h.

References m_vid.

        {
            m_vid = Array<OneD, int>(nconn);
        }

SeparateRegions()

void Nektar::FieldUtils::Iso::SeparateRegions	(	std::vector< std::shared_ptr< Iso > > &	iso,
		int	minsize,
		bool	verbose
	)

Definition at line 1075 of file ProcessIsoContour.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::PrintProgressbar(), WARNINGL0, and Vmath::Zero().

{
    int i,j,k,id;
    Array<OneD, vector<int> >vertcon(m_nvert);
    list<int> tocheck;

    Array<OneD, bool> viddone(m_nvert,false);

    // make list of connecting tris around each vertex
    for(i = 0; i < m_ntris; ++i)
    {
        for(j = 0; j < 3; ++j)
        {
            vertcon[m_vid[3*i+j]].push_back(i);
        }
    }

    Array<OneD, int> vidregion(m_nvert,-1);

    int nregions = -1;


    // check all points are assigned to a region
    int progcnt  = -1;
    for(k = 0; k < m_nvert; ++k)
    {
        if(verbose)
        {
            progcnt = LibUtilities::PrintProgressbar(k,m_nvert,"Separating regions",progcnt);
        }

        if(vidregion[k] == -1)
        {
            vidregion[k] = ++nregions;

            // find all elmts around this.. vertex  that need to be checked
            for(auto &ipt : vertcon[k])
            {
                for(i = 0; i < 3; ++i)
                {
                    if(vidregion[id = m_vid[3*ipt+i]] == -1)
                    {
                        tocheck.push_back(id);
                        vidregion[id] = nregions;
                    }
                }
            }
            viddone[k] = 1;

            // check all other neighbouring vertices
            while(tocheck.size())
            {
                auto cid = tocheck.begin();
                while(cid != tocheck.end())
                {
                    if(!viddone[*cid])
                    {
                        for(auto &ipt : vertcon[*cid])
                        {
                            for(i = 0; i < 3; ++i)
                            {
                                if(vidregion[id = m_vid[3*ipt+i]] == -1)
                                {
                                    tocheck.push_back(id);
                                    vidregion[id] = nregions;
                                }
                            }
                        }
                        viddone[*cid] = 1;
                        ++cid;
                        tocheck.pop_front();
                    }
                }
            }
        }
    }
    nregions++;


    Array<OneD, int> nvert(nregions,0);
    Array<OneD, int> nelmt(nregions,0);

    // count nverts
    for(i = 0; i < m_nvert; ++i)
    {
        nvert[vidregion[i]] +=1;
    }

    // count nelmts
    for(i = 0; i < m_ntris; ++i)
    {
        nelmt[vidregion[m_vid[3*i]]] +=1;
    }

    Array<OneD, int> vidmap(m_nvert);
    // generate new list of isocontour
    for(int n = 0; n < nregions; ++n)
    {
        if(nelmt[n] > minsize)
        {
            int nfields = m_fields.size();
            IsoSharedPtr iso = MemoryManager<Iso>::AllocateSharedPtr(nfields);

            iso->m_ntris = nelmt[n];
            iso->m_vid = Array<OneD, int>(3*nelmt[n]);

            iso->m_nvert = nvert[n];
            iso->m_x.resize(nvert[n]);
            iso->m_y.resize(nvert[n]);
            iso->m_z.resize(nvert[n]);

            iso->m_fields.resize(nfields);
            for(i = 0; i < nfields; ++i)
            {
                iso->m_fields[i].resize(nvert[n]);
            }


            int cnt = 0;
            // generate vid map;
            Vmath::Zero(m_nvert,vidmap,1);
            for(i = 0; i < m_nvert; ++i)
            {
                if(vidregion[i] == n)
                {
                    vidmap[i] = cnt++;
                }
            }

            cnt = 0;
            for(i = 0; i < m_ntris; ++i)
            {
                if(vidregion[m_vid[3*i]] == n)
                {
                    for(j = 0; j < 3; ++j)
                    {
                        iso->m_vid[3*cnt+j] = vidmap[m_vid[3*i+j]];

                        iso->m_x[vidmap[m_vid[3*i+j]]] = m_x[m_vid[3*i+j]];
                        iso->m_y[vidmap[m_vid[3*i+j]]] = m_y[m_vid[3*i+j]];
                        iso->m_z[vidmap[m_vid[3*i+j]]] = m_z[m_vid[3*i+j]];

                        for(k = 0; k < nfields; ++k)
                        {
                            iso->m_fields[k][vidmap[m_vid[3*i+j]]] = m_fields[k][m_vid[3*i+j]];
                        }
                    }
                    cnt++;
                }
            }

            WARNINGL0(cnt == nelmt[n],"Number of elements do not match");
            sep_iso.push_back(iso);
        }
    }
}

SetFields()

void Nektar::FieldUtils::Iso::SetFields ( const int  loc, const Array< OneD, Array< OneD, NekDouble > > &  intfields, const int  j  )

inline

Definition at line 75 of file ProcessIsoContour.h.

References CG_Iterations::loc, m_fields, m_x, m_y, and m_z.

        {
            m_x[loc] = intfields[0][j];
            m_y[loc] = intfields[1][j];
            m_z[loc] = intfields[2][j];

            for(int i = 0; i < intfields.num_elements()-3; ++i)
            {
                m_fields[i][loc] = intfields[i+3][j];
            }
        }

SetNTris()

void Nektar::FieldUtils::Iso::SetNTris ( int  n )

inline

Definition at line 70 of file ProcessIsoContour.h.

References m_ntris.

        {
            m_ntris = n;
        }

SetNVert()

void Nektar::FieldUtils::Iso::SetNVert ( int  n )

inline

Definition at line 60 of file ProcessIsoContour.h.

References m_nvert.

        {
            m_nvert = n;
        }

SetVId()

void Nektar::FieldUtils::Iso::SetVId ( int  i, int  j  )

inline

Definition at line 134 of file ProcessIsoContour.h.

References m_vid.

        {
            m_vid[i] = j;
        }

SetX()

void Nektar::FieldUtils::Iso::SetX ( int  loc, NekDouble  val  )

inline

Definition at line 94 of file ProcessIsoContour.h.

References CG_Iterations::loc, and m_x.

        {
            m_x[loc] = val;
        }

SetY()

void Nektar::FieldUtils::Iso::SetY ( int  loc, NekDouble  val  )

inline

Definition at line 99 of file ProcessIsoContour.h.

References CG_Iterations::loc, and m_y.

        {
            m_y[loc] = val;
        }

SetZ()

void Nektar::FieldUtils::Iso::SetZ ( int  loc, NekDouble  val  )

inline

Definition at line 104 of file ProcessIsoContour.h.

References CG_Iterations::loc, and m_z.

        {
            m_z[loc] = val;
        }

Smooth()

void Nektar::FieldUtils::Iso::Smooth	(	int	n_iter,
		NekDouble	lambda,
		NekDouble	mu
	)

Definition at line 971 of file ProcessIsoContour.cpp.

{
    int   iter,i,j,k;
    NekDouble del_v[3];
    NekDouble w;
    Array<OneD, NekDouble>  xtemp, ytemp, ztemp;
    vector< vector<int> > adj,vertcon;
    vector< vector<NekDouble > >  wght;

    // determine elements around each vertex
    vertcon.resize(m_nvert);
    for(i = 0; i < m_ntris; ++i)
    {
        for(j = 0; j < 3; ++j)
        {
            vertcon[m_vid[3*i+j]].push_back(i);
        }
    }

    // determine vertices and weights around each vertex
    adj.resize(m_nvert);
    wght.resize(m_nvert);

    for(i =0; i < m_nvert; ++i)
    {
        // loop over surrounding elements
        for(auto &ipt : vertcon[i])
        {
            for(j = 0; j < 3; ++j)
            {
                // make sure not adding own vertex
                if(m_vid[3*ipt+j] != i)
                {
                    // check to see if vertex has already been added
                    for(k = 0; k < adj[i].size(); ++k)
                    {
                        if(adj[i][k] == m_vid[3*ipt+j])
                        {
                            break;
                        }
                    }

                    if(k == adj[i].size())
                    {
                        adj[i].push_back(m_vid[3*ipt+j]);
                    }
                }
            }
        }
    }

    // Currently set weights up as even distribution
    for(i =0; i < m_nvert; ++i)
    {
        w = 1.0/((NekDouble)adj[i].size());
        for(j = 0; j < adj[i].size(); ++j)
        {
            wght[i].push_back(w);
        }
    }


    xtemp = Array<OneD, NekDouble>(m_nvert);
    ytemp = Array<OneD, NekDouble>(m_nvert);
    ztemp = Array<OneD, NekDouble>(m_nvert);

    // smooth each point
    for (iter=0;iter<n_iter;iter++)
    {
        // compute first weighted average
        for(i=0;i< m_nvert;++i)
        {
            del_v[0] = del_v[1] = del_v[2] = 0.0;
            for(j = 0; j < adj[i].size(); ++j)
            {
                del_v[0] =  del_v[0] + (m_x[adj[i][j]]-m_x[i])*wght[i][j];
                del_v[1] =  del_v[1] + (m_y[adj[i][j]]-m_y[i])*wght[i][j];
                del_v[2] =  del_v[2] + (m_z[adj[i][j]]-m_z[i])*wght[i][j];
            }

            xtemp[i] = m_x[i] + del_v[0] * lambda;
            ytemp[i] = m_y[i] + del_v[1] * lambda;
            ztemp[i] = m_z[i] + del_v[2] * lambda;
        }

        // compute second weighted average
        for(i=0;i< m_nvert;++i)
        {
            del_v[0] = del_v[1] = del_v[2] = 0;
            for(j = 0; j < adj[i].size(); ++j)
            {
                del_v[0] =  del_v[0] + (xtemp[adj[i][j]]-xtemp[i])*wght[i][j];
                del_v[1] =  del_v[1] + (ytemp[adj[i][j]]-ytemp[i])*wght[i][j];
                del_v[2] =  del_v[2] + (ztemp[adj[i][j]]-ztemp[i])*wght[i][j];
            }

            m_x[i] = xtemp[i] + del_v[0] * mu;
            m_y[i] = ytemp[i] + del_v[1] * mu;
            m_z[i] = ztemp[i] + del_v[2] * mu;
        }
    }
}

Member Data Documentation

m_condensed

bool Nektar::FieldUtils::Iso::m_condensed

private

Definition at line 175 of file ProcessIsoContour.h.

Referenced by Iso().

m_fields

std::vector<std::vector<NekDouble> > Nektar::FieldUtils::Iso::m_fields

private

Definition at line 181 of file ProcessIsoContour.h.

Referenced by GetFields(), Iso(), ResizeFields(), and SetFields().

m_ntris

int Nektar::FieldUtils::Iso::m_ntris

private

Definition at line 177 of file ProcessIsoContour.h.

Referenced by GetNTris(), and SetNTris().

m_nvert

int Nektar::FieldUtils::Iso::m_nvert

private

Definition at line 176 of file ProcessIsoContour.h.

Referenced by GetNVert(), Iso(), ResizeFields(), and SetNVert().

m_vid

Array<OneD, int> Nektar::FieldUtils::Iso::m_vid

private

Definition at line 182 of file ProcessIsoContour.h.

Referenced by GetVId(), ResizeVId(), and SetVId().

m_x

std::vector<NekDouble> Nektar::FieldUtils::Iso::m_x

private

Definition at line 178 of file ProcessIsoContour.h.

Referenced by GetX(), Iso(), Nektar::FieldUtils::IsoVertex::IsoVertex(), ResizeFields(), SetFields(), and SetX().

m_y

std::vector<NekDouble> Nektar::FieldUtils::Iso::m_y

private

Definition at line 179 of file ProcessIsoContour.h.

Referenced by GetY(), Iso(), Nektar::FieldUtils::IsoVertex::IsoVertex(), ResizeFields(), SetFields(), and SetY().

m_z

std::vector<NekDouble> Nektar::FieldUtils::Iso::m_z

private

Definition at line 180 of file ProcessIsoContour.h.

Referenced by GetZ(), Iso(), Nektar::FieldUtils::IsoVertex::IsoVertex(), ResizeFields(), SetFields(), and SetZ().