Nektar::Utilities::Iso Class Reference

#include <ProcessIsoContour.h>

Collaboration diagram for Nektar::Utilities::Iso:

Public Member Functions
void	condense (void)
void	globalcondense (vector< boost::shared_ptr< Iso > > &iso)
void	smooth (int n_iter, NekDouble lambda, NekDouble mu)
int	get_nvert (void)
void	set_nvert (int n)
int	get_ntris (void)
void	set_ntris (int n)
void	set_fields (const int loc, const Array< OneD, Array< OneD, NekDouble > > &intfields, const int j)
NekDouble	get_fields (const int i, const int j)
void	set_x (int loc, NekDouble val)
void	set_y (int loc, NekDouble val)
void	set_z (int loc, NekDouble val)
NekDouble	get_x (int loc)
NekDouble	get_y (int loc)
NekDouble	get_z (int loc)
int	get_vid (int i)
void	resize_vid (int nconn)
void	set_vid (int i, int j)
void	resize_fields (int size)
	Iso (int nfields)
	~Iso (void)

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

Detailed Description

Definition at line 47 of file ProcessIsoContour.h.

Constructor & Destructor Documentation

Nektar::Utilities::Iso::Iso ( int  nfields )

inline

Definition at line 153 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 10000 long
        m_x.resize(10000);
        m_y.resize(10000);
        m_z.resize(10000);
        for(int i = 0; i < m_fields.size(); ++i)
        {
            m_fields[i].resize(10000);
        }
    };

Nektar::Utilities::Iso::~Iso ( void  )

inline

Definition at line 168 of file ProcessIsoContour.h.

    {
    }

Member Function Documentation

void Nektar::Utilities::Iso::condense

(

void

)

Definition at line 524 of file ProcessIsoContour.cpp.

References Nektar::StdRegions::find(), Nektar::iterator, m_condensed, m_fields, Nektar::Utilities::IsoVertex::m_fields, Nektar::Utilities::IsoVertex::m_id, m_ntris, m_nvert, m_vid, m_x, Nektar::Utilities::IsoVertex::m_x, m_y, Nektar::Utilities::IsoVertex::m_y, m_z, and Nektar::Utilities::IsoVertex::m_z.

{
    register int i,j,cnt;
    IsoVertex v;
    vector<IsoVertex> vert;
    vector<IsoVertex>::iterator pt;

    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];

            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];
        }
    }
}

NekDouble Nektar::Utilities::Iso::get_fields ( const int  i, const int  j  )

inline

Definition at line 88 of file ProcessIsoContour.h.

References m_fields.

    {
        return m_fields[i][j];
    }

int Nektar::Utilities::Iso::get_ntris ( void  )

inline

Definition at line 64 of file ProcessIsoContour.h.

References m_ntris.

    {
        return m_ntris;
    }

int Nektar::Utilities::Iso::get_nvert ( void  )

inline

Definition at line 54 of file ProcessIsoContour.h.

References m_nvert.

    {
        return m_nvert;
    }

int Nektar::Utilities::Iso::get_vid ( int  i )

inline

Definition at line 122 of file ProcessIsoContour.h.

References m_vid.

    {
        return m_vid[i];
    }

NekDouble Nektar::Utilities::Iso::get_x ( int  loc )

inline

Definition at line 107 of file ProcessIsoContour.h.

References m_x.

    {
        return m_x[loc];
    }

NekDouble Nektar::Utilities::Iso::get_y ( int  loc )

inline

Definition at line 112 of file ProcessIsoContour.h.

References m_y.

    {
        return m_y[loc];
    }

NekDouble Nektar::Utilities::Iso::get_z ( int  loc )

inline

Definition at line 117 of file ProcessIsoContour.h.

References m_z.

    {
        return m_z[loc];
    }

void Nektar::Utilities::Iso::globalcondense

(

vector< boost::shared_ptr< Iso > > &

iso

)

Definition at line 659 of file ProcessIsoContour.cpp.

References m_condensed, m_fields, m_ntris, m_nvert, m_vid, m_x, m_y, m_z, and Nektar::Utilities::same().

