#include <ProcessIsoContour.h>

Collaboration diagram for Nektar::Utilities::Iso:

Public Member Functions
void	condense (void)

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

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

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

172 {

173 }

Member Function Documentation

void Nektar::Utilities::Iso::condense ( void )

Definition at line 642 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, Nektar::Utilities::IsoVertex::m_z, and Nektar::NekMeshUtils::vert.

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

References m_fields.

         {
             return m_fields[i][j];
         }

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

inline

Definition at line 66 of file ProcessIsoContour.h.

References m_ntris.

         {
             return m_ntris;
         }

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

inline

Definition at line 56 of file ProcessIsoContour.h.

References m_nvert.

         {
             return m_nvert;
         }

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

inline

Definition at line 125 of file ProcessIsoContour.h.

References m_vid.

         {
             return m_vid[i];
         }

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

inline

Definition at line 110 of file ProcessIsoContour.h.

References m_x.

         {
             return m_x[loc];
         }

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

inline

Definition at line 115 of file ProcessIsoContour.h.

References m_y.

         {
             return m_y[loc];
         }

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

inline

Definition at line 120 of file ProcessIsoContour.h.

References m_z.

         {
             return m_z[loc];
         }

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

Definition at line 777 of file ProcessIsoContour.cpp.

References m_condensed, m_fields, m_ntris, m_nvert, m_vid, m_x, m_y, m_z, Nektar::LibUtilities::PrintProgressbar(), 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; 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;
     int cnt = 0;
     // count up amount of checking to be done
     NekDouble totiso = 0; 
     for(i = 0; i < niso; ++i)
     {
         totiso += isocon[i].size();
     }
 
 
     if(verbose)
     {
         cout << "Progress Bar totiso: " << totiso << endl;
     }
     for(i = 0; i < niso; ++i)
     {
         for(n = 0; n < isocon[i].size(); ++n, ++cnt)
         {
             
             if(verbose && totiso >= 40)
             {
                 LibUtilities::PrintProgressbar(cnt,totiso,"Condensing verts");
             }
 
             int con = isocon[i][n];
             for(id1 = 0; id1 < iso[i]->m_nvert; ++id1)
             {
 
                 if(verbose && totiso < 40)
                 {
                     LibUtilities::PrintProgressbar(id1,iso[i]->m_nvert,"isocon");
                 }
 
                 int start  = 0; 
                 if(con == i)
                 {
                     start = id1+1;
                 }
                 for(id2 = start; id2 < iso[con]->m_nvert; ++id2)
                 {
                     
                     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++;
             }
         }
     }
     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];
             }
         }
     }
     cout << endl;
 }

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

inline

Definition at line 140 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::Utilities::Iso::resize_vid ( int nconn )

inline

Definition at line 130 of file ProcessIsoContour.h.

References m_vid.

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

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

Definition at line 1089 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
         for(k = 0; k < m_nvert; ++k)
         {
             if(verbose)
             {
                 LibUtilities::PrintProgressbar(k,m_nvert,"Separating regions");
             }
 
             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::Utilities::Iso::set_fields	(	const int	loc,
		const Array< OneD, Array< OneD, NekDouble > > &	intfields,
		const int	j
	)

inline

Definition at line 76 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 71 of file ProcessIsoContour.h.

References m_ntris.

         {
             m_ntris = n;
         }

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

inline

Definition at line 61 of file ProcessIsoContour.h.

References m_nvert.

         {
             m_nvert = n;
         }

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

inline

Definition at line 135 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 95 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 100 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 105 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 994 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;
         }
     }
 }