Commit e9e76af7 authored by Michael Turner's avatar Michael Turner
Browse files

made octree unordered set, element count broken, possibly normal alignment is off

parent 966a7933
......@@ -78,6 +78,7 @@ void CADSystem::SmallFeatureAnalysis(NekDouble min)
{
cout << endl;
cout << "Small feature testing:" << endl;
cout << "\tWARNING" << endl;
for(int i = 0; i < lens.size(); i++)
{
cout << "\tCurve: " << lens[i].first << " has length: " << lens[i].second << endl;
......@@ -88,15 +89,6 @@ void CADSystem::SmallFeatureAnalysis(NekDouble min)
cout << "\tThese curve's length are less than or close to minDelta" << endl;
cout << "\tIf these featurs are not planar the mesh generator may struggle ";
cout << "to produce a high-order mesh due to high curvature" << endl;
cout << "\tWould you like the Octree system to adjust accordingly, this may ";
cout << "compramise mesh smoothness (Y/N)" << endl;
string ans;
cin >> ans;
if(boost::iequals(ans,"Y"))
smallfeatreAdjust = true;
else
smallfeatreAdjust = false;
}
}
......
......@@ -43,7 +43,7 @@ namespace Nektar
namespace MeshUtils
{
Octant::Octant(OctantSharedPtr p, Array<OneD, NekDouble> dir) : m_parent(p)
Octant::Octant(int i, OctantSharedPtr p, Array<OneD, NekDouble> dir) : m_id(i), m_parent(p)
{
//initialise variables to defualt states
m_leaf = true;
......@@ -128,14 +128,13 @@ Octant::Octant(OctantSharedPtr p, Array<OneD, NekDouble> dir) : m_parent(p)
}
//constructor for the master octant
Octant::Octant(NekDouble x, NekDouble y, NekDouble z, NekDouble dx,
Octant::Octant(int i, NekDouble x, NekDouble y, NekDouble z, NekDouble dx,
const vector<CurvaturePointSharedPtr> &cplist)
: m_hd(dx)
: m_id(i), m_hd(dx)
{
//initialise variables to defualt states
m_leaf = false;
m_needToDivide = true;
m_deltaSet = true;
m_numValidPoints = 0;
m_delta = -1;
......@@ -158,109 +157,9 @@ Octant::Octant(NekDouble x, NekDouble y, NekDouble z, NekDouble dx,
m_location = 2;
}
void Octant::CreateNeighbourList(OctantSet Octants)
{
//clear old list
DeleteNeighbourList();
OctantSet::iterator it;
//look over all octants and consider if they are neigbours
for(it = Octants.begin(); it != Octants.end(); it++)
{
OctantSharedPtr oct = *it;
if(oct == boost::shared_ptr<Octant>(this))
continue;
if(oct->IsLeaf())
{
//work out the max distance between the two octants if they were
//joined corner to corner
NekDouble rmax = DiagonalDim() + oct->DiagonalDim();
NekDouble ractual = Distance(oct);
//if the actucal distance between them is greater that this
//this octant should not be considered, massive speed up
if(ractual > 1.1*rmax)
{
continue;
}
//check overlapping in 2-D for all three dimensions
if(fabs(FX(-1) - oct->FX(+1)) < 1E-6 ||
fabs(FX(+1) - oct->FX(-1)) < 1E-6 )
{
//check yz rects
bool Cond1 = false;
bool Cond2 = false;
bool Cond3 = false;
bool Cond4 = false;
if(FY(+1) < oct->FY(-1))
Cond1 = true;
if(FY(-1) > oct->FY(+1))
Cond1 = true;
if(FZ(+1) < oct->FZ(-1))
Cond1 = true;
if(FZ(-1) > oct->FZ(+1))
Cond1 = true;
if(!Cond1 && !Cond2 && !Cond3 && !Cond4)
{
m_neighbourList.push_back(oct);
}
}
else if(fabs(FY(-1) - oct->FY(+1)) <1E-6 ||
fabs(FY(+1) - oct->FY(-1)) <1E-6)
{
//check xz rects
bool Cond1 = false;
bool Cond2 = false;
bool Cond3 = false;
bool Cond4 = false;
if(FX(+1) < oct->FX(-1))
Cond1 = true;
if(FX(-1) > oct->FX(+1))
Cond1 = true;
if(FZ(+1) < oct->FZ(-1))
Cond1 = true;
if(FZ(-1) > oct->FZ(+1))
Cond1 = true;
if(!Cond1 && !Cond2 && !Cond3 && !Cond4)
{
m_neighbourList.push_back(oct);
}
}
else if(fabs(FZ(-1) - oct->FZ(+1)) <1E-6 ||
fabs(FZ(+1) - oct->FZ(-1)) <1E-6)
{
//check xy rects
bool Cond1 = false;
bool Cond2 = false;
bool Cond3 = false;
bool Cond4 = false;
if(FX(+1) < oct->FX(-1))
Cond1 = true;
if(FX(-1) > oct->FX(+1))
Cond1 = true;
if(FY(+1) < oct->FY(-1))
Cond1 = true;
if(FY(-1) > oct->FY(+1))
Cond1 = true;
if(!Cond1 && !Cond2 && !Cond3 && !Cond4)
{
m_neighbourList.push_back(oct);
}
}
}
}
}
bool operator==(OctantSharedPtr const &p1, OctantSharedPtr const &p2)
{
Array<OneD, NekDouble> loc1 = p1->GetLoc();
Array<OneD, NekDouble> loc2 = p2->GetLoc();
if(loc1[0] == loc2[0] && loc1[1] == loc2[1] && loc1[2] == loc2[2])
if(p1->GetId() == p2->GetId())
return true;
else
return false;
......
......@@ -51,7 +51,7 @@ namespace MeshUtils
class Octant; //have to forward declare the class for the sharedptr
typedef boost::shared_ptr<Octant> OctantSharedPtr;
typedef boost::unordered_set<OctantSharedPtr> OctantSet;
typedef std::set<OctantSharedPtr> OctantSet;
/**
* @brief this class contains the infomration and methods for individal octants
......@@ -65,17 +65,19 @@ class Octant
/**
* @brief Defualt constructor
*/
Octant(OctantSharedPtr p, Array<OneD, NekDouble> dir);
Octant(int i, OctantSharedPtr p, Array<OneD, NekDouble> dir);
//master constructor
Octant(NekDouble x, NekDouble y, NekDouble z, NekDouble dx,
Octant(int i, NekDouble x, NekDouble y, NekDouble z, NekDouble dx,
const std::vector<CurvaturePointSharedPtr> &cplist);
/**
* @brief scans over all octants in the octantlist and finds neighouring
* octants
*/
void CreateNeighbourList(OctantSet OctantList);
void ClearNeigbourList(){m_neighbourList.clear();}
void SetNeigbourList(std::vector<OctantSharedPtr> const &l)
{m_neighbourList = l;}
int GetId(){return m_id;}
/**
* @brief get boolean on whether the octant needs to divide based on
......@@ -93,11 +95,6 @@ class Octant
*/
bool IsLeaf(){return m_leaf;}
/**
* @brief alter leaf status
*/
void SetLeaf(bool l){m_leaf = l;}
/**
* @brief get boolean on whether the octant has curvature sample points
* within its volume
......@@ -191,12 +188,7 @@ class Octant
/**
* @brief set the list of child octants
*/
void SetChildren(Array<OneD, OctantSharedPtr> i){m_children = i;}
/**
* @brief clear neigbour list
*/
void DeleteNeighbourList(){m_neighbourList.clear();}
void SetChildren(Array<OneD, OctantSharedPtr> i){m_children = i; m_leaf=false;}
/**
* @brief get the parent of this octant
......@@ -247,6 +239,8 @@ class Octant
private:
///id
int m_id;
/// leaf identifer
bool m_leaf;
/// parent id
......
......@@ -146,8 +146,8 @@ void Octree::Build()
//make master octant based on the bounding box of the domain
OctantSharedPtr newOctant =
MemoryManager<Octant>::AllocateSharedPtr
((BoundingBox[1]+BoundingBox[0])/2,
MemoryManager<Octant>::AllocateSharedPtr(0,
(BoundingBox[1]+BoundingBox[0])/2,
(BoundingBox[3]+BoundingBox[2])/2,
(BoundingBox[5]+BoundingBox[4])/2, maxdim, m_cpList);
......@@ -156,31 +156,10 @@ void Octree::Build()
m_totNotDividing=0;
if(m_masteroct->GetDivide())
{
m_masteroct->SetLeaf(false);
InitialSubDivide(m_masteroct);
}
InitialSubDivide(m_masteroct);
OctantSet::iterator it;
if(m_verbose)
{
int ct=0;
int maxLevel=0;
for(it = Octants.begin(); it != Octants.end(); it++)
{
OctantSharedPtr oct = *it;
if(oct->IsLeaf())
ct++;
if(oct->GetLevel() > maxLevel)
maxLevel = oct->GetLevel();
}
cout << "No. octant leaves: " << ct <<
", Max octree level: " << maxLevel << endl;
}
int ct = 0;
for(it = Octants.begin(); it != Octants.end(); it++, ct++)
{
......@@ -192,25 +171,15 @@ void Octree::Build()
}
if(oct->IsLeaf())
{
oct->CreateNeighbourList(Octants);
AssignNeigbours(oct);
}
}
if(m_verbose)
cout << endl;
SubDivideByLevel();
//memory is good down to hear
ct = 0;
if(m_verbose)
{
for(it = Octants.begin(); it != Octants.end(); it++)
{
if((*it)->IsLeaf()){ct++;}
}
cout << "\tNew Stats: ";
cout << "No. octant leaves: " << ct << endl;
cout << "\tSmoothing across the geometry surface" << endl;
}
SubDivideByLevel();
SmoothSurfaceOctants();
......@@ -236,7 +205,7 @@ void Octree::Build()
if(m_verbose)
{
int elem = CountElemt();
cout << endl<< "\tPredicted mesh: " << elem << " elements" << endl;
printf("\tPredicted mesh: %.2d elements\n", elem);
}
exit(-1);
......@@ -248,8 +217,6 @@ void Octree::InitialSubDivide(OctantSharedPtr parent)
//if that child also needs sub dividing, call this function recursively
Array<OneD, OctantSharedPtr> children(8);
Array<OneD, NekDouble> parentloc = parent->GetLoc();
for(int i = 0; i < 8; i++)
{
//set up x,y,z ordering of the 8 octants
......@@ -296,16 +263,16 @@ void Octree::InitialSubDivide(OctantSharedPtr parent)
}
OctantSharedPtr newOctant = MemoryManager<Octant>::AllocateSharedPtr
(parent, dir);
(Octants.size(), parent, dir);
Octants.insert(newOctant);
pair<OctantSet::iterator, bool> test = Octants.insert(newOctant);
ASSERTL0(test.second,"insertion should not fail");
children[i] = newOctant;
if(children[i]->GetDivide())
{
if(children[i]->DX() / 2.0 > m_minDelta)
{
children[i]->SetLeaf(false);
InitialSubDivide(children[i]);
}
}
......@@ -319,19 +286,23 @@ void Octree::SubDivideByLevel()
{
//until all subdivision ceases, evaluate each octant and subdivide if
//the neigbour levels are not smooth.
int ct=0;
bool brk;
int j = 0;
int imax=0;
OctantSet::iterator it;
do
{
ct=0;
brk = false;
j = 0;
OctantSharedPtr oct;
Array<OneD, OctantSharedPtr> newoctants;
for(it = Octants.begin(); it != Octants.end(); it++, j++)
{
OctantSharedPtr oct = *it;
if(j > imax)
imax = j;
oct = *it;
if(oct->IsLeaf())
{
vector<OctantSharedPtr> nList = oct->GetNeighbourList();
......@@ -340,36 +311,50 @@ void Octree::SubDivideByLevel()
{
if(nList[k]->GetLevel() - oct->GetLevel() > 1)
{
ct+=1;
if(j > imax)
imax = j;
SubDivideLevel(oct);
brk = true;
newoctants = SubDivideLevel(oct);
break;
}
}
if(ct > 0)
if(brk)
break;
}
nList.clear();
if(m_verbose)
{
LibUtilities::PrintProgressbar(imax, Octants.size(),
"\tSubdivide by level\t");
}
}
if(m_verbose)
if(brk)
{
LibUtilities::PrintProgressbar(imax, Octants.size(),
"\tSubdivide by level\t");
oct->SetChildren(newoctants);
for(int i = 0; i < 8; i++)
{
pair<OctantSet::iterator, bool> test = Octants.insert(newoctants[i]);
ASSERTL0(test.second,"insertion should not fail");
}
for(int i = 0; i < 8; i++)
{
AssignNeigbours(newoctants[i]);
}
vector<OctantSharedPtr> nList = oct->GetNeighbourList();
for(int i = 0; i < nList.size(); i++)
{
AssignNeigbours(nList[i]);
}
}
}while(ct>0);
}while(brk);
if(m_verbose)
cout <<endl;
}
void Octree::SubDivideLevel(OctantSharedPtr parent)
Array<OneD, OctantSharedPtr> Octree::SubDivideLevel(OctantSharedPtr parent)
{
//create 8 child octants in turn for octant parent
//after creation, re-evaluate all nessercary neighbour lists
parent->SetLeaf(false);
Array<OneD, OctantSharedPtr> children(8);
......@@ -418,25 +403,25 @@ void Octree::SubDivideLevel(OctantSharedPtr parent)
}
OctantSharedPtr newOctant = MemoryManager<Octant>::AllocateSharedPtr
(parent, dir);
Octants.insert(newOctant);
(Octants.size()+i, parent, dir);
children[i] = newOctant;
}
parent->SetChildren(children);
return children;
/*parent->SetChildren(children);
for(int i = 0; i < 8; i++)
{
children[i]->CreateNeighbourList(Octants);
AssignNeigbours(children[i]);
}
//need to revaluate the neighbour list of all the neighbours of the parent
vector<OctantSharedPtr> nList = parent->GetNeighbourList();
for(int i = 0; i < nList.size(); i++)
{
nList[i]->CreateNeighbourList(Octants);
}
AssignNeigbours(nList[i]);
}*/
}
void Octree::SmoothSurfaceOctants()
......@@ -625,26 +610,19 @@ void Octree::PropagateDomain()
if(dot < 0.0)
{
oct->SetLocation(3); //outside the domain
oct->SetLocation(1); //outside the domain
}
else
{
oct->SetLocation(1); //inside
oct->SetLocation(3); //inside
}
}
else
{
int newloc = known[0]->GetLocation();
for(int j = 1; j < known.size(); j++)
{
if(known[j]->GetLocation() != newloc)
{
ASSERTL0(false,"conflicting locations which"
"should be the same");
}
}
oct->SetLocation(newloc);
//doing this may cause issues, needs testing and fixing
}
}
known.clear();
......@@ -720,47 +698,48 @@ int Octree::CountElemt()
//by considering the volume of a tet evaluate the number of elements in the
//mesh
Array<OneD, NekDouble> box = m_cad->GetBoundingBox();
NekDouble total = 0.0;
NekDouble altot = 0.0;
OctantSet::iterator it;
int ctcor = 0;
int ctinc = 0;
int c2 = 0;
/*for(int i = 0 ; i < OctantList.size(); i++)
for(it = Octants.begin(); it != Octants.end(); it++)
{
if(OctantList[i]->GetLeaf() && OctantList[i]->HasPoints())
{
NekDouble xmin, xmax, ymin, ymax, zmin, zmax;
xmin = max(box[0], OctantList[i]->FX(-1.0));
xmax = min(box[1], OctantList[i]->FX(+1.0));
ymin = max(box[2], OctantList[i]->FY(-1.0));
ymax = min(box[3], OctantList[i]->FY(+1.0));
zmin = max(box[4], OctantList[i]->FZ(-1.0));
zmax = min(box[5], OctantList[i]->FZ(+1.0));
ASSERTL0(xmin < xmax, "error");
ASSERTL0(ymin < ymax, "error");
ASSERTL0(zmin < zmax, "error");
NekDouble voloverlap = (xmax - xmin)*(ymax - ymin)*(zmax - zmin);
NekDouble volumeTet = OctantList[i]->GetDelta()*
OctantList[i]->GetDelta()*
OctantList[i]->GetDelta()/6.0/sqrt(2.0);
total += voloverlap/volumeTet;
}
else if(OctantList[i]->GetLeaf())
OctantSharedPtr oct = *it;
if(oct->IsLeaf())
{
Array<OneD, NekDouble> loc = OctantList[i]->GetLoc();
if(m_cad->InsideShape(loc))
if(oct->GetLocation() == 1)
{
NekDouble volumeTet = OctantList[i]->GetDelta()*
OctantList[i]->GetDelta()*
OctantList[i]->GetDelta()/6.0/sqrt(2.0);
total += 8.0*oct->DX()*oct->DX()*oct->DX() /
(oct->GetDelta()*oct->GetDelta()*oct->GetDelta()/6.0/sqrt(2));
NekDouble volumeOct = OctantList[i]->DX()*OctantList[i]->DX()*
OctantList[i]->DX()*8.0;
Array<OneD, NekDouble> loc = oct->GetLoc();
if(!m_cad->InsideShape(loc))
{
ctinc++;
}
else
{
ctcor++;
}
}
if(oct->GetLocation() == 2)
c2++;
total += volumeOct/volumeTet;
Array<OneD, NekDouble> loc = oct->GetLoc();
if(m_cad->InsideShape(loc) || oct->GetLocation() == 2)
{
altot += 8.0*oct->DX()*oct->DX()*oct->DX() /
(oct->GetDelta()*oct->GetDelta()*oct->GetDelta()/6.0/sqrt(2));
}
}
}*/
}
cout << ctcor << " correct " << ctinc << " incorcect " << c2 << " in 2 loc" << endl;
cout << "alternative total " << int(altot) << endl;
return int(total);
}
......@@ -914,5 +893,112 @@ NekDouble Octree::ddx(OctantSharedPtr i, OctantSharedPtr j)
return fabs(i->GetDelta() - j->GetDelta()) / i->Distance(j);
}
void Octree::AssignNeigbours(OctantSharedPtr const &o)
{
//clear old list
o->ClearNeigbourList();
OctantSet::iterator it;
vector<OctantSharedPtr> nlist;
int eqhit = 0;
//look over all octants and consider if they are neigbours
for(it = Octants.begin(); it != Octants.end(); it++)
{
OctantSharedPtr oct = *it;
if(oct->IsLeaf())
{
if(oct == o)
{
eqhit++;
continue;
}
//work out the max distance between the two octants if they were
//joined corner to corner
NekDouble rmax = o->DiagonalDim() + oct->DiagonalDim();
NekDouble ractual = o->Distance(oct);
//if the actucal distance between them is greater that this
//this octant should not be considered, massive speed up
if(ractual > 1.1*rmax)
{