Commit e77133a4 authored by Julia's avatar Julia

stufff

parent 77ee6956
......@@ -2,52 +2,116 @@
#define FILTERTYPE_H
namespace Nektar
{
namespace LSIAC
{
namespace SIACUtilities
{
//! This enum specifies the kind of filter to be used by the object.
/*! The user can specify any type of smoothing or derivative filter.
* All the properties of memeber depends on this enum. This enum cannot be changed after the object has been initialized.
* A new enum needs to be added to the list to create any new type of filter.
*/
enum FilterType
{
eNONE, /**< enum value 1 */
eSYM_2kp1, /**< enum value 2, Will apply filter only when symmetric filter is possible. */
eSYM_4kp1,
eSYM,
eSYM_2kp1_1SIDED_2kp1, /**< enum value 3, Simplest form of SIAC and oneSided filter. */
eSYM_2kp1_1SIDED_2kp2, /**< enum value 4, Xli 2k+1 CentralBspl and 1 GeneralBSlp only at borders appropriately*/
eSYM_2kp1_1SIDED_XLI,
eSYM_2kp1_1SIDED_4kp1, /**< enum value 5, SRV's 4k+1 CentralBspl at boundaries.*/
eSYM_2kp1_1SIDED_SRV, /**< enum value 5, SRV's 4k+1 CentralBspl at boundaries.*/
eSYM_1SIDED_BASIC,
eSYM_1SIDED_XLI,
eSYM_1SIDED,
eSYM_DER_2kp1,
eSYM_DER,
eSYM_DER_1SIDED,
eSYM_DER_2kp1_1SIDED_2kp1, /**< enum value 6, Simple Derivative filter*/
eSYM_DER_2kp1_1SIDED_2kp2, /**< enum value 7, SRV's Derivative filter */
eSYM_DER_2kp1_1SIDED_4kp1, /**< enum value 8, Xli's Derivative fitler */
eSYM_DER_1SIDED_SRV,
eSYM_DER_1SIDED_XLi,
};
struct SMOOTHIE_Status {
int count;
int cancelled;
int SOURCE;
int TAG;
int ERROR;
};
}
}
}
{
namespace LSIAC
{
namespace SIACUtilities
{
//! This enum specifies the kind of filter to be used by the object.
/*! The user can specify any type of smoothing or derivative filter.
* All the properties of memeber depends on this enum. This enum cannot be changed after the object has been initialized.
* A new enum needs to be added to the list to create any new type of filter.
*/
enum FilterType
{
eNONE, /// 0
ePURE,
ePURE_2kp1,
eSTD, /// 1
eSTD_1SIDED_XLI, /// 2
eSTD_2kp1, /// 3
eSTD_2kp1_1SIDED_XLI,/// 4
eSTD_2kp1_1SIDED_4kp1, /// 5
eSTD_2kp1_1SIDED_SRV, /// 6
eDER,/// 7
eDER_1SIDED, ///8
eDER_1SIDED_XLI, /// 9
eDER_2kp1, /// 10
eDER_2kp1_1SIDED, /// 11
eDER_2kp1_1SIDED_XLI, /// 12
eDER_2kp1_1SIDED_4kp1, /// 13
eDER_2kp1_1SIDED_SRV, /// 14
eDER_LOWER_ORDER,
eDER_2kp1_LOWER_ORDER
};
inline void write_filter_name(FilterType filter,string &name)
{
switch(filter)
{
case eNONE:
name="NONE";
return;
case ePURE:
name="PURE";
return;
case ePURE_2kp1:
name="PURE_2kp1";
return;
case eSTD:
name="STD";
return;
case eSTD_1SIDED_XLI:
name="STD_1SIDED_XLI";
return;
case eSTD_2kp1:
name="STD_2kp1";
return;
case eSTD_2kp1_1SIDED_XLI:
name="STD_2kp1_1SIDED_XLI";
return;
case eSTD_2kp1_1SIDED_4kp1:
name="STD_2kp1_1SIDED_4kp1";
return;
case eSTD_2kp1_1SIDED_SRV:
name="STD_2kp1_1SIDED_SRV";
return;
case eDER:
name="DER";
return;
case eDER_1SIDED:
name="DER_1SIDED";
return;
case eDER_1SIDED_XLI:
name="DER_1SIDED_XLI";
return;
case eDER_2kp1:
name="DER_2kp1";
return;
case eDER_2kp1_1SIDED:
name="DER_2kp1_1SIDED";
return;
case eDER_2kp1_1SIDED_XLI:
name="DER_2kp1_1SIDED_XLI";
return;
case eDER_2kp1_1SIDED_4kp1:
name="DER_2kp1_1SIDED_4kp1";
return;
case eDER_2kp1_1SIDED_SRV:
name="DER_2kp1_1SIDED_SRV";
return;
case eDER_LOWER_ORDER:
name="DER_LOWER_ORDER";
return;
case eDER_2kp1_LOWER_ORDER:
name="DER_2kp1_LOWER_ORDER";
return;
}
}
struct SMOOTHIE_Status {
int count;
int cancelled;
int SOURCE;
int TAG;
int ERROR;
};
}
}
}
#endif
......@@ -7,7 +7,6 @@ namespace Nektar
bool HandleNekMesh::EvaluateAt(const Array<OneD,Array<OneD,NekDouble> >point,
const int gID, const int eID, Array<OneD,NekDouble> &values,int varNum)
{
if (!m_expansions[varNum]->GetPhysState())
{
m_expansions[varNum]->BwdTrans(m_expansions[varNum]->GetCoeffs(),m_expansions[varNum]->UpdatePhys());
......@@ -22,9 +21,51 @@ namespace Nektar
return true;
}
bool HandleNekMesh::GetListOfGIDs(const Array<OneD,NekDouble> Pos, Array<OneD,NekDouble> &direction, vector<double> t_breaks, Array<OneD,int> &t_GIDs,
Array<OneD, int> &t_EIDs)
{
t_GIDs=Array<OneD,int>(t_breaks.size());
t_EIDs=Array<OneD,int>(t_breaks.size());
Array<OneD,NekDouble> locCoord(m_meshdim,0.0);
for( int i =0; i < t_breaks.size()-1;i++)
{
t_GIDs[i] = -1; t_EIDs[i] = -1;
NekDouble t_break = (t_breaks[i] + t_breaks[i+1])/2.0;
for(int d=0;d<m_meshdim;++d)
locCoord[d] = Pos[d]+t_break*direction[d];
if (m_useRTree)
{
t_GIDs[i] = GetExpansionIndexUsingRTree(locCoord);
//t_GIDs[i] = m_expansions[0]->GetExpIndex(locCoord,TOLERENCE);
}else{
t_GIDs[i] = m_expansions[0]->GetExpIndex(locCoord,TOLERENCE);
}
t_EIDs[i] = t_GIDs[i];
if (t_GIDs[i] <0)
{
// cout<<" LOCAL COORD "<<locCoord[0]<<" "<<locCoord[1]<<endl;
for(int d=0;d<m_meshdim;++d)
locCoord[d] -= TOLERENCE*direction[d];
if (m_useRTree)
{
t_GIDs[i] = GetExpansionIndexUsingRTree(locCoord);
//t_GIDs[i] = m_expansions[0]->GetExpIndex(locCoord,TOLERENCE);
}
else
t_GIDs[i] = m_expansions[0]->GetExpIndex(locCoord,TOLERENCE);
if (t_GIDs[i] <0)
{
assert( t_GIDs[i] >=0 && "Will fail down the line");
}
}
}
return true;
}
void HandleNekMesh::IntersectWithEdges2(vector<double> &tvalT,Array<OneD,NekDouble> x_bar, Array<OneD,NekDouble> direction,
double t1,double t2)
bool HandleNekMesh::IntersectWithEdges(vector<double> &tvalT,Array<OneD,NekDouble> x_bar, Array<OneD,NekDouble> direction,
double t1,double t2)
{
tvalT.clear();
Array<OneD,NekDouble> p1(m_meshdim), p2(m_meshdim),intersectionPoint(m_meshdim);
......@@ -40,7 +81,7 @@ namespace Nektar
// pick a direction not zero.
for (int i =0; i<m_meshdim; i++)
{
if (std::abs(direction[i]) >TOLERENCE)
if (std::abs(direction[i]) >PTTOL)
{
dirID = i;
break;
......@@ -54,9 +95,13 @@ namespace Nektar
int element_id=m_expansions[0]->GetExpIndex(intersectionPoint,NekConstants::kGeomFactorsTol);
do
{
if(it >itMax)
assert( it > itMax && " I THINK I GOT STUCK IN INTERESECTION SCAN INSIDE INTERESECTWITHEDGES2");
if(it >itMax){
//cout<<" HOLA "<<endl;
assert( it > itMax && " I THINK I GOT STUCK IN INTERESECTION SCAN INSIDE INTERESECTWITHEDGES2");
return false;}
++it;
LocalRegions::ExpansionSharedPtr e = m_expansions[0]->GetExp(element_id); // pointer at the element
// cout<<" ID "<<element_id<<endl;
LocalRegions::Expansion2DSharedPtr e2d = boost::dynamic_pointer_cast<LocalRegions::Expansion2D> (e); // 2D expansion for edges
for(int edgeCount = 0; edgeCount < e2d->GetNedges() ; ++edgeCount)
{
......@@ -72,12 +117,12 @@ namespace Nektar
Array<OneD,NekDouble> v0coords(m_meshdim),v1coords(m_meshdim);
vert0->GetCoords(v0coords);
vert1->GetCoords(v1coords);
// cout.precision(16);
// cout<<"\n\n COMPARING SEGMENT "<<p1[0]<<" "<<p1[1]<<" WITH "<<p2[0]<<" "<<p2[1]<<endl;
// cout<<" WITH MESH VERTEX SEGMENT "<<v0coords[0]<<" "<<v0coords[1]<<" WITH "<<v1coords[0]<<" "<<v1coords[1]<<endl;
// cout.precision(16);
// cout<<"\n\n COMPARING SEGMENT "<<p1[0]<<" "<<p1[1]<<" WITH "<<p2[0]<<" "<<p2[1]<<endl;
// cout<<" WITH MESH VERTEX SEGMENT "<<v0coords[0]<<" "<<v0coords[1]<<" WITH "<<v1coords[0]<<" "<<v1coords[1]<<endl;
if(intersect(p1,p2,v0coords,v1coords,intersectionPoint,i2))
{
// cout<<" INTERSECTION POINT "<<intersectionPoint[0]<<" "<<intersectionPoint[1]<<endl;
// cout<<" INTERSECTION POINT "<<intersectionPoint[0]<<" "<<intersectionPoint[1]<<endl;
edgeCount=e2d->GetNedges();
NekDouble t =(intersectionPoint[dirID]-x_bar[dirID])/direction[dirID];
tvalT.push_back(t);
......@@ -87,7 +132,7 @@ namespace Nektar
}
for(int i=0;i<m_meshdim;++i)
{
intersectionPoint[i]+=TOLERENCE*direction[i];
intersectionPoint[i]+=PTTOL*direction[i];
}
if(i2.size() > 0)
{
......@@ -95,30 +140,33 @@ namespace Nektar
tvalT.push_back(t);
for(int i=0;i<m_meshdim;++i)
{
intersectionPoint[i]=i2[i]+direction[i]*(TOLERENCE);
intersectionPoint[i]=i2[i]+direction[i]*(PTTOL);
}
// cout<<" ANOTHER !"<<intersectionPoint[0]<<"\t"<<intersectionPoint[1]<<endl;
// cout<<" ANOTHER !"<<intersectionPoint[0]<<"\t"<<intersectionPoint[1]<<endl;
}
element_id=m_expansions[0]->GetExpIndex(intersectionPoint,NekConstants::kGeomFactorsTol);
if(element_id == -1 )
element_id = lastPointId;
if(fabs(t-t2)<TTOL)
element_id=lastPointId;
// cout<<" INTERSECTION POINT "<<intersectionPoint[0]<<" "<<intersectionPoint[1]<<endl;
// cout<<" INTERSECTION POINT "<<intersectionPoint[0]<<" "<<intersectionPoint[1]<<endl;
}
}
}
}
while(lastPointId != element_id);
tvalT.push_back(t1);
tvalT.push_back(t2);
std::sort(tvalT.begin(),tvalT.end());
std::vector<NekDouble>::iterator ittt;
ittt = std::unique(tvalT.begin(),tvalT.end(),this->compare2NekDoublesH);
tvalT.resize( std::distance(tvalT.begin(),ittt));
return true;
}
void HandleNekMesh::IntersectWithEdges ( const SpatialDomains::SegGeomMap &segMap, const SpatialDomains::PointGeomMap &pointMap,
const Array<OneD,NekDouble> &direction, const Array<OneD,NekDouble> &point, const NekDouble t1, const NekDouble t2,
vector<NekDouble> &tvalT)
void HandleNekMesh::IntersectWithAllEdges ( const SpatialDomains::SegGeomMap &segMap, const SpatialDomains::PointGeomMap &pointMap,
const Array<OneD,NekDouble> &direction, const Array<OneD,NekDouble> &point, const NekDouble t1, const NekDouble t2,
vector<NekDouble> &tvalT)
{
tvalT.clear();
Array<OneD,NekDouble> p1(m_meshdim), p2(m_meshdim),intersectionPoint (m_meshdim,0.);
......@@ -159,6 +207,8 @@ namespace Nektar
}
}
// sort and keep them unique
tvalT.push_back(t1);
tvalT.push_back(t2);
std::sort(tvalT.begin(),tvalT.end());
std::vector<NekDouble>::iterator ittt;
ittt = std::unique(tvalT.begin(),tvalT.end(),this->compare2NekDoublesH);
......@@ -166,15 +216,109 @@ namespace Nektar
}
bool HandleNekMesh::threeD_intersect( Array<OneD,NekDouble> &p1, Array<OneD,NekDouble> &p2,
Array<OneD,NekDouble> &q1, Array<OneD,NekDouble> &q2,
Array<OneD,NekDouble>&i1, vector<NekDouble>&i2 )
{
// Assuming all points are 3D.
// 1. compute r and s ; r= P2-P1; s = Q2-Q1;
i2.clear();
NekDouble t0,t1;
Array<OneD,NekDouble> r(threeD,0.), s(threeD,0.) , pMq(threeD,0.);
sub_Math(p2,p1,r,threeD); sub_Math(q2,q1,s,threeD); sub_Math(p1,q1,pMq,threeD);
Array<OneD,NekDouble> rCs = threeDcross_Math(r,s);
Array<OneD,NekDouble> pMq_Cr = threeDcross_Math(pMq,r);
if( (norm2_Math(rCs,threeD) < TOLERENCE) && (norm2_Math(pMq_Cr,threeD) < TOLERENCE) ) // colinear segments. We will get two intersection points.
{
t0 = -1.0*dot_Math(pMq,r,threeD)/dot_Math(r,r,threeD);
Array<OneD,NekDouble> q2Mp(threeD,0.);
sub_Math(q2,p1,q2Mp,threeD);
t1 = dot_Math(q2Mp,r,threeD)/dot_Math(r,r,threeD);
if (t0 > t1)
{
NekDouble temp = t0;
t0 = t1;
t1 = temp;
}
if (t0 <-TOLERENCE && t1 <-TOLERENCE) // no intersection
return false;
else if (t0<-TOLERENCE && t1<1.0+TOLERENCE) // p0 is intersection point and t1.
{
i1=p1;
for(int d=0;d<p1.num_elements();++d)
i2.push_back( p1[d] + t1*r[d] );
return true;
}
else if (t0<-TOLERENCE && t1>1.0+TOLERENCE) // p0 and p1 are intersection point
{
i1=p1;
for(int d=0;d<threeD;++d)
i2.push_back( p2[d] );
return true;
}
else if ( t0< 1.0+TOLERENCE && t1 > 1.0+TOLERENCE) // t0 is intersectino point and endpoint
{
for(int d=0;d<p1.num_elements();++d)
{
i1[d]=p1[d] + t0*r[d];
i2.push_back( p2[d] );
}
return true;
}
else
return false;
}
if( (norm2_Math(rCs,threeD) < TOLERENCE) && (norm2_Math(pMq_Cr,threeD) > TOLERENCE) ) // Skew lines
{
return false;
}
if( norm2_Math(rCs,threeD) > TOLERENCE) // Proper intersection: get single point
{
Array<OneD,NekDouble> qMp(threeD);
sub_Math(q1,p1,qMp,threeD);
Array<OneD,NekDouble> qMp_Cs = threeDcross_Math( qMp,s );
Array<OneD,NekDouble> qMp_Cr = threeDcross_Math( qMp,r );
rCs= threeDcross_Math(r,s);
NekDouble t = std::sqrt( norm2_Math(qMp_Cs,threeD)/norm2_Math(rCs,threeD) );
NekDouble terr = (std::abs(qMp_Cs[0] -t*rCs[0]) +
std::abs(qMp_Cs[1]-t*rCs[1]) + std::abs(qMp_Cs[2] - t*rCs[2]) );
NekDouble u = std::sqrt( norm2_Math(qMp_Cr,threeD)/norm2_Math(rCs,threeD) );
NekDouble uerr = (std::abs(qMp_Cr[0] -u*rCs[0]) +
std::abs(qMp_Cr[1]-u*rCs[1]) + std::abs(qMp_Cr[2] - u*rCs[2]) );
if ( (t > -TOLERENCE && t < 1+TOLERENCE) &&(u> -TOLERENCE && u < 1+TOLERENCE) &&
(terr<TOLERENCE) && (uerr<TOLERENCE) )
{
for(int d=0;d<p1.num_elements();++d)
i1[d]=p1[d] + t*r[d];
return true;
}
else{
return false;
}
}
return true;
}
bool HandleNekMesh::filterFits(const Array<OneD,NekDouble> point, Array<OneD,NekDouble> &direction,Array<OneD,Array<OneD,NekDouble> >&rotation_axis,int Ntries,NekDouble &tmin,NekDouble &tmax)
bool HandleNekMesh::filterFits(const Array<OneD,NekDouble> point, Array<OneD,NekDouble> &direction,NekDouble &tmin_sym,NekDouble &tmax_sym,NekDouble &tmin_one,NekDouble &tmax_one)
{
Array<OneD,Array<OneD,NekDouble> > axis(1);
axis[0]=direction;
return filterFits(point,direction,axis,1,tmin_sym,tmax_sym,tmin_one,tmax_one) ;
}
bool HandleNekMesh::filterFits(const Array<OneD,NekDouble> point, Array<OneD,NekDouble> &direction,Array<OneD,Array<OneD,NekDouble> >&rotation_axis,int Ntries,NekDouble &tmin_sym,NekDouble &tmax_sym,NekDouble &tmin_one,NekDouble &tmax_one)
{
bool support_fits=false;
int good_rotation=0;
double tmin_aux,tmax_aux;
for(int n=0;n<Ntries;++n)
{
support_fits=checkSupport(point,rotation_axis[n],tmin,tmax,false);
// cout<<" ROTATION AXIS "<<rotation_axis[n][0]<<" "<<rotation_axis[n][1]<<" "<<rotation_axis[n][2]<<endl;
support_fits=checkSupport(point,rotation_axis[n],tmin_sym,tmax_sym,false);
if(support_fits)
{
good_rotation = n;
......@@ -185,8 +329,8 @@ namespace Nektar
{
for(int n=0;n<Ntries;++n)
{
tmin_aux = tmin;
tmax_aux = tmax;
tmin_aux = tmin_one;
tmax_aux = tmax_one;
support_fits=checkSupport(point,rotation_axis[n],tmin_aux,tmax_aux,true);
if(support_fits)
{
......@@ -194,8 +338,8 @@ namespace Nektar
break;
}
}
tmin = tmin_aux;
tmax = tmax_aux;
tmin_one= tmin_aux;
tmax_one = tmax_aux;
}
if(support_fits)
{
......@@ -207,22 +351,29 @@ namespace Nektar
return false;
}
bool HandleNekMesh::checkSupport(const Array<OneD,NekDouble> point, Array<OneD,NekDouble> &direction, NekDouble &tmin,NekDouble &tmax,bool try_shifting)
{
Array<OneD,NekDouble> left_limit(m_meshdim),right_limit(m_meshdim);
cout<<" dimension "<<m_meshdim<<endl;
cout<<" TMIN "<<tmin<<" TMAX "<<tmax<<endl;
if(fabs(tmin-tmax)<TTOL)
return false;
// cout<<" COORDINATES TO TEST "<<endl;
for(int i=0;i<m_meshdim;++i)
{
left_limit[i]=point[i]+tmin*direction[i];
right_limit[i]=point[i]+tmax*direction[i];
cout<<" LEFT "<<left_limit[i]<<" RIGHT "<<right_limit[i]<<"\t DIR "<<direction[i]<<endl;
// cout<<" LEFT "<<left_limit[i]<<"\t right "<<right_limit[i]<<endl;
}
int lId,rId;
if (m_useRTree)
{
lId = GetExpansionIndexUsingRTree(left_limit);
rId = GetExpansionIndexUsingRTree( right_limit);
}
else
{
lId=m_expansions[0]->GetExpIndex(left_limit,NekConstants::kGeomFactorsTol);
rId=m_expansions[0]->GetExpIndex(right_limit,NekConstants::kGeomFactorsTol);
}
int lId=m_expansions[0]->GetExpIndex(left_limit,NekConstants::kGeomFactorsTol);
int rId=m_expansions[0]->GetExpIndex(right_limit,NekConstants::kGeomFactorsTol);
if(!try_shifting && (lId == -1 || rId == -1 ))
return false;
else if(try_shifting && lId == -1 )
......@@ -233,7 +384,7 @@ namespace Nektar
tmax-=tmin;
if(!checkSupport(point,direction,tmin,tmax,false))
{
tmin=0.;
tmin=0;
tmax=tmax_aux-tmin_aux;
return(checkSupport(point,direction,tmin,tmax,false));
}
......@@ -254,8 +405,6 @@ namespace Nektar
return true;
}
}
}
......
This diff is collapsed.
This diff is collapsed.
......@@ -3,57 +3,42 @@
#include "HandleNekMesh.h"
/// Handles Nektar 1D meshes.
/**
/**
*/
namespace Nektar
{
namespace LSIAC
{
class HandleNekMesh1D;
typedef boost::shared_ptr<HandleNekMesh1D> HandleNekMesh1DSharedPtr;
class HandleNekMesh1D: public HandleNekMesh
{
public:
HandleNekMesh1D( LibUtilities::SessionReaderSharedPtr sessionPtr):HandleNekMesh(sessionPtr)
{
// cout << "into HandleNekMesh1D constructor" << endl;
// m_expansions.push_back(MemoryManager<MultiRegions::ContField1D>::
// AllocateSharedPtr(m_session,m_graph, m_session->GetVariable(0)) );
// cout << "Expansion Type: " << m_expansions[0]->GetExpType() << endl;
}
protected:
virtual bool v_GetBreakPts( const NekDouble xcen_offset,const NekDouble ycen_offset,
const NekDouble zcen_offset, Array<OneD,NekDouble> &direction, const NekDouble t_offset_min,
const NekDouble t_offset_max, vector<NekDouble> &xPos, vector<NekDouble> &yPos,
vector<NekDouble> &zPos, vector<NekDouble> &tPos );
virtual bool v_CanTRangebeApplied( const NekDouble PtsX, const NekDouble PtsY, const NekDouble PtsZ,
const NekDouble scaling, const NekDouble tmin, const NekDouble tmax,
NekDouble &tminUpdate, NekDouble &tmaxUpdate);
virtual bool v_CanTRangebeApplied( const NekDouble PtsX, const NekDouble PtsY, const NekDouble PtsZ,
const Array<OneD,NekDouble> &direction, const NekDouble tmin, const NekDouble tmax,
NekDouble &meshShift);
virtual bool v_LoadData(string Filename, vector<string> &variables);
virtual bool v_LoadMesh(string var);
virtual bool v_EvaluateAt(const Array<OneD,NekDouble> &xPos,const Array<OneD,NekDouble> &yPos,
const Array<OneD,NekDouble> &zPos, const int gID, const int eID,
Array<OneD,NekDouble> &values,int varNum=0);
virtual bool v_GetListOfGIDs( const NekDouble xPos, const NekDouble yPos, const NekDouble zPos,
Array<OneD,NekDouble> &direction, const vector<NekDouble> t_breaks, vector<int> &t_GIDs,
vector<int> & t_EIDs) const;
virtual void v_LoadExpListIntoRTree()
{
assert(false && "Not implemented, should not use it.");
}
};
}
}
{
namespace LSIAC
{
class HandleNekMesh1D;
typedef boost::shared_ptr<HandleNekMesh1D> HandleNekMesh1DSharedPtr;
class HandleNekMesh1D: public HandleNekMesh
{
public:
HandleNekMesh1D( LibUtilities::SessionReaderSharedPtr sessionPtr):HandleNekMesh(sessionPtr)
{
m_useRTree=false;
// cout << "into HandleNekMesh1D constructor" << endl;
// m_expansions.push_back(MemoryManager<MultiRegions::ContField1D>::
// AllocateSharedPtr(m_session,m_graph, m_session->GetVariable(0)) );
// cout << "Expansion Type: " << m_expansions[0]->GetExpType() << endl;
}
protected:
virtual bool v_GetBreakPts( const NekDouble xcen_offset,const NekDouble ycen_offset,
const NekDouble zcen_offset, Array<OneD,NekDouble> &direction, const NekDouble t_offset_min,
const NekDouble t_offset_max, vector<NekDouble> &xPos, vector<NekDouble> &yPos,
vector<NekDouble> &zPos, vector<NekDouble> &tPos );
virtual bool v_LoadData(string Filename, vector<string> &variables);
virtual bool v_LoadMesh(string var);
virtual int v_GetExpansionIndexUsingRTree(Array<OneD,NekDouble> locCoord);
};
}