{
    int    i,j,n;
    int    nvert,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);

    // identify which iso are connected by at least one point;
    // find min x,y,z and max x,y,z and see if overlap to select
    // which zones should be connected
    {
        vector<Array<OneD, NekDouble> > sph(niso);
        Array<OneD, NekDouble> rng(6);
        for(i = 0; i < niso; ++i)
        {
            sph[i] = Array<OneD, NekDouble>(4);

            // find max and min of isocontour
            rng[0] = rng[3] = iso[i]->m_x[0];
            rng[1] = rng[4] = iso[i]->m_x[1];
            rng[2] = rng[5] = iso[i]->m_x[2];

            for(id1 = 1; id1 < iso[i]->m_nvert;++id1)
            {
                rng[0] = min(rng[0],iso[i]->m_x[i]);
                rng[1] = min(rng[1],iso[i]->m_y[i]);
                rng[2] = min(rng[2],iso[i]->m_z[i]);

                rng[3] = max(rng[3],iso[i]->m_x[i]);
                rng[4] = max(rng[4],iso[i]->m_y[i]);
                rng[5] = max(rng[5],iso[i]->m_z[i]);
            }

            // centroid
            sph[i][0] = (rng[3]+rng[0])/2.0;
            sph[i][1] = (rng[4]+rng[1])/2.0;
            sph[i][2] = (rng[5]+rng[2])/2.0;

            // radius;
            sph[i][3] = sqrt((rng[3]-rng[0])*(rng[3]-rng[0]) +
                             (rng[4]-rng[1])*(rng[4]-rng[1]) +
                             (rng[5]-rng[2])*(rng[5]-rng[2]));
        }

        for(i = 0; i < niso; ++i)
        {
            for(j = i+1; j < niso; ++j)
            {
                NekDouble diff=sqrt((sph[i][0]-sph[j][0])*(sph[i][0]-sph[j][0])+
                          (sph[i][1]-sph[j][1])*(sph[i][1]-sph[j][1])+
                          (sph[i][2]-sph[j][2])*(sph[i][2]-sph[j][2]));

                // if centroid is closer than added radii
                if(diff < sph[i][3] + sph[j][3])
                {
                    isocon[i].push_back(j);
                }
            }
        }
    }


    for(i = 0; i < niso; ++i)
    {
        vidmap[i] = Array<OneD, int>(iso[i]->m_nvert,-1);
    }
    nvert = 0;
    // identify which vertices are connected to tolerance
    cout << "Matching Vertices [" <<flush;
    int cnt = 0;
    for(i = 0; i < niso; ++i)
    {
        for(n = 0; n < isocon[i].size(); ++n)
        {
            int con = isocon[i][n];
            for(id1 = 0; id1 < iso[i]->m_nvert;++id1)
            {
                for(id2 = 0; id2 < iso[con]->m_nvert;++id2,++cnt)
                {
                    if(cnt%1000000 == 0)
                    {
                        cout <<"." <<flush;
                    }

                    //if((vidmap[con][id2] != -1)&&(vidmap[i][id1] != -1))
                    {
                        if(same(iso[i]->m_x[id1],  iso[i]->m_y[id1],
                                iso[i]->m_z[id1],  iso[con]->m_x[id2],
                                iso[con]->m_y[id2],iso[con]->m_z[id2]))
                        {
                            if((vidmap[i][id1] == -1) &&
                               (vidmap[con][id2] != -1))
                            {
                                vidmap[i][id1] = vidmap[con][id2];
                            }
                            else if((vidmap[con][id2] == -1) &&
                                    (vidmap[i][id1] != -1))
                            {
                                vidmap[con][id2] = vidmap[i][id1];
                            }
                            else if((vidmap[con][id2] == -1) &&
                                    (vidmap[i][id1] == -1))
                            {
                                vidmap[i][id1] = vidmap[con][id2] = nvert++;
                            }
                        }
                    }
                }
            }
        }

        for(id1 = 0; id1 < iso[i]->m_nvert;++id1)
        {
            if(vidmap[i][id1] == -1)
            {
                vidmap[i][id1] = nvert++;
            }
        }
    }
    cout <<"]"<<endl;
    m_nvert = nvert;

    nelmt = 0;
    // reset m_vid;
    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] = vidmap[n][iso[n]->m_vid[3*i+j]];
            }
        }
    }

    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);
    }


    // reset coordinate and fields.
    for(n = 0; n < niso; ++n)
    {
        for(i = 0; i < iso[n]->m_nvert; ++i)
        {
            m_x[vidmap[n][i]] = iso[n]->m_x[i];
            m_y[vidmap[n][i]] = iso[n]->m_y[i];
            m_z[vidmap[n][i]] = iso[n]->m_z[i];

            for(j = 0; j < m_fields.size(); ++j)
            {
                m_fields[j][vidmap[n][i]] = iso[n]->m_fields[j][i];
            }
        }
    }
}

void Nektar::Utilities::Iso::resize_fields ( int  size )

inline

Definition at line 137 of file ProcessIsoContour.h.

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

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

        }
        m_nvert = size;
    }

void Nektar::Utilities::Iso::resize_vid ( int  nconn )

inline

Definition at line 127 of file ProcessIsoContour.h.

References m_vid.

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

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

inline

Definition at line 74 of file ProcessIsoContour.h.

References 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];
        }
    }

void Nektar::Utilities::Iso::set_ntris ( int  n )

inline

Definition at line 69 of file ProcessIsoContour.h.

References m_ntris.

    {
        m_ntris = n;
    }

