#include <ProcessIsoContour.h>

Collaboration diagram for Nektar::FieldUtils::Iso:

Public Member Functions
void	Condense (void)

void	GlobalCondense (vector< boost::shared_ptr< Iso > > &iso, bool verbose)

void	SeparateRegions (vector< boost::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

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 58 of file ProcessIsoContour.h.

Constructor & Destructor Documentation

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

inline

Definition at line 167 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::FieldUtils::Iso::~Iso ( void )

inline

Definition at line 182 of file ProcessIsoContour.h.

183 {

184 }

Member Function Documentation

void Nektar::FieldUtils::Iso::Condense ( void )

Definition at line 686 of file ProcessIsoContour.cpp.

References Nektar::StdRegions::find(), Nektar::iterator, m_condensed, m_fields, Nektar::FieldUtils::IsoVertex::m_fields, Nektar::FieldUtils::IsoVertex::m_id, m_ntris, m_nvert, m_vid, m_x, Nektar::FieldUtils::IsoVertex::m_x, m_y, Nektar::FieldUtils::IsoVertex::m_y, m_z, and Nektar::FieldUtils::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::FieldUtils::Iso::GetFields	(	const int	i,
		const int	j
	)

inline

Definition at line 101 of file ProcessIsoContour.h.

References m_fields.

         {
             return m_fields[i][j];
         }

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

inline

Definition at line 77 of file ProcessIsoContour.h.

References m_ntris.

         {
             return m_ntris;
         }

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

inline

Definition at line 67 of file ProcessIsoContour.h.

References m_nvert.

         {
             return m_nvert;
         }

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

inline

Definition at line 136 of file ProcessIsoContour.h.

References m_vid.

         {
             return m_vid[i];
         }

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

inline

Definition at line 121 of file ProcessIsoContour.h.

References m_x.

         {
             return m_x[loc];
         }

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

inline

Definition at line 126 of file ProcessIsoContour.h.

References m_y.

         {
             return m_y[loc];
         }

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

inline

Definition at line 131 of file ProcessIsoContour.h.

References m_z.

         {
             return m_z[loc];
         }

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

Definition at line 793 of file ProcessIsoContour.cpp.

References Nektar::iterator, m_condensed, m_fields, m_ntris, m_nvert, m_vid, m_x, m_y, m_z, Nektar::LibUtilities::PrintProgressbar(), and Nektar::FieldUtils::SQ_PNT_TOL.

 {
     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 << "Process building tree ..." << endl;
     }
 
     //Build tree
     bgi::rtree<PointPair, bgi::rstar<16> > rtree;
     rtree.insert(inPoints.begin(), inPoints.end());
 
     //Find neipghbours
     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;
             set<int>::iterator it;
             int new_index = -1;
             for(id1 = 0; id1 < result.size(); ++id1)
             {
                 NekDouble dist = bg::distance(queryPoint, result[id1].first);
                 if(dist*dist<SQ_PNT_TOL) // 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(it = samept.begin(); it != samept.end(); ++it)
             {
                 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];
         }
     }
 }

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

inline

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

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

inline

Definition at line 141 of file ProcessIsoContour.h.

References m_vid.

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

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

Definition at line 1094 of file ProcessIsoContour.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::iterator, m_fields, m_ntris, m_nvert, m_vid, m_x, m_y, m_z, Nektar::LibUtilities::PrintProgressbar(), WARNINGL0, and Vmath::Zero().

 {
     int i,j,k,id;
     Array<OneD, vector<int> >vertcon(m_nvert);
     vector<int>::iterator ipt;
     list<int> tocheck;
     list<int>::iterator cid;
 
     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(ipt = vertcon[k].begin(); ipt != vertcon[k].end(); ++ipt)
             {
                 for(i = 0; i < 3; ++i)
                 {
                     if(vidregion[id = m_vid[3*(ipt[0])+i]] == -1)
                     {
                         tocheck.push_back(id);
                         vidregion[id] = nregions;
                     }
                 }
             }
             viddone[k] = 1;
 
             // check all other neighbouring vertices
             while(tocheck.size())
             {
                 cid = tocheck.begin();
                 while(cid != tocheck.end())
                 {
                     if(!viddone[*cid])
                     {
                         for(ipt = vertcon[*cid].begin(); ipt != vertcon[*cid].end(); ++ipt)
                         {
                             for(i = 0; i < 3; ++i)
                             {
                                 if(vidregion[id = m_vid[3*(ipt[0])+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);
         }
     }
 }

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

inline

Definition at line 87 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::FieldUtils::Iso::SetNTris ( int n )

inline

Definition at line 82 of file ProcessIsoContour.h.

References m_ntris.

         {
             m_ntris = n;
         }

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

inline

Definition at line 72 of file ProcessIsoContour.h.

References m_nvert.

         {
             m_nvert = n;
         }

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

inline

Definition at line 146 of file ProcessIsoContour.h.

References m_vid.

         {
             m_vid[i] = j;
         }

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

inline

Definition at line 106 of file ProcessIsoContour.h.

References m_x.

         {
             m_x[loc] = val;
         }

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

inline

Definition at line 111 of file ProcessIsoContour.h.

References m_y.

         {
             m_y[loc] = val;
         }

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

inline

Definition at line 116 of file ProcessIsoContour.h.

References m_z.

         {
             m_z[loc] = val;
         }

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

Definition at line 988 of file ProcessIsoContour.cpp.

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

 {
     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;
     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 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(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(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;
         }
     }
 }