Commit 8d610a03 authored by Michael Turner's avatar Michael Turner
Browse files

full variable bl thickness

parents b4bf156e b2fce300
......@@ -83,6 +83,11 @@ v4.4.0
(!712)
- 2D to 3D mesh extrusion module (!715)
- Add new two-dimensional mesher from NACA code or step file (!720)
- Fix inverted boundary layer in 2D (!736)
- More sensible element sizing with boudnary layers in 2D (!736)
- Change variable names in mcf file to make more sense (!736)
- Fix issues in varopti module so that in can be compiled without meshgen on
(!736)
**FieldConvert:**
- Move all modules to a new library, FieldUtils, to support post-processing
......
......@@ -100,19 +100,19 @@ class ThreadJob
{
public:
/// Base constructor
ThreadJob();
LIB_UTILITIES_EXPORT ThreadJob();
/// Base destructor.
virtual ~ThreadJob();
LIB_UTILITIES_EXPORT virtual ~ThreadJob();
/**
* This method will be called when the task is loaded
* onto a worker thread and is ready to run. When Run
* has finished this instance will be destructed.
*/
virtual void Run() = 0;
LIB_UTILITIES_EXPORT virtual void Run() = 0;
/// Set number of worker threads.
void SetWorkerNum(unsigned int num);
LIB_UTILITIES_EXPORT void SetWorkerNum(unsigned int num);
protected:
/**
......@@ -168,7 +168,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
{
public:
/// Destructor.
virtual ~ThreadManager();
LIB_UTILITIES_EXPORT virtual ~ThreadManager();
/**
* @brief Pass a list of tasklets to the master queue.
* @param joblist Vector of ThreadJob pointers.
......@@ -181,7 +181,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
*
* @see SchedType
*/
virtual void QueueJobs(std::vector<ThreadJob*>& joblist) = 0;
LIB_UTILITIES_EXPORT virtual void QueueJobs(std::vector<ThreadJob*>& joblist) = 0;
/**
* @brief Pass a single job to the master queue.
* @param job A pointer to a ThreadJob subclass.
......@@ -190,14 +190,14 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* issue then suspend the workers with SetNumWorkers(0) until the jobs
* are queued.
*/
virtual void QueueJob(ThreadJob* job) = 0;
LIB_UTILITIES_EXPORT virtual void QueueJob(ThreadJob* job) = 0;
/**
* @brief Return the number of active workers.
*
* Active workers are threads that are either running jobs
* or are waiting for jobs to be queued.
*/
virtual unsigned int GetNumWorkers() = 0;
LIB_UTILITIES_EXPORT virtual unsigned int GetNumWorkers() = 0;
/**
* @brief Returns the worker number of the executing thread.
*
......@@ -216,7 +216,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
*
* Returns 0 if called by non-thread.
*/
virtual unsigned int GetWorkerNum() = 0;
LIB_UTILITIES_EXPORT virtual unsigned int GetWorkerNum() = 0;
/**
* @brief Sets the number of active workers.
* @param num The number of active workers.
......@@ -227,18 +227,18 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* If num is greater than the maximum allowed number of active workers,
* then the maximum value will be used instead.
*/
virtual void SetNumWorkers(const unsigned int num) = 0;
LIB_UTILITIES_EXPORT virtual void SetNumWorkers(const unsigned int num) = 0;
/**
* @brief Sets the number of active workers to the maximum.
*
* Sets the number of active workers to the maximum available.
*/
virtual void SetNumWorkers() = 0;
LIB_UTILITIES_EXPORT virtual void SetNumWorkers() = 0;
/**
* @brief Gets the maximum available number of threads.
* @return The maximum number of workers.
*/
virtual unsigned int GetMaxNumWorkers() = 0;
LIB_UTILITIES_EXPORT virtual unsigned int GetMaxNumWorkers() = 0;
/**
* @brief Waits until all queued jobs are finished.
*
......@@ -261,7 +261,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* number of worker threads, implementations should increase the number
* of active workers by 1 on entering Wait().
*/
virtual void Wait() = 0;
LIB_UTILITIES_EXPORT virtual void Wait() = 0;
/**
* @brief Controls how many jobs are sent to each worker at a time.
*
......@@ -271,17 +271,17 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* @see SchedType
* @see SetSchedType()
*/
virtual void SetChunkSize(unsigned int chnk) = 0;
LIB_UTILITIES_EXPORT virtual void SetChunkSize(unsigned int chnk) = 0;
/**
* @brief Sets the current scheduling algorithm.
* @see SetChunkSize()
*/
virtual void SetSchedType(SchedType s) = 0;
LIB_UTILITIES_EXPORT virtual void SetSchedType(SchedType s) = 0;
/**
* @brief Indicates whether the code is in a worker thread or not.
* @return True if the caller is in a worker thread.
*/
virtual bool InThread() = 0;
LIB_UTILITIES_EXPORT virtual bool InThread() = 0;
/**
* @brief A calling threads holds until all active threads call this
* method.
......@@ -294,16 +294,16 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* is altered after a thread has called this method. It is only safe to
* call SetNumWorkers() when no threads are holding.
*/
virtual void Hold() = 0;
LIB_UTILITIES_EXPORT virtual void Hold() = 0;
/**
* @brief Returns a description of the type of threading.
*
* E.g. "Threading with Boost"
*/
virtual const std::string& GetType() const = 0;
LIB_UTILITIES_EXPORT virtual const std::string& GetType() const = 0;
/// ThreadManager implementation.
virtual bool IsInitialised();
LIB_UTILITIES_EXPORT virtual bool IsInitialised();
inline int GetThrFromPartition(int pPartition)
{
......@@ -354,17 +354,17 @@ class ThreadMaster
THREADMANAGER_MAX
};
/// Constructor
ThreadMaster();
LIB_UTILITIES_EXPORT ThreadMaster();
/// Destructor
~ThreadMaster();
LIB_UTILITIES_EXPORT ~ThreadMaster();
/// Sets what ThreadManagers will be created in CreateInstance.
void SetThreadingType(const std::string &p_type);
LIB_UTILITIES_EXPORT void SetThreadingType(const std::string &p_type);
/// Gets the ThreadManager associated with string s.
ThreadManagerSharedPtr& GetInstance(const ThreadManagerName t);
LIB_UTILITIES_EXPORT ThreadManagerSharedPtr& GetInstance(const ThreadManagerName t);
/// Creates an instance of a ThreadManager (which one is determined by
/// a previous call to SetThreadingType) and associates it with
/// the string s.
ThreadManagerSharedPtr CreateInstance(const ThreadManagerName t,
LIB_UTILITIES_EXPORT ThreadManagerSharedPtr CreateInstance(const ThreadManagerName t,
unsigned int nThr);
};
......
......@@ -54,7 +54,7 @@ Generator2D::Generator2D(MeshSharedPtr m) : ProcessModule(m)
m_config["blcurves"] =
ConfigOption(false, "0", "Generate parallelograms on these curves");
m_config["blthick"] =
ConfigOption(false, "0", "Parallelogram layer thickness");
ConfigOption(false, "0.0", "Parallelogram layer thickness");
}
Generator2D::~Generator2D()
......@@ -71,6 +71,9 @@ void Generator2D::Process()
m_mesh->m_numNodes = m_mesh->m_cad->GetNumVerts();
m_thickness_ID =
m_thickness.DefineFunction("x y z", m_config["blthick"].as<string>());
// linear mesh all curves
for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
{
......@@ -80,12 +83,25 @@ void Generator2D::Process()
"Curve progress");
}
m_curvemeshes[i] =
MemoryManager<CurveMesh>::AllocateSharedPtr(i, m_mesh);
vector<unsigned int>::iterator f =
find(m_blCurves.begin(), m_blCurves.end(), i);
if (f == m_blCurves.end())
{
m_curvemeshes[i] =
MemoryManager<CurveMesh>::AllocateSharedPtr(i, m_mesh);
}
else
{
m_curvemeshes[i] = MemoryManager<CurveMesh>::AllocateSharedPtr(
i, m_mesh, m_config["blthick"].as<string>());
}
m_curvemeshes[i]->Mesh();
}
////////////////////////////////////////
EdgeSet::iterator it;
for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end(); it++)
{
......@@ -108,6 +124,14 @@ void Generator2D::Process()
m_mesh->m_element[1].push_back(E2);
}
for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
{
m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
i, m_mesh, m_curvemeshes, 100);
m_facemeshes[i]->OrientateCurves();
}
if (m_config["blcurves"].beenSet)
{
// we need to do the boundary layer generation in a face by face basis
......@@ -120,10 +144,6 @@ void Generator2D::Process()
// orienation
for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
{
m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
i, m_mesh, m_curvemeshes, 100);
m_facemeshes[i]->OrientateCurves();
MakeBL(i, m_facemeshes[i]->GetEdges());
}
}
......@@ -136,7 +156,9 @@ void Generator2D::Process()
}
// linear mesh all surfaces
for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
map<int, FaceMeshSharedPtr>::iterator fit;
int i = 1;
for (fit = m_facemeshes.begin(); fit != m_facemeshes.end(); fit++)
{
if (m_mesh->m_verbose)
{
......@@ -144,10 +166,16 @@ void Generator2D::Process()
"Face progress");
}
m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
i, m_mesh, m_curvemeshes, 100);
if (m_config["blcurves"].beenSet)
{
// for bl surfaces orientate curves needs to be run again to
// push nodes to the edges of the system
fit->second->ResetCurvemeshes(m_curvemeshes);
fit->second->OrientateCurves();
}
m_facemeshes[i]->Mesh();
fit->second->Mesh();
i++;
}
ProcessVertices();
......@@ -169,8 +197,6 @@ void Generator2D::MakeBLPrep()
// identify the nodes which will become the boundary layer.
ParseUtils::GenerateSeqVector(m_config["blcurves"].as<string>().c_str(),
m_blCurves);
m_thickness_ID =
m_thickness.DefineFunction("x y z", m_config["blthick"].as<string>());
for (vector<unsigned>::iterator it = m_blCurves.begin();
it != m_blCurves.end(); ++it)
......@@ -184,14 +210,14 @@ void Generator2D::MakeBLPrep()
}
}
void Generator2D::MakeBL(int faceid, vector<EdgeLoop> e)
void Generator2D::MakeBL(int faceid, vector<EdgeLoopSharedPtr> e)
{
map<int, int> edgeToOrient;
for (vector<EdgeLoop>::iterator lit = e.begin(); lit != e.end(); ++lit)
for (int j = 0; j < e.size(); j++)
{
for (int i = 0; i < lit->edges.size(); ++i)
for (int i = 0; i < e[j]->edges.size(); ++i)
{
edgeToOrient[lit->edges[i]->GetId()] = lit->edgeo[i];
edgeToOrient[e[j]->edges[i]->GetId()] = e[j]->edgeo[i];
}
}
......
......@@ -69,7 +69,7 @@ private:
void MakeBLPrep();
void MakeBL(int faceid, std::vector<EdgeLoop> e);
void MakeBL(int faceid, std::vector<EdgeLoopSharedPtr> e);
void Report();
/// map of individual surface meshes from parametric surfaces
......
......@@ -72,7 +72,7 @@ public:
/**
* @brief Get the loop structures which bound the cad surface
*/
std::vector<EdgeLoop> GetEdges()
std::vector<EdgeLoopSharedPtr> GetEdges()
{
return m_edges;
}
......@@ -164,7 +164,7 @@ protected:
/// normal
bool m_correctNormal;
/// List of bounding edges in loops with orientation.
std::vector<EdgeLoop> m_edges;
std::vector<EdgeLoopSharedPtr> m_edges;
/// Function which tests the the value of uv used is within the surface
virtual void Test(Array<OneD, NekDouble> uv) = 0;
......
......@@ -70,6 +70,8 @@ struct EdgeLoop
NekDouble area;
};
typedef boost::shared_ptr<EdgeLoop> EdgeLoopSharedPtr;
/**
* @brief Base class for CAD interface system.
*
......
......@@ -45,7 +45,8 @@ namespace NekMeshUtils
std::string CADSurfOCE::key = GetCADSurfFactory().RegisterCreatorFunction(
"oce", CADSurfOCE::create, "CADSurfOCE");
void CADSurfOCE::Initialise(int i, TopoDS_Shape in, vector<EdgeLoop> ein)
void CADSurfOCE::Initialise(int i, TopoDS_Shape in,
vector<EdgeLoopSharedPtr> ein)
{
// this bit of code changes the units of the cad from mm opencascade
// defualt to m
......
......@@ -63,7 +63,7 @@ public:
{
}
void Initialise(int i, TopoDS_Shape in, std::vector<EdgeLoop> ein);
void Initialise(int i, TopoDS_Shape in, std::vector<EdgeLoopSharedPtr> ein);
virtual Array<OneD, NekDouble> GetBounds();
virtual Array<OneD, NekDouble> N (Array<OneD, NekDouble> uv);
......
......@@ -190,12 +190,12 @@ bool CADSystemOCE::LoadCAD()
}
}
vector<EdgeLoop> edgeloops;
vector<EdgeLoopSharedPtr> edgeloops;
// now we acutally analyse the loops for cad building
for (int j = 1; j <= mapOfWires.Extent(); j++)
{
EdgeLoop edgeloop;
EdgeLoopSharedPtr edgeloop = boost::shared_ptr<EdgeLoop>(new EdgeLoop);
TopoDS_Shape wire = mapOfWires.FindKey(j);
......@@ -217,7 +217,7 @@ bool CADSystemOCE::LoadCAD()
Array<OneD, NekDouble> cen(2);
cen[0] = p2.X();
cen[1] = p2.Y();
edgeloop.center = cen;
edgeloop->center = cen;
BRepTools_WireExplorer exp;
......@@ -230,8 +230,8 @@ bool CADSystemOCE::LoadCAD()
if (mapOfEdges.Contains(edge))
{
int e = mapOfEdges.FindIndex(edge);
edgeloop.edges.push_back(m_curves[e]);
edgeloop.edgeo.push_back(exp.Orientation());
edgeloop->edges.push_back(m_curves[e]);
edgeloop->edgeo.push_back(exp.Orientation());
adjsurfmap[e].push_back(i);
}
......@@ -310,7 +310,8 @@ void CADSystemOCE::AddCurve(int i, TopoDS_Shape in, int fv, int lv)
m_curves[i]->SetVert(vs);
}
void CADSystemOCE::AddSurf(int i, TopoDS_Shape in, vector<EdgeLoop> ein)
void CADSystemOCE::AddSurf(int i, TopoDS_Shape in,
vector<EdgeLoopSharedPtr> ein)
{
CADSurfSharedPtr newSurf = GetCADSurfFactory().CreateInstance(key);
boost::static_pointer_cast<CADSurfOCE>(newSurf)->Initialise(i, in, ein);
......@@ -324,7 +325,7 @@ void CADSystemOCE::AddSurf(int i, TopoDS_Shape in, vector<EdgeLoop> ein)
int tote = 0;
for (int i = 0; i < ein.size(); i++)
{
tote += ein[i].edges.size();
tote += ein[i]->edges.size();
}
ASSERTL0(tote != 1, "cannot handle periodic curves");
......
......@@ -71,7 +71,7 @@ private:
/// Function to add curve to CADSystem::m_curves.
void AddCurve(int i, TopoDS_Shape in, int fv, int lv);
/// Function to add surface to CADSystem::m_surfs.
void AddSurf(int i, TopoDS_Shape in, std::vector<EdgeLoop> ein);
void AddSurf(int i, TopoDS_Shape in, std::vector<EdgeLoopSharedPtr> ein);
TopoDS_Shape BuildNACA(std::string naca);
/// OCC master object
......
......@@ -101,11 +101,11 @@ ADD_DEFINITIONS(-DNEKMESHUTILS_EXPORTS)
ADD_NEKTAR_LIBRARY(NekMeshUtils lib ${NEKTAR_LIBRARY_TYPE}
${NEKMESHUTILS_SOURCES} ${NEKMESHUTILS_HEADERS})
TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PUBLIC LocalRegions)
TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PUBLIC SpatialDomains)
IF(NEKTAR_USE_MESHGEN)
TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PRIVATE ${TETGEN_LIBRARY})
TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PUBLIC ${OCC_LIBRARIES})
TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PRIVATE ${OCC_LIBRARIES})
ADD_DEPENDENCIES(NekMeshUtils oce-0.17 tetgen-1.5)
......
......@@ -100,6 +100,12 @@ void Octree::Process()
}
}
NekDouble Octree::QueryR(Array<OneD, NekDouble> loc)
{
NekDouble d = Query(loc);
return d/2.0/(sqrt(m_eps * (2.0 - m_eps)));
}
NekDouble Octree::Query(Array<OneD, NekDouble> loc)
{
// starting at master octant 0 move through succsesive m_octants which
......
......@@ -127,6 +127,8 @@ public:
*/
NekDouble Query(Array<OneD, NekDouble> loc);
NekDouble QueryR(Array<OneD, NekDouble> loc);
/**
* @brief returns the miminum spacing in the octree (for meshing purposes)
*
......
......@@ -33,8 +33,8 @@
//
////////////////////////////////////////////////////////////////////////////////
#include <NekMeshUtils/SurfaceMeshing/CurveMesh.h>
#include <NekMeshUtils/Octree/Octree.h>
#include <NekMeshUtils/SurfaceMeshing/CurveMesh.h>
using namespace std;
namespace Nektar
......@@ -44,12 +44,13 @@ namespace NekMeshUtils
void CurveMesh::Mesh()
{
//this algorithm is mostly based on the work in chapter 19
// this algorithm is mostly based on the work in chapter 19
m_bounds = m_cadcurve->Bounds();
m_curvelength = m_cadcurve->GetTotLength();
m_numSamplePoints = int(m_curvelength / m_mesh->m_octree->GetMinDelta()) + 5;
ds = m_curvelength / (m_numSamplePoints - 1);
m_bounds = m_cadcurve->Bounds();
m_curvelength = m_cadcurve->GetTotLength();
m_numSamplePoints =
int(m_curvelength / m_mesh->m_octree->GetMinDelta()) + 5;
ds = m_curvelength / (m_numSamplePoints - 1);
GetSampleFunction();
......@@ -90,7 +91,7 @@ void CurveMesh::Mesh()
iterationcounter++;
NekDouble rhs = EvaluateDS(ski) / Ae * (EvaluatePS(ski) - k);
lastSki = ski;
ski = ski - rhs;
ski = ski - rhs;
if (abs(lastSki - ski) < 1E-8)
{
iterate = false;
......@@ -107,11 +108,11 @@ void CurveMesh::Mesh()
Array<OneD, NekDouble> loc;
vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();
vector<CADSurfSharedPtr> s = m_cadcurve->GetAdjSurf();
vector<CADSurfSharedPtr> s = m_cadcurve->GetAdjSurf();
ASSERTL0(s.size() == (m_mesh->m_cad->Is2D()) ? 1 : 2, "invalid curve");
NodeSharedPtr n = verts[0]->GetNode();
t = m_bounds[0];
t = m_bounds[0];
n->SetCADCurve(m_id, m_cadcurve, t);
loc = n->GetLoc();
for (int j = 0; j < s.size(); j++)
......@@ -172,9 +173,10 @@ void CurveMesh::Mesh()
m_meshedges.push_back(e);
}
if(m_mesh->m_verbose)
if (m_mesh->m_verbose)
{
cout << "\r "
cout << "\r "
" "
" ";
cout << scientific << "\r\t\tCurve " << m_id << endl
<< "\t\t\tLength: " << m_curvelength << endl
......@@ -309,7 +311,19 @@ void CurveMesh::GetSampleFunction()
loc = m_cadcurve->P(t);
dsti[0] = m_mesh->m_octree->Query(loc);
NekDouble ts =
m_bl.Evaluate(m_blID, loc[0], loc[1], loc[2], 0.0);
if (ts > 0.0)
{
NekDouble R = m_mesh->m_octree->QueryR(loc);
dsti[0] = R / (R + ts) * m_mesh->m_octree->Query(loc);
}
else
{
dsti[0] = m_mesh->m_octree->Query(loc);
}
dsti[2] = t;
m_dst[i] = dsti;
......
......@@ -38,11 +38,12 @@
#include <boost/shared_ptr.hpp>
#include <NekMeshUtils/MeshElements/Mesh.h>
#include <NekMeshUtils/CADSystem/CADVert.h>
#include <NekMeshUtils/CADSystem/CADCurve.h>
#include <NekMeshUtils/CADSystem/CADVert.h>
#include <NekMeshUtils/MeshElements/Mesh.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Interpreter/AnalyticExpressionEvaluator.hpp>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
/**
......@@ -61,9 +62,9 @@ public:
/**
* @brief default constructor
*/
CurveMesh(int id, MeshSharedPtr m)
: m_id(id), m_mesh(m)
CurveMesh(int id, MeshSharedPtr m) : m_id(id), m_mesh(m)
{
m_blID = m_bl.DefineFunction("x y z", "0.0");
m_cadcurve = m_mesh->m_cad->GetCurve(m_id);
}
......@@ -73,6 +74,17 @@ public:
CurveMesh(int id, MeshSharedPtr m, std::vector<NodeSharedPtr> n)
: m_id(id), m_mesh(m), m_meshpoints(n)
{
m_blID = m_bl.DefineFunction("x y z", "0.0");
m_cadcurve = m_mesh->m_cad->GetCurve(m_id);
}
/**
* @brief alternative constructor which use the octree in a different way
*/
CurveMesh(int id, MeshSharedPtr m, std::string expr)
: m_id(id), m_mesh(m)
{
m_blID = m_bl.DefineFunction("x y z", expr);
m_cadcurve = m_mesh->m_cad->GetCurve(m_id);
}
......@@ -176,6 +188,8 @@ private:
MeshSharedPtr m_mesh;
/// ids of the mesh nodes
std::vector<NodeSharedPtr> m_meshpoints;
LibUtilities::AnalyticExpressionEvaluator m_bl;
int m_blID;
};
typedef boost::shared_ptr<CurveMesh> CurveMeshSharedPtr;
......
......@@ -49,9 +49,9 @@ bool FaceMesh::ValidateCurves()
vector<int> curvesInSurface;
for(int i = 0; i < m_edgeloops.size(); i++)
{