void Nektar::Utilities::Iso::set_nvert ( int  n )

inline

Definition at line 59 of file ProcessIsoContour.h.

References m_nvert.

    {
        m_nvert = n;
    }

void Nektar::Utilities::Iso::set_vid ( int  i, int  j  )

inline

Definition at line 132 of file ProcessIsoContour.h.

References m_vid.

    {
        m_vid[i] = j;
    }

void Nektar::Utilities::Iso::set_x ( int  loc, NekDouble  val  )

inline

Definition at line 92 of file ProcessIsoContour.h.

References m_x.

    {
        m_x[loc] = val;
    }

void Nektar::Utilities::Iso::set_y ( int  loc, NekDouble  val  )

inline

Definition at line 97 of file ProcessIsoContour.h.

References m_y.

    {
        m_y[loc] = val;
    }

void Nektar::Utilities::Iso::set_z ( int  loc, NekDouble  val  )

inline

Definition at line 102 of file ProcessIsoContour.h.

References m_z.

    {
        m_z[loc] = val;
    }

void Nektar::Utilities::Iso::smooth	(	int	n_iter,
		NekDouble	lambda,
		NekDouble	mu
	)

Definition at line 853 of file ProcessIsoContour.cpp.

References Nektar::iterator, m_ntris, m_nvert, m_vid, m_x, m_y, and m_z.

{
    int   iter,i,j;
    NekDouble del_v[3];
    NekDouble w;
    Array<OneD, NekDouble>  xtemp, ytemp, ztemp;
    vector< vector<int> > adj,vertcon;
    vector<int>::iterator iad;
    vector<int>::iterator ipt;

    // 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 around each vertex
    adj.resize(m_nvert);

    for(i =0; i < m_nvert; ++i)
    {
        for(ipt = vertcon[i].begin(); ipt != vertcon[i].end(); ++ipt)
        {
            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(iad = adj[i].begin(); iad != adj[i].end();++iad)
                    {
                        if(*iad == (m_vid[3*(*ipt)+j])) break;
                    }

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

    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)
        {
            w = 1.0/(NekDouble)(adj[i].size());

            del_v[0] = del_v[1] = del_v[2] = 0.0;

            for(iad = adj[i].begin(); iad != adj[i].end(); ++iad)
            {
                del_v[0] =  del_v[0] + (m_x[*iad]-m_x[i])*w;
                del_v[1] =  del_v[1] + (m_y[*iad]-m_y[i])*w;
                del_v[2] =  del_v[2] + (m_z[*iad]-m_z[i])*w;
            }

            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)
        {

            w = 1.0/(NekDouble)(adj[i].size());
            del_v[0] = del_v[1] = del_v[2] = 0;

            for(iad = adj[i].begin(); iad != adj[i].end(); ++iad)
            {
                del_v[0] =  del_v[0] + (m_x[*iad]-m_x[i])*w;
                del_v[1] =  del_v[1] + (m_y[*iad]-m_y[i])*w;
                del_v[2] =  del_v[2] + (m_z[*iad]-m_z[i])*w;
            }

            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

bool Nektar::Utilities::Iso::m_condensed

private

Definition at line 173 of file ProcessIsoContour.h.

Referenced by condense(), globalcondense(), and Iso().

vector<vector<NekDouble> > Nektar::Utilities::Iso::m_fields

private

Definition at line 179 of file ProcessIsoContour.h.

Referenced by condense(), get_fields(), globalcondense(), Iso(), resize_fields(), and set_fields().

int Nektar::Utilities::Iso::m_ntris

private

Definition at line 175 of file ProcessIsoContour.h.

Referenced by condense(), get_ntris(), globalcondense(), set_ntris(), and smooth().

int Nektar::Utilities::Iso::m_nvert

private

Definition at line 174 of file ProcessIsoContour.h.

Referenced by condense(), get_nvert(), globalcondense(), Iso(), resize_fields(), set_nvert(), and smooth().

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

private

Definition at line 180 of file ProcessIsoContour.h.

Referenced by condense(), get_vid(), globalcondense(), resize_vid(), set_vid(), and smooth().

vector<NekDouble> Nektar::Utilities::Iso::m_x

private

Definition at line 176 of file ProcessIsoContour.h.

Referenced by condense(), get_x(), globalcondense(), Iso(), resize_fields(), set_fields(), set_x(), and smooth().

vector<NekDouble> Nektar::Utilities::Iso::m_y

private

Definition at line 177 of file ProcessIsoContour.h.

Referenced by condense(), get_y(), globalcondense(), Iso(), resize_fields(), set_fields(), set_y(), and smooth().

vector<NekDouble> Nektar::Utilities::Iso::m_z

private

Definition at line 178 of file ProcessIsoContour.h.

Referenced by condense(), get_z(), globalcondense(), Iso(), resize_fields(), set_fields(), set_z(), and smooth().