Commit 17d9ef66 authored by Michael Turner's avatar Michael Turner
Browse files

code clean up and commenting, few fixes

parent ac784984
......@@ -81,7 +81,9 @@ void BLMesh::Mesh()
{
if((e[j]->m_n1->GetNumCadCurve() > 0 && e[j]->m_n2->GetNumCadCurve() > 0) ||
(!(e[j]->m_n1->GetNumCadCurve() > 0) && !(e[j]->m_n2->GetNumCadCurve() > 0))) //if both or none are on curve skip
continue;
{
continue;
}
if(e[j]->m_n1->GetNumCadCurve() > 0)
{
......@@ -166,7 +168,7 @@ void BLMesh::Mesh()
ElmtConfig conf(LibUtilities::ePrism,1,false,false);
vector<int> tags;
tags.push_back(11); //need to fix again
tags.push_back(1); //all prisms are comp 1
ElementSharedPtr E = GetElementFactory().
CreateInstance(LibUtilities::ePrism, conf, pn, tags);
......@@ -187,7 +189,7 @@ void BLMesh::Mesh()
if(!(F == f[j])) //only two triangle faces so if its not this one, this is the psudo surfaces
{
surftopriface[ptri[i]->GetId()] = f[j];
m_surftopriface[ptri[i]->GetId()] = f[j];
}
}
}
......
......@@ -57,11 +57,9 @@ public:
/**
*@brief default constructor
*/
BLMesh(MeshSharedPtr m, const std::vector<unsigned int> &bl,
const std::vector<unsigned int> &sym,
const NekDouble &b,
std::map<int, FaceSharedPtr> &stp)
: m_mesh(m), m_blsurfs(bl), m_symsurfs(sym), m_bl(b), surftopriface(stp)
BLMesh(MeshSharedPtr m, std::vector<unsigned int> bls,
std::vector<unsigned int> syms, NekDouble b) :
m_mesh(m), m_blsurfs(bls), m_symsurfs(syms), m_bl(b)
{
};
......@@ -71,15 +69,26 @@ public:
*/
void Mesh();
/**
* @brief Get the map of surface element id to pseudo surface prism face
*/
std::map<int, FaceSharedPtr> GetSurfToPri()
{
return m_surftopriface;
}
private:
/// print stuff to screen?
/// mesh object containing surface mesh
MeshSharedPtr m_mesh;
/// list of surfaces onto which boundary layers are placed
std::vector<unsigned int> m_blsurfs;
/// list of symetry surfaces
std::vector<unsigned int> m_symsurfs;
/// thickness of the boundary layer
NekDouble m_bl;
std::map<int, FaceSharedPtr>& surftopriface;
/// map from surface element id to opposite face of prism
std::map<int, FaceSharedPtr> m_surftopriface;
};
typedef boost::shared_ptr<BLMesh> BLMeshSharedPtr;
......
......@@ -33,10 +33,6 @@
//
////////////////////////////////////////////////////////////////////////////////
#include <string>
#include <sstream>
#include <limits>
#include <boost/algorithm/string.hpp>
#include <boost/filesystem.hpp>
......@@ -59,49 +55,6 @@ void CADSystem::Report()
cout << endl << "CAD report:" << endl;
cout << "\tCAD has: " << m_curves.size() << " curves." << endl;
cout << "\tCAD has: " << m_surfs.size() << " surfaces." << endl;
if(m_surfs.size() == 1 && m_surfs[1]->IsPlane())
{
cout << "\tCAD is 2D" << endl;
m_2d = true;
}
else
{
m_2d = false;
}
}
void CADSystem::SmallFeatureAnalysis(NekDouble min)
{
vector<pair<int,NekDouble> >lens;
map<int, CADCurveSharedPtr>::iterator it;
for(it = m_curves.begin(); it != m_curves.end(); it++)
{
ASSERTL0(it->second->GetTotLength() > 1e-5,
"curve too small for meshing");
if(it->second->GetTotLength() < 5.0 * min)
{
lens.push_back(pair<int,NekDouble>(
it->second->GetID(),it->second->GetTotLength()));
}
}
if(lens.size()>0)
{
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;
}
if(lens.size() == 1)
cout << "\tThis curve's length is less than or close to minDelta" << endl;
else
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;
}
}
Array<OneD, NekDouble> CADSystem::GetBoundingBox()
......
......@@ -77,7 +77,7 @@ public:
std::string GetName();
/**
* @brief Initialises CAD and makes surface and curve maps.
* @brief Initialises CAD and makes surface, curve and vertex maps.
*
* @return true if completed successfully
*/
......@@ -99,7 +99,7 @@ public:
Array<OneD, NekDouble> GetBoundingBox();
/**
* @brief Return number of surfaces.
* @brief Get the number of surfaces.
*/
int GetNumSurf()
{
......@@ -107,7 +107,7 @@ public:
}
/**
* @brief Return number of curves.
* @brief Get the number of curves.
*/
int GetNumCurve()
{
......@@ -115,9 +115,9 @@ public:
}
/**
* @brief Gets curve type from map.
* @brief Gets a curve from the map.
*/
const CADCurveSharedPtr GetCurve(int i)
CADCurveSharedPtr GetCurve(int i)
{
std::map<int,CADCurveSharedPtr>::iterator
search = m_curves.find(i);
......@@ -127,7 +127,7 @@ public:
}
/**
* @brief Gets suface from map.
* @brief Gets a suface from the map.
*/
CADSurfSharedPtr GetSurf(int i)
{
......@@ -139,7 +139,7 @@ public:
}
/**
* @brief Gets map of all verts
* @brief Gets map of all vertices
*/
std::map<int, CADVertSharedPtr> GetVerts()
{
......@@ -159,18 +159,13 @@ public:
*/
bool InsideShape(Array<OneD, NekDouble> loc);
/**
* @brief tests for features which are smaller that the specified mesh limit
*/
void SmallFeatureAnalysis(NekDouble min);
private:
/// Private function to add curve to CADSystem::m_verts.
/// function to add curve to CADSystem::m_verts.
void AddVert(int i, TopoDS_Shape in);
/// Private function to add curve to CADSystem::m_curves.
/// function to add curve to CADSystem::m_curves.
void AddCurve(int i, TopoDS_Shape in, int fv, int lv);
/// Private function to add surface to CADSystem::m_surfs.
/// function to add surface to CADSystem::m_surfs.
void AddSurf(int i, TopoDS_Shape in, std::vector<EdgeLoop> ein);
/// Name of cad file to be opened, including file extension.
std::string m_name;
......@@ -178,12 +173,10 @@ private:
std::map<int, CADCurveSharedPtr> m_curves;
/// map of surfaces
std::map<int, CADSurfSharedPtr> m_surfs;
/// list of edge end vertices
/// map of vertices
std::map<int, CADVertSharedPtr> m_verts;
/// occ master object
TopoDS_Shape shape;
bool m_2d;
};
typedef boost::shared_ptr<CADSystem> CADSystemSharedPtr;
......
......@@ -85,7 +85,6 @@ void TetGenInterface::InitialMesh(map<int, NodeSharedPtr> tgidton,
}
tetrahedralize("pYzqQ", &surface, &output);
cout << output.numberofpoints << " " << output.numberoftetrahedra << endl;
}
......@@ -115,7 +114,6 @@ void TetGenInterface::RefineMesh(std::map<int, NekDouble> delta)
}
tetrahedralize("pYrmzqQO2/7o/120", &input, &output);
cout << output.numberofpoints << " " << output.numberoftetrahedra << endl;
}
......
......@@ -87,6 +87,8 @@ namespace NekMeshUtils
unsigned int m_spaceDim;
/// a order tag to aid output, a bit of a hack
unsigned int m_nummode;
///
unsigned int m_numcomp;
/// List of mesh nodes.
std::vector<NodeSharedPtr> m_node;
/// Set of element vertices.
......
......@@ -33,8 +33,6 @@
//
////////////////////////////////////////////////////////////////////////////////
#include <limits>
#include <NekMeshUtils/Octree/Octant.h>
using namespace std;
......@@ -117,7 +115,7 @@ Octant::Octant(int i, OctantSharedPtr p, Array<OneD, OctantFace> dir) : m_id(i),
//look over the curvature point list provided by the parent,
//firstly look to see if it is in the new octant and if so
//add it to the conserdation of the delta specification
for(int i = 0; i<CurvaturePointList.size(); i++)
for(int i = 0; i < CurvaturePointList.size(); i++)
{
Array<OneD, NekDouble> cploc = CurvaturePointList[i]->GetLoc();
if(!(cploc[0] > m_loc[0] + m_hd ||
......@@ -127,19 +125,19 @@ Octant::Octant(int i, OctantSharedPtr p, Array<OneD, OctantFace> dir) : m_id(i),
cploc[2] > m_loc[2] + m_hd ||
cploc[2] < m_loc[2] - m_hd ))
{
m_localCPIDList.push_back(CurvaturePointList[i]);
m_localCPList.push_back(CurvaturePointList[i]);
if(CurvaturePointList[i]->IsValid())
{
if(CurvaturePointList[i]->GetDelta()>maxDif)
if(CurvaturePointList[i]->GetDelta() > maxDif)
{
maxDif = CurvaturePointList[i]->GetDelta();
}
if(CurvaturePointList[i]->GetDelta()<minDif)
if(CurvaturePointList[i]->GetDelta() < minDif)
{
minDif = CurvaturePointList[i]->GetDelta();
}
m_numValidPoints++;
}
}
......@@ -186,11 +184,11 @@ Octant::Octant(int i, NekDouble x, NekDouble y, NekDouble z, NekDouble dx,
m_loc[1] = y;
m_loc[2] = z;
m_localCPIDList = cplist;
m_localCPList = cplist;
for(int i = 0; i < m_localCPIDList.size(); i++)
for(int i = 0; i < m_localCPList.size(); i++)
{
if(m_localCPIDList[i]->IsValid())
if(m_localCPList[i]->IsValid())
{
m_numValidPoints++;
}
......
......@@ -33,7 +33,6 @@
//
////////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_MESHUTILS_OCTREE_OCTANT_H
#define NEKTAR_MESHUTILS_OCTREE_OCTANT_H
......@@ -42,8 +41,6 @@
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <boost/unordered_set.hpp>
namespace Nektar
{
namespace NekMeshUtils
......@@ -69,7 +66,6 @@ enum OctantLocation
class Octant; //have to forward declare the class for the sharedptr
typedef boost::shared_ptr<Octant> OctantSharedPtr;
typedef std::set<OctantSharedPtr> OctantSet;
/**
* @brief this class contains the infomration and methods for individal octants
......@@ -85,12 +81,21 @@ class Octant
*/
Octant(int i, OctantSharedPtr p, Array<OneD, OctantFace> dir);
//master constructor
/**
* @brief constructor for master octant
*/
Octant(int i, NekDouble x, NekDouble y, NekDouble z, NekDouble dx,
const std::vector<CurvaturePointSharedPtr> &cplist);
/**
* @brief Subdivide the octant
*/
void Subdivide(OctantSharedPtr p, int &numoct);
/**
* @brief Recursive method which gets a list of the leaf octants
* Moves down levels if octant is not a leaf
*/
void CompileLeaves(std::vector<OctantSharedPtr> &Octants)
{
for(int i = 0; i < 8; i++)
......@@ -106,42 +111,64 @@ class Octant
}
}
/// Get the Id of the octant
int GetId()
{
return m_id;
}
/**
* @brief Get the location of the center of the octant
*/
Array<OneD, NekDouble> GetLoc()
{
return m_loc;
}
/**
* @brief Get the octants half dimension
*/
NekDouble DX()
{
return m_hd;
}
/**
* @brief Get the list of curvature points that are within this octant
*/
std::vector<CurvaturePointSharedPtr> GetCPList()
{
return m_localCPIDList;
return m_localCPList;
}
/**
* @brief Get the number of points in the octants cp list
*/
int NumCurvePoint()
{
return m_localCPIDList.size();
return m_localCPList.size();
}
/**
* @brief Get the number of valid cp points in the octants list
*/
int NumValidCurvePoint()
{
return m_numValidPoints;
}
/**
* @brief Set the value for delta for this octant
*/
void SetDelta(NekDouble d)
{
m_delta.first = true;
m_delta.second = d;
}
/**
* @brief Get value of delta
*/
NekDouble GetDelta()
{
ASSERTL0(m_delta.first, "Tried to acsess delta of octant"
......@@ -149,16 +176,25 @@ class Octant
return m_delta.second;
}
/**
* @brief Set the children of this octant
*/
void SetChildren(Array<OneD, OctantSharedPtr> c)
{
m_children = c;
}
/**
* @brief Get whether the octant is a leaf or not
*/
bool IsLeaf()
{
return m_leaf;
}
/**
* @brief Get the far dimension in a given direction f
*/
NekDouble FX(OctantFace f)
{
switch (f)
......@@ -184,6 +220,9 @@ class Octant
}
}
/**
* @brief Remove a neigbour from this octants list
*/
void RemoveNeigbour(int id, OctantFace f);
void SetNeigbour(OctantSharedPtr o, OctantFace f)
......@@ -191,16 +230,25 @@ class Octant
m_neigbours[f].push_back(o);
}
/**
* @brief Get the map of neigbours
*/
std::map<OctantFace, std::vector<OctantSharedPtr> > GetNeigbours()
{
return m_neigbours;
}
/**
* @brief Get whether this octants needs to divide based on the curvature sampling
*/
bool NeedDivide()
{
return m_needToDivide;
}
/**
* @brief Get the distance from this octant to another
*/
NekDouble Distance(OctantSharedPtr o)
{
Array<OneD, NekDouble> loc = o->GetLoc();
......@@ -209,27 +257,42 @@ class Octant
(loc[2] - m_loc[2])*(loc[2] - m_loc[2]));
}
/**
* @brief Get whether a value of delta has been set or not
*/
bool IsDeltaKnown()
{
return m_delta.first;
}
/**
* @brief set the location of the octant with respect to the CAD
*/
void SetLocation(OctantLocation l)
{
m_location = l;
}
/**
* @brief Get the first cp point in the list
*/
CurvaturePointSharedPtr GetCPPoint()
{
ASSERTL0(m_localCPIDList.size() > 0, "tried to get cp point where there is none");
return m_localCPIDList[0];
ASSERTL0(m_localCPList.size() > 0, "tried to get cp point where there is none");
return m_localCPList[0];
}
/**
* @brief Get the location of the octant with respect to the geometry
*/
OctantLocation GetLocation()
{
return m_location;
}
/**
* @brief Get a specific child of this octant
*/
OctantSharedPtr GetChild(int q)
{
return m_children[q];
......@@ -250,7 +313,8 @@ class Octant
/// half dimension of the octant
NekDouble m_hd;
/// curvature sampling point list
std::vector<CurvaturePointSharedPtr> m_localCPIDList;
std::vector<CurvaturePointSharedPtr> m_localCPList;
/// number of valid cp points
int m_numValidPoints;
/// mesh sizing parameter
std::pair<bool, NekDouble> m_delta;
......
......@@ -33,10 +33,8 @@
//
////////////////////////////////////////////////////////////////////////////////
#include <algorithm>
#include <limits>
#include <LibUtilities/BasicUtils/Progressbar.hpp>
//#include <algorithm>
//#include <limits>
#include <NekMeshUtils/Octree/Octree.h>
#include <NekMeshUtils/CADSystem/CADSurf.h>
......@@ -49,7 +47,7 @@ namespace NekMeshUtils
NekDouble Octree::Query(Array<OneD, NekDouble> loc)
{
//starting at master octant 0 move through succsesive octants which contain
//starting at master octant 0 move through succsesive m_octants which contain
//the point loc until a leaf is found
OctantSharedPtr n = m_masteroct;
int quad;
......@@ -123,72 +121,53 @@ NekDouble Octree::Query(Array<OneD, NekDouble> loc)
return n->GetDelta();
}
vector<Array<OneD, Array<OneD, NekDouble> > > Octree::GetOctantVerts()
void Octree::GetOctreeMesh(MeshSharedPtr m)
{
vector<Array<OneD, Array<OneD, NekDouble> > > out;
for(int i = 0; i < Octants.size(); i++)
for(int i = 0; i < m_octants.size(); i++)
{
/*if(Octants[i]->GetLocation() != eOnBoundary)
/*if(m_octants[i]->GetLocation() != eOnBoundary)
{
continue;
}*/
Array<OneD, Array<OneD, NekDouble> > oct(8);
Array<OneD, NekDouble> p1(3);
p1[0] = Octants[i]->FX(eBack);
p1[1] = Octants[i]->FX(eDown);
p1[2] = Octants[i]->FX(eRight);
oct[0] = p1;
Array<OneD, NekDouble> p2(3);
p2[0] = Octants[i]->FX(eForward);
p2[1] = Octants[i]->FX(eDown);
p2[2] = Octants[i]->FX(eRight);
oct[1] = p2;
Array<OneD, NekDouble> p3(3);
p3[0] = Octants[i]->FX(eForward);
p3[1] = Octants[i]->FX(eUp);
p3[2] = Octants[i]->FX(eRight);
oct[2] = p3;
Array<OneD, NekDouble> p4(3);
p4[0] = Octants[i]->FX(eBack);
p4[1] = Octants[i]->FX(eUp);
p4[2] = Octants[i]->FX(eRight);
oct[3] = p4;
Array<OneD, NekDouble> p5(3);
p5[0] = Octants[i]->FX(eBack);
p5[1] = Octants[i]->FX(eDown);
p5[2] = Octants[i]->FX(eLeft);
oct[4] = p5;
Array<OneD, NekDouble> p6(3);
p6[0] = Octants[i]->FX(eForward);
p6[1] = Octants[i]->FX(eDown);
p6[2] = Octants[i]->FX(eLeft);
oct[5] = p6;
Array<OneD, NekDouble> p7(3);
p7[0] = Octants[i]->FX(eForward);