Commit 083816a2 authored by Michael Turner's avatar Michael Turner

fix for input nek, needs fixing properly, massive code slow down

parent 458f4dbb
......@@ -39,7 +39,7 @@ using namespace std;
namespace Nektar
{
namespace LibUtilities
namespace MeshUtils
{
CADCurve::CADCurve(int i, TopoDS_Shape in) : m_ID(i)
......
......@@ -48,7 +48,7 @@
namespace Nektar
{
namespace LibUtilities
namespace MeshUtils
{
/**
......
......@@ -40,7 +40,7 @@
using namespace std;
namespace Nektar {
namespace LibUtilities {
namespace MeshUtils {
CADSurf::CADSurf(int i, TopoDS_Shape in, vector<EdgeLoop> ein) : m_ID(i), m_edges(ein)
{
......
......@@ -47,7 +47,7 @@
namespace Nektar
{
namespace LibUtilities
namespace MeshUtils
{
class CADCurve;
......
......@@ -46,7 +46,7 @@ using namespace std;
namespace Nektar
{
namespace LibUtilities
namespace MeshUtils
{
string CADSystem::GetName()
......
......@@ -50,7 +50,7 @@
namespace Nektar
{
namespace LibUtilities
namespace MeshUtils
{
/**
......
......@@ -45,7 +45,7 @@
namespace Nektar
{
namespace LibUtilities
namespace MeshUtils
{
/**
......@@ -73,6 +73,15 @@ public:
m_occVert = BRep_Tool::Pnt(TopoDS::Vertex(in));
}
Array<OneD, NekDouble> GetLoc()
{
Array<OneD, NekDouble> out(3);
out[0] = m_occVert.X(); out[1] = m_occVert.Y(); out[2] = m_occVert.Z();
return out;
}
int GetId(){return m_ID;}
private:
/// ID of the vert.
int m_ID;
......
......@@ -6,8 +6,9 @@ IF(NEKTAR_USE_MESH)
# SurfaceMeshing/CurveMesh.cpp
# Octree/Octant.cpp
# Octree/Octree.cpp
# SurfaceMeshing/SurfaceMesh.cpp
# SurfaceMeshing/SurfaceMeshing.cpp
# SurfaceMeshing/SurfaceMeshMesh.cpp
# SurfaceMeshing/SurfaceMeshHOMesh.cpp
# SurfaceMeshing/FaceMesh.cpp
# SurfaceMeshing/OptimiseFunctions.cpp
# SurfaceMeshing/OptimiseUtils.cpp
# ExtLibInterface/TriangleInterface.cpp
......@@ -41,7 +42,7 @@ IF(NEKTAR_USE_MESH)
# Octree/Octant.h
# Octree/Octree.h
# SurfaceMeshing/SurfaceMesh.h
# SurfaceMeshing/SurfaceMeshing.h
# SurfaceMeshing/FaceMesh.h
# ExtLibInterface/TriangleInterface.h
# ExtLibInterface/TetGenInterface.h
# TetMeshing/TetMesh.h
......@@ -67,8 +68,7 @@ ENDIF(NEKTAR_USE_OCC)
ADD_NEKTAR_LIBRARY(MeshUtils lib ${NEKTAR_LIBRARY_TYPE}
${MESHUTILS_SOURCES} ${MESHUTILS_HEADERS})
TARGET_LINK_LIBRARIES(MeshUtils LINK_PUBLIC LibUtilities LocalRegions Collections
MultiRegions SpatialDomains StdRegions SolverUtils)
TARGET_LINK_LIBRARIES(MeshUtils LINK_PUBLIC SolverUtils)
IF(NEKTAR_USE_MESH)
TARGET_LINK_LIBRARIES(MeshUtils LINK_PRIVATE ${TRIANGLE_LIBRARY})
......
......@@ -54,6 +54,10 @@ namespace MeshUtils
LibUtilities::PointsType pCurveType)
: m_n1(pVertex1), m_n2(pVertex2), m_edgeNodes(pEdgeNodes),
m_curveType(pCurveType), m_geom() {}
/// Creates a new linear edge.
Edge(NodeSharedPtr pVertex1, NodeSharedPtr pVertex2)
: m_n1(pVertex1), m_n2(pVertex2), m_edgeNodes(),
m_curveType(), m_geom() {}
/// Copies an existing edge.
Edge(const Edge& pSrc)
: m_n1(pSrc.m_n1), m_n2(pSrc.m_n2), m_edgeNodes(pSrc.m_edgeNodes),
......
......@@ -102,6 +102,7 @@ namespace MeshUtils
LibUtilities::PointsType m_faceCurveType;
};
bool operator==(ElmtConfig const &c1, ElmtConfig const &c2);
/**
* @brief Base class for element definitions.
......
......@@ -87,6 +87,10 @@ namespace MeshUtils
unsigned int m_spaceDim;
/// a order tag to aid output, a bit of a hack
unsigned int m_nummode;
///list of nodes used in meshing routines before its a full mesh
std::vector<NodeSharedPtr> m_meshnode;
///list of edges used in meshing routines before its a full mesh
std::vector<EdgeSharedPtr> m_meshedge;
/// List of mesh nodes.
std::vector<NodeSharedPtr> m_node;
/// Set of element vertices.
......
......@@ -63,6 +63,11 @@ namespace Nektar
return Type::Instance();
}
bool operator==(ElmtConfig const &c1, ElmtConfig const &c2)
{
return (c1.m_e == c2.m_e && c1.m_order == c2.m_order);
}
Element::Element(ElmtConfig pConf,
unsigned int pNumNodes,
unsigned int pGotNodes)
......
......@@ -150,6 +150,36 @@ namespace MeshUtils
m_z*pSrc.m_x-m_x*pSrc.m_z, m_x*pSrc.m_y-m_y*pSrc.m_x);
}
void SetCADCurve(int i, NekDouble t)
{
CADCurve[i] = t;
}
void SetCADSurf(int i, Array<OneD, NekDouble> uv)
{
CADSurf[i] = uv;
}
NekDouble GetCADCurve(int i)
{
std::map<int, NekDouble>::iterator search =
CADCurve.find(i);
ASSERTL0(search->first == i, "node not on this curve");
return search->second;
}
Array<OneD, NekDouble> GetCADSurf(int i)
{
//I dont know why I ahev to do this to get it to work
//this really needs bound checking
std::map<int, Array<OneD, NekDouble> >::iterator search =
CADSurf.find(i);
ASSERTL0(search->first == i,"surface not found");
return search->second;
}
/// Generate a %SpatialDomains::PointGeom for this node.
SpatialDomains::PointGeomSharedPtr GetGeom(int coordDim)
{
......@@ -168,8 +198,10 @@ namespace MeshUtils
NekDouble m_y;
/// Z-coordinate.
NekDouble m_z;
/// Mesh Gen ID
int m_mid;
///list of cadcurves the node lies on
std::map<int, NekDouble> CADCurve;
///list of cadsurfs the node lies on
std::map<int, Array<OneD, NekDouble> > CADSurf;
private:
SpatialDomains::PointGeomSharedPtr m_geom;
......
......@@ -33,105 +33,107 @@
//
////////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_MESHUTILS_OCTREE_CURAVTUREPOINT_H
#define NEKTAR_MESHUTILS_OCTREE_CURAVTUREPOINT_H
#ifndef MESHUTILS_OCTREE_CURAVTUREPOINT
#define MESHUTILS_OCTREE_CURAVTUREPOINT
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
namespace Nektar {
namespace MeshUtils {
namespace Nektar
{
namespace MeshUtils
{
/**
* @brief class for a curvature samlping Point
*/
class CurvaturePoint {
public:
friend class MemoryManager<CurvaturePoint>;
/**
* @brief constructor for a valid point (has radius of curvature)
*/
CurvaturePoint(Array<OneD, NekDouble> l, NekDouble R,
Array<OneD, NekDouble> n) : m_loc(l), m_norm(n),
m_radius(R)
{
m_valid = true;
}
/**
* @brief constructor for a invalid point
*/
CurvaturePoint(Array<OneD, NekDouble> l,
Array<OneD, NekDouble> n) : m_loc(l), m_norm(n)
class CurvaturePoint
{
public:
friend class MemoryManager<CurvaturePoint>;
/**
* @brief constructor for a valid point (has radius of curvature)
*/
CurvaturePoint(Array<OneD, NekDouble> l, NekDouble R,
Array<OneD, NekDouble> n) : m_loc(l), m_norm(n),
m_radius(R)
{
m_valid = true;
}
/**
* @brief constructor for a invalid point
*/
CurvaturePoint(Array<OneD, NekDouble> l,
Array<OneD, NekDouble> n) : m_loc(l), m_norm(n)
{
m_delta = -1;
m_valid = false;
}
/**
* @brief calculate mesh spacing delta from radius of curvature
*/
void Process(NekDouble min, NekDouble max, NekDouble eps)
{
if(m_valid)
{
m_delta = -1;
m_valid = false;
}
/**
* @brief calculate mesh spacing delta from radius of curvature
*/
void Process(NekDouble min, NekDouble max, NekDouble eps)
{
if(m_valid)
{
m_delta = 2.0*m_radius*sqrt(eps*(2.0-eps));
if(m_delta>max)
{
m_delta = max;
}
if(m_delta<min)
{
m_delta = min;
}
}
}
/**
* @brief return bool on whether point is valid
*/
bool IsValid(){return m_valid;}
m_delta = 2.0*m_radius*sqrt(eps*(2.0-eps));
/**
* @brief get mesh spacing paramter
*/
NekDouble GetDelta()
{
if(m_valid)
if(m_delta>max)
{
return m_delta;
m_delta = max;
}
else
if(m_delta<min)
{
return -1;
m_delta = min;
}
}
/**
* @brief get location of point
*/
Array<OneD, NekDouble> GetLoc(){return m_loc;}
/**
* @brief get normal vector
*/
Array<OneD, NekDouble> GetNormal(){return m_norm;}
private:
/// x,y,z location
Array<OneD, NekDouble> m_loc;
///normal vector of surface at point
Array<OneD, NekDouble> m_norm;
/// radius of curvature at point
NekDouble m_radius;
/// mesh spacing parameter at point
NekDouble m_delta;
/// valid point or not
bool m_valid;
}
/**
* @brief return bool on whether point is valid
*/
bool IsValid(){return m_valid;}
/**
* @brief get mesh spacing paramter
*/
NekDouble GetDelta()
{
if(m_valid)
{
return m_delta;
}
else
{
return -1;
}
}
/**
* @brief get location of point
*/
Array<OneD, NekDouble> GetLoc(){return m_loc;}
/**
* @brief get normal vector
*/
Array<OneD, NekDouble> GetNormal(){return m_norm;}
private:
/// x,y,z location
Array<OneD, NekDouble> m_loc;
///normal vector of surface at point
Array<OneD, NekDouble> m_norm;
/// radius of curvature at point
NekDouble m_radius;
/// mesh spacing parameter at point
NekDouble m_delta;
/// valid point or not
bool m_valid;
};
typedef boost::shared_ptr<CurvaturePoint> CurvaturePointSharedPtr;
......
......@@ -169,13 +169,8 @@ void Octree::Modify(Array<OneD, NekDouble> loc, NekDouble delta)
OctantList[n]->SetDelta(delta);
}
void Octree::Build(const NekDouble min, const NekDouble max,
const NekDouble eps)
void Octree::Build()
{
m_minDelta = min;
m_maxDelta = max;
m_eps = eps;
BoundingBox = m_cad->GetBoundingBox();
if(m_verbose)
......@@ -746,7 +741,7 @@ void Octree::CompileCuravturePointList()
for(int i = 1; i <= m_cad->GetNumSurf(); i++)
{
LibUtilities::CADSurfSharedPtr surf = m_cad->GetSurf(i);
CADSurfSharedPtr surf = m_cad->GetSurf(i);
Array<OneD, NekDouble> bounds = surf->GetBounds();
//to figure out the amount of curvature sampling to be conducted on
......@@ -818,15 +813,15 @@ void Octree::CompileCuravturePointList()
//these are the acutal number of sample points in each parametric
//direction
int nu = ceil(DeltaU/m_minDelta)*40*1.2;
int nv = ceil(DeltaV/m_minDelta)*40*1.2;
int nu = ceil(DeltaU/m_minDelta)*40;
int nv = ceil(DeltaV/m_minDelta)*40;
for(int j = 0; j < nu; j++)
{
for(int k = 0; k < nv; k++)
{
Array<OneD, NekDouble> uv(2);
uv[0] = (bounds[1]- bounds[0])/(nu-1)*j + bounds[0];
uv[0] = (bounds[1]-bounds[0])/(nu-1)*j + bounds[0];
uv[1] = (bounds[3]-bounds[2])/(nv-1)*k + bounds[2];
//this prevents round off error at the end of the surface
......@@ -873,17 +868,16 @@ void Octree::CompileCuravturePointList()
if(kv[0] != 0 || kv[1] != 0)
{
CurvaturePointSharedPtr newCPoint =
MemoryManager<CurvaturePoint>::AllocateSharedPtr
(surf->P(uv),
1.0/(kv[0] > kv[1] ? kv[0] : kv[1]),
N);
MemoryManager<CurvaturePoint>::AllocateSharedPtr
(surf->P(uv), 1.0/(kv[0] > kv[1] ? kv[0] : kv[1]), N);
m_cpList.push_back(newCPoint);
}else
{
CurvaturePointSharedPtr newCPoint =
MemoryManager<CurvaturePoint>::AllocateSharedPtr
(surf->P(uv), N);
MemoryManager<CurvaturePoint>::AllocateSharedPtr
(surf->P(uv), N);
m_cpList.push_back(newCPoint);
}
}
......
......@@ -34,8 +34,8 @@
////////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_MESHUTILS_OCTREE_OCTREE_H
#define NEKTAR_MESHUTILS_OCTREE_OCTREE_H
#ifndef NEKTAR_MESHUTILS_OCTREE_OCTREE
#define NEKTAR_MESHUTILS_OCTREE_OCTREE
#include <boost/shared_ptr.hpp>
......@@ -46,8 +46,10 @@
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
namespace Nektar {
namespace MeshUtils {
namespace Nektar
{
namespace MeshUtils
{
/**
* @brief class for octree
......@@ -57,127 +59,125 @@ namespace MeshUtils {
*/
class Octree
{
public:
friend class MemoryManager<Octree>;
/**
* @brief Defualt constructor
*
* @param cad CAD object
* @param ver bool verbose
*/
Octree(const LibUtilities::CADSystemSharedPtr &cad,
const bool ver) : m_cad(cad), m_verbose(ver)
{
};
/**
* @brief executes octree building routines
*
* @param min minimum delta to be found in the octree
* @param max maximum delta to be found in the Octree
* @param eps curvature sensivity parameter
*/
void Build(const NekDouble min, const NekDouble max,
const NekDouble eps);
/**
* @brief once constructed queryies the octree based on x,y,z location
* to get a mesh spacing
*
* @param loc array of x,y,z
* @return mesh spacing parameter
* @todo improve this algorithm for both robustness and completness,
* functions just fine regardless
*/
NekDouble Query(Array<OneD, NekDouble> loc);
void Modify(Array<OneD, NekDouble> loc, NekDouble delta);
/**
* @brief returns the miminum spacing in the octree (for meshing purposes)
*
* @return miminum delta in octree
*/
NekDouble GetMinDelta(){return m_minDelta;}
/**
* @brief Smooths specification over all octants to a
* gradation criteria
*/
void SmoothAllOctants();
private:
/**
* @brief gets an optimum number of curvature sampling points and
* calculates the curavture at these points
*/
void CompileCuravturePointList();
/**
* @brief Recursive alorithm which divides and creates new octants based
* on the geometry
*/
void InitialSubDivide(int parent);
/**
* @brief Recursive alorithm which subdivides octants so that the
* neighbours differ by no more that one level
*/
void SubDivideByLevel();
/**
* @brief Subdivision step for smoothoctants()
*/
void SubDivideLevel(int parent);
/**
* @brief Smooths specification over the surface octants to a
* gradation criteria
*/
void SmoothSurfaceOctants();
/**
* @brief takes the mesh specification from surface octants and
* progates that through the domain so all octants have a specification
* using gradiation crieteria
*/
void PropagateDomain();
/**
* @brief estimates the number of elements to be creted in the mesh
*/
int CountElemt();
/**
* @brief Calculates the difference in delta divided by the difference
* in location between two octants i and j
*/
NekDouble ddx(int i, int j);
/// minimum delta in the octree
NekDouble m_minDelta;
/// maximum delta in the octree
NekDouble m_maxDelta;
/// curavture sensivity paramter
NekDouble m_eps;
/// cad object
LibUtilities::CADSystemSharedPtr m_cad;
/// verbose output?
bool m_verbose;
/// max and min dimensions of the domain 6 varibles