Commit b0dc56a3 authored by Dave Moxey's avatar Dave Moxey

Merge branch 'feature/madrid-tutorial' into 'master'

meshing update for madrid

This updates a number of the 2D meshing routines to be more useful for the tutorial in madrid. I.e removing unhelpful IO
Also has some minor updates to the documentation.

See merge request !726
parents 689f5dd8 4ce0d8f5
...@@ -24,16 +24,16 @@ if(NOT DEFINED OCE_DIR) ...@@ -24,16 +24,16 @@ if(NOT DEFINED OCE_DIR)
# Check for OSX needs to come first because UNIX evaluates to true on OSX # Check for OSX needs to come first because UNIX evaluates to true on OSX
if(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") if(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
if(DEFINED MACPORTS_PREFIX) if(DEFINED MACPORTS_PREFIX)
find_package(OCE 0.17 QUIET HINTS ${MACPORTS_PREFIX}/Library/Frameworks) find_package(OCE 0.15 QUIET HINTS ${MACPORTS_PREFIX}/Library/Frameworks)
elseif(DEFINED HOMEBREW_PREFIX) elseif(DEFINED HOMEBREW_PREFIX)
find_package(OCE 0.17 QUIET HINTS ${HOMEBREW_PREFIX}/Cellar/oce/*) find_package(OCE 0.15 QUIET HINTS ${HOMEBREW_PREFIX}/Cellar/oce/*)
endif() endif()
elseif(UNIX) elseif(UNIX)
set(OCE_DIR "/usr/local/share/cmake/") set(OCE_DIR "/usr/local/share/cmake/")
endif() endif()
endif() endif()
find_package(OCE 0.17 QUIET) find_package(OCE 0.15 QUIET)
if(OCE_FOUND) if(OCE_FOUND)
message(STATUS "-- OpenCASCADE Community Edition has been found.") message(STATUS "-- OpenCASCADE Community Edition has been found.")
set(OCC_INCLUDE_DIR ${OCE_INCLUDE_DIRS}) set(OCC_INCLUDE_DIR ${OCE_INCLUDE_DIRS})
......
...@@ -83,7 +83,7 @@ ...@@ -83,7 +83,7 @@
pages = {293-301}, pages = {293-301},
year = {1996}, year = {1996},
} }
@article{CoRaNa98, @article{CoRaNa98,
author = {M. Courtemanche\, R. J. Ramirez and S. Nattel}, author = {M. Courtemanche\, R. J. Ramirez and S. Nattel},
title = {Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model}, title = {Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model},
...@@ -103,7 +103,7 @@ ...@@ -103,7 +103,7 @@
pages = {1501-1526}, pages = {1501-1526},
year = {1991}, year = {1991},
} }
@article{TuPa06, @article{TuPa06,
author = {K. H. W. J. ten Tusscher and A. V. Panfilov}, author = {K. H. W. J. ten Tusscher and A. V. Panfilov},
title = {Alternans and spiral breakup in a human ventricular tissue model}, title = {Alternans and spiral breakup in a human ventricular tissue model},
...@@ -123,7 +123,7 @@ ...@@ -123,7 +123,7 @@
pages = {4331-51}, pages = {4331-51},
year = {2011}, year = {2011},
} }
@article{ShKa96, @article{ShKa96,
title={Tetrahedral< i> hp</i> Finite Elements: Algorithms and Flow Simulations}, title={Tetrahedral< i> hp</i> Finite Elements: Algorithms and Flow Simulations},
author={Sherwin, SJ and Karniadakis, G Em}, author={Sherwin, SJ and Karniadakis, G Em},
...@@ -378,9 +378,9 @@ year={2011} ...@@ -378,9 +378,9 @@ year={2011}
} }
@article{GuSh03, @article{GuSh03,
Author="J.L. Guermond and J. Shen", Author="J.L. Guermond and J. Shen",
title="Velocity-correction projection methods for incompressible flows", title="Velocity-correction projection methods for incompressible flows",
journal="SIAM J. Numer.\ Anal.", journal="SIAM J. Numer.\ Anal.",
volume=41, volume=41,
pages = "112--134", pages = "112--134",
year=2003 year=2003
...@@ -460,3 +460,18 @@ year={2011} ...@@ -460,3 +460,18 @@ year={2011}
pages = {1079-1097}, pages = {1079-1097},
} }
@inproceedings{TuPeMo16,
abstract = {The generation of sufficiently high quality unstructured high-order meshes remains a significant obstacle in the adoption of high-order methods. However, there is little consensus on which approach is the most robust, fastest and produces the 'best' meshes. We aim to provide a route to investigate this question, by examining popular high-order mesh generation methods in the context of an efficient variational framework for the generation of curvilinear meshes. By considering previous works in a variational form, we are able to compare their characteristics and study their robustness. Alongside a description of the theory and practical implementation details, including an efficient multi-threading parallelisation strategy, we demonstrate the effectiveness of the framework, showing how it can be used for both mesh quality optimisation and untangling of invalid meshes.},
author = {Turner, M and Peir{\'{o}}, J and Moxey, D},
booktitle = {25th International Meshing Roundtable},
doi = {10.1016/j.proeng.2016.11.069},
file = {:Users/mike/Downloads/1-s2.0-S1877705816333781-main.pdf:pdf},
issn = {18777058},
keywords = {energy functional,high-order mesh generation,numerical optimization,variational mesh generation},
pages = {340--352},
title = {{A Variational Framework for High-Order Mesh Generation}},
url = {www.elsevier.com/locate/procedia%5Cnhttp://linkinghub.elsevier.com/retrieve/pii/S1877705816333781},
volume = {163},
year = {2016}
}
This diff is collapsed.
This diff is collapsed.
...@@ -121,7 +121,7 @@ void Generator2D::Process() ...@@ -121,7 +121,7 @@ void Generator2D::Process()
for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++) for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
{ {
m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr( m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
i, m_mesh, m_curvemeshes, m_mesh->m_cad->GetNumSurf() > 100); i, m_mesh, m_curvemeshes, 100);
m_facemeshes[i]->OrientateCurves(); m_facemeshes[i]->OrientateCurves();
MakeBL(i, m_facemeshes[i]->GetEdges()); MakeBL(i, m_facemeshes[i]->GetEdges());
...@@ -145,7 +145,7 @@ void Generator2D::Process() ...@@ -145,7 +145,7 @@ void Generator2D::Process()
} }
m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr( m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
i, m_mesh, m_curvemeshes, m_mesh->m_cad->GetNumSurf() > 100); i, m_mesh, m_curvemeshes, 100);
m_facemeshes[i]->Mesh(); m_facemeshes[i]->Mesh();
} }
...@@ -282,7 +282,7 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoop> e) ...@@ -282,7 +282,7 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoop> e)
ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false); ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
vector<int> tags; vector<int> tags;
tags.push_back(102); tags.push_back(101);
ElementSharedPtr E = GetElementFactory().CreateInstance( ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::eQuadrilateral, conf, qns, tags); LibUtilities::eQuadrilateral, conf, qns, tags);
......
...@@ -111,6 +111,11 @@ public: ...@@ -111,6 +111,11 @@ public:
return m_2d; return m_2d;
} }
void SetNACA(std::string i)
{
m_naca = i;
}
/** /**
* @brief Initialises CAD and makes surface, curve and vertex maps. * @brief Initialises CAD and makes surface, curve and vertex maps.
* *
...@@ -203,6 +208,7 @@ protected: ...@@ -203,6 +208,7 @@ protected:
std::map<int, CADVertSharedPtr> m_verts; std::map<int, CADVertSharedPtr> m_verts;
bool m_2d; bool m_2d;
std::string m_naca;
}; };
typedef boost::shared_ptr<CADSystem> CADSystemSharedPtr; typedef boost::shared_ptr<CADSystem> CADSystemSharedPtr;
......
...@@ -32,6 +32,7 @@ ...@@ -32,6 +32,7 @@
// Description: cad object methods. // Description: cad object methods.
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/BasicUtils/ParseUtils.hpp>
#include <NekMeshUtils/CADSystem/OCE/CADSystemOCE.h> #include <NekMeshUtils/CADSystem/OCE/CADSystemOCE.h>
#include <NekMeshUtils/CADSystem/OCE/CADVertOCE.h> #include <NekMeshUtils/CADSystem/OCE/CADVertOCE.h>
...@@ -51,8 +52,9 @@ std::string CADSystemOCE::key = GetEngineFactory().RegisterCreatorFunction( ...@@ -51,8 +52,9 @@ std::string CADSystemOCE::key = GetEngineFactory().RegisterCreatorFunction(
bool CADSystemOCE::LoadCAD() bool CADSystemOCE::LoadCAD()
{ {
if (m_name.find('.') != std::string::npos) if (m_naca.size() == 0)
{ {
//not a naca profile behave normally
// Takes step file and makes OpenCascade shape // Takes step file and makes OpenCascade shape
STEPControl_Reader reader; STEPControl_Reader reader;
reader = STEPControl_Reader(); reader = STEPControl_Reader();
...@@ -67,8 +69,6 @@ bool CADSystemOCE::LoadCAD() ...@@ -67,8 +69,6 @@ bool CADSystemOCE::LoadCAD()
} }
else else
{ {
cout << m_name << " is not a STEP file, assuming it is "
<< "a 4 digit NACA code" << endl;
shape = BuildNACA(m_name); shape = BuildNACA(m_name);
} }
...@@ -361,6 +361,9 @@ Array<OneD, NekDouble> CADSystemOCE::GetBoundingBox() ...@@ -361,6 +361,9 @@ Array<OneD, NekDouble> CADSystemOCE::GetBoundingBox()
TopoDS_Shape CADSystemOCE::BuildNACA(string naca) TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
{ {
ASSERTL0(naca.length() == 4, "not a 4 digit code"); ASSERTL0(naca.length() == 4, "not a 4 digit code");
vector<NekDouble> data;
ParseUtils::GenerateUnOrderedVector(m_naca.c_str(), data);
ASSERTL0(data.size() == 5, "not a vaild domain");
int n = boost::lexical_cast<int>(naca); int n = boost::lexical_cast<int>(naca);
NekDouble T = (n%100) / 100.0; NekDouble T = (n%100) / 100.0;
...@@ -442,17 +445,26 @@ TopoDS_Shape CADSystemOCE::BuildNACA(string naca) ...@@ -442,17 +445,26 @@ TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
BRepBuilderAPI_MakeWire aeroWireBuilder(aeroEdge, TeEdge); BRepBuilderAPI_MakeWire aeroWireBuilder(aeroEdge, TeEdge);
TopoDS_Wire aeroWire = aeroWireBuilder.Wire(); TopoDS_Wire aeroWire = aeroWireBuilder.Wire();
BRepBuilderAPI_MakeEdge domInlBuilder(gp_Pnt(-2000.0,-2000.0,0.0), gp_Trsf transform;
gp_Pnt(-2000.0,2000.0,0.0)); gp_Ax1 rotAx(gp_Pnt(500.0,0.0,0.0),gp_Dir(gp_Vec(0.0,0.0,-1.0)));
transform.SetRotation(rotAx, data[4]/180.0*M_PI);
TopLoc_Location mv(transform);
aeroWire.Move(mv);
BRepBuilderAPI_MakeEdge domInlBuilder(gp_Pnt(data[0]*1000.0,data[1]*1000.0,0.0),
gp_Pnt(data[0]*1000.0,data[3]*1000.0,0.0));
TopoDS_Edge inlEdge = domInlBuilder.Edge(); TopoDS_Edge inlEdge = domInlBuilder.Edge();
BRepBuilderAPI_MakeEdge domTopBuilder(gp_Pnt(-2000.0,2000.0,0.0),
gp_Pnt(5000.0,2000.0,0.0)); BRepBuilderAPI_MakeEdge domTopBuilder(gp_Pnt(data[0]*1000.0,data[3]*1000.0,0.0),
gp_Pnt(data[2]*1000.0,data[3]*1000.0,0.0));
TopoDS_Edge topEdge = domTopBuilder.Edge(); TopoDS_Edge topEdge = domTopBuilder.Edge();
BRepBuilderAPI_MakeEdge domOutBuilder(gp_Pnt(5000.0,2000.0,0.0),
gp_Pnt(5000.0,-2000.0,0.0)); BRepBuilderAPI_MakeEdge domOutBuilder(gp_Pnt(data[2]*1000.0,data[3]*1000.0,0.0),
gp_Pnt(data[2]*1000.0,data[1]*1000.0,0.0));
TopoDS_Edge outEdge = domOutBuilder.Edge(); TopoDS_Edge outEdge = domOutBuilder.Edge();
BRepBuilderAPI_MakeEdge domBotBuilder(gp_Pnt(5000.0,-2000.0,0.0),
gp_Pnt(-2000.0,-2000.0,0.0)); BRepBuilderAPI_MakeEdge domBotBuilder(gp_Pnt(data[2]*1000.0,data[1]*1000.0,0.0),
gp_Pnt(data[0]*1000.0,data[1]*1000.0,0.0));
TopoDS_Edge botEdge = domBotBuilder.Edge(); TopoDS_Edge botEdge = domBotBuilder.Edge();
BRepBuilderAPI_MakeWire domWireBuilder(inlEdge, topEdge, outEdge, botEdge); BRepBuilderAPI_MakeWire domWireBuilder(inlEdge, topEdge, outEdge, botEdge);
......
...@@ -78,5 +78,6 @@ ...@@ -78,5 +78,6 @@
#include <BRepBuilderAPI_MakeWire.hxx> #include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepBuilderAPI_MakeFace.hxx> #include <BRepBuilderAPI_MakeFace.hxx>
#include <STEPControl_Writer.hxx> #include <STEPControl_Writer.hxx>
#include <gp_Ax1.hxx>
#endif #endif
...@@ -53,6 +53,8 @@ ProcessLoadCAD::ProcessLoadCAD(MeshSharedPtr m) : ProcessModule(m) ...@@ -53,6 +53,8 @@ ProcessLoadCAD::ProcessLoadCAD(MeshSharedPtr m) : ProcessModule(m)
ConfigOption(false, "", "Generate prisms on these surfs"); ConfigOption(false, "", "Generate prisms on these surfs");
m_config["2D"] = m_config["2D"] =
ConfigOption(true, "", "allow 2d loading"); ConfigOption(true, "", "allow 2d loading");
m_config["NACA"] =
ConfigOption(false, "", "naca domain");
} }
ProcessLoadCAD::~ProcessLoadCAD() ProcessLoadCAD::~ProcessLoadCAD()
...@@ -75,6 +77,11 @@ void ProcessLoadCAD::Process() ...@@ -75,6 +77,11 @@ void ProcessLoadCAD::Process()
m_mesh->m_cad->Set2D(); m_mesh->m_cad->Set2D();
} }
if(m_config["NACA"].beenSet)
{
m_mesh->m_cad->SetNACA(m_config["NACA"].as<string>());
}
ASSERTL0(m_mesh->m_cad->LoadCAD(), "Failed to load CAD"); ASSERTL0(m_mesh->m_cad->LoadCAD(), "Failed to load CAD");
if (m_mesh->m_verbose) if (m_mesh->m_verbose)
......
...@@ -66,16 +66,14 @@ Line::Line(ElmtConfig pConf, ...@@ -66,16 +66,14 @@ Line::Line(ElmtConfig pConf,
{ {
m_vertex.push_back(pNodeList[i]); m_vertex.push_back(pNodeList[i]);
} }
vector<NodeSharedPtr> edgeNodes;
if (m_conf.m_order > 1) if (m_conf.m_order > 1)
{ {
for (int j = 0; j < n; ++j) for (int j = 0; j < n; ++j)
{ {
edgeNodes.push_back(pNodeList[2 + j]); m_volumeNodes.push_back(pNodeList[2 + j]);
} }
} }
m_edge.push_back(boost::shared_ptr<Edge>(new Edge(
pNodeList[0], pNodeList[1], edgeNodes, m_conf.m_edgeCurveType)));
} }
SpatialDomains::GeometrySharedPtr Line::GetGeom(int coordDim) SpatialDomains::GeometrySharedPtr Line::GetGeom(int coordDim)
......
...@@ -104,6 +104,18 @@ NekDouble Octree::Query(Array<OneD, NekDouble> loc) ...@@ -104,6 +104,18 @@ NekDouble Octree::Query(Array<OneD, NekDouble> loc)
{ {
// starting at master octant 0 move through succsesive m_octants which // starting at master octant 0 move through succsesive m_octants which
// contain the point loc until a leaf is found // contain the point loc until a leaf is found
//first search through sourcepoints
NekDouble tmp = numeric_limits<double>::max();
for(int i = 0; i < m_lsources.size(); i++)
{
if(m_lsources[i].withinRange(loc))
{
tmp = min(m_lsources[i].delta,tmp);
}
}
OctantSharedPtr n = m_masteroct; OctantSharedPtr n = m_masteroct;
int quad; int quad;
...@@ -173,7 +185,7 @@ NekDouble Octree::Query(Array<OneD, NekDouble> loc) ...@@ -173,7 +185,7 @@ NekDouble Octree::Query(Array<OneD, NekDouble> loc)
found = true; found = true;
} }
} }
return n->GetDelta(); return min(n->GetDelta(),tmp);
} }
void Octree::WriteOctree(string nm) void Octree::WriteOctree(string nm)
...@@ -468,7 +480,7 @@ void Octree::SmoothSurfaceOctants() ...@@ -468,7 +480,7 @@ void Octree::SmoothSurfaceOctants()
{ {
if (it->second[j]->IsDeltaKnown() && if (it->second[j]->IsDeltaKnown() &&
it->second[j]->GetDelta() < oct->GetDelta() && it->second[j]->GetDelta() < oct->GetDelta() &&
ddx(oct, it->second[j]) > 0.1) ddx(oct, it->second[j]) > 0.2)
{ {
check.push_back(it->second[j]); check.push_back(it->second[j]);
} }
...@@ -485,9 +497,9 @@ void Octree::SmoothSurfaceOctants() ...@@ -485,9 +497,9 @@ void Octree::SmoothSurfaceOctants()
{ {
NekDouble r = oct->Distance(check[j]); NekDouble r = oct->Distance(check[j]);
if (0.099 * r + check[j]->GetDelta() < deltaSM) if (0.199 * r + check[j]->GetDelta() < deltaSM)
{ {
deltaSM = 0.099 * r + check[j]->GetDelta(); deltaSM = 0.199 * r + check[j]->GetDelta();
} }
} }
oct->SetDelta(deltaSM); oct->SetDelta(deltaSM);
...@@ -538,9 +550,9 @@ void Octree::PropagateDomain() ...@@ -538,9 +550,9 @@ void Octree::PropagateDomain()
{ {
NekDouble r = oct->Distance(known[j]); NekDouble r = oct->Distance(known[j]);
if (0.14 * r + known[j]->GetDelta() < m_maxDelta) if (0.199 * r + known[j]->GetDelta() < m_maxDelta)
{ {
deltaPrime.push_back(0.14 * r + deltaPrime.push_back(0.199 * r +
known[j]->GetDelta()); known[j]->GetDelta());
} }
else else
...@@ -817,59 +829,6 @@ int Octree::CountElemt() ...@@ -817,59 +829,6 @@ int Octree::CountElemt()
return int(total); return int(total);
} }
//struct to assist in the creation of linesources in the code
struct linesource
{
Array<OneD, NekDouble> x1, x2;
NekDouble R, delta;
linesource(Array<OneD, NekDouble> p1,
Array<OneD, NekDouble> p2,
NekDouble r,
NekDouble d)
: x1(p1), x2(p2), R(r), delta(d)
{
}
bool withinRange(Array<OneD, NekDouble> p)
{
Array<OneD, NekDouble> Le(3), Re(3), s(3);
for (int i = 0; i < 3; i++)
{
Le[i] = p[i] - x1[i];
Re[i] = p[i] - x2[i];
s[i] = x2[i] - x1[i];
}
Array<OneD, NekDouble> dev(3);
dev[0] = Le[1] * Re[2] - Re[1] * Le[2];
dev[1] = Le[0] * Re[2] - Re[0] * Le[2];
dev[2] = Le[0] * Re[1] - Re[0] * Le[1];
NekDouble dist =
sqrt(dev[0] * dev[0] + dev[1] * dev[1] + dev[2] * dev[2]) /
sqrt(s[0] * s[0] + s[1] * s[1] + s[2] * s[2]);
NekDouble t = -1.0 * ((x1[0] - p[0]) * s[0] + (x1[1] - p[1]) * s[1] +
(x1[1] - p[1]) * s[1]) /
Length() / Length();
if (dist < R && !(t > 1) && !(t < 0))
{
return true;
}
else
{
return false;
}
}
NekDouble Length()
{
return sqrt((x1[0] - x2[0]) * (x1[0] - x2[0]) +
(x1[1] - x2[1]) * (x1[1] - x2[1]) +
(x1[2] - x2[2]) * (x1[2] - x2[2]));
}
};
void Octree::CompileSourcePointList() void Octree::CompileSourcePointList()
{ {
...@@ -991,8 +950,8 @@ void Octree::CompileSourcePointList() ...@@ -991,8 +950,8 @@ void Octree::CompileSourcePointList()
// these are the acutal number of sample points in each parametric // these are the acutal number of sample points in each parametric
// direction // direction
int nu = ceil(DeltaU / m_minDelta) * 40; int nu = ceil(DeltaU / m_minDelta) * 2;
int nv = ceil(DeltaV / m_minDelta) * 40; int nv = ceil(DeltaV / m_minDelta) * 2;
for (int j = 0; j < nu; j++) for (int j = 0; j < nu; j++)
{ {
...@@ -1057,7 +1016,6 @@ void Octree::CompileSourcePointList() ...@@ -1057,7 +1016,6 @@ void Octree::CompileSourcePointList()
cout << "\t\tModifying based on refinement lines" << endl; cout << "\t\tModifying based on refinement lines" << endl;
} }
// now deal with the user defined spacing // now deal with the user defined spacing
vector<linesource> lsources;
vector<string> lines; vector<string> lines;
boost::split(lines, m_refinement, boost::is_any_of(":")); boost::split(lines, m_refinement, boost::is_any_of(":"));
...@@ -1076,16 +1034,16 @@ void Octree::CompileSourcePointList() ...@@ -1076,16 +1034,16 @@ void Octree::CompileSourcePointList()
x2[1] = data[4]; x2[1] = data[4];
x2[2] = data[5]; x2[2] = data[5];
lsources.push_back(linesource(x1, x2, data[6], data[7])); m_lsources.push_back(linesource(x1, x2, data[6], data[7]));
} }
// this takes any existing sourcepoints within the influence range // this takes any existing sourcepoints within the influence range
// and modifies them // and modifies them
for (int i = 0; i < m_SPList.size(); i++) /*for (int i = 0; i < m_SPList.size(); i++)
{ {
for (int j = 0; j < lsources.size(); j++) for (int j = 0; j < m_lsources.size(); j++)
{ {
if (lsources[j].withinRange(m_SPList[i]->GetLoc())) if (m_lsources[j].withinRange(m_SPList[i]->GetLoc()))
{ {
if(m_SPList[i]->GetType() == ePBoundary) if(m_SPList[i]->GetType() == ePBoundary)
{ {
...@@ -1096,10 +1054,10 @@ void Octree::CompileSourcePointList() ...@@ -1096,10 +1054,10 @@ void Octree::CompileSourcePointList()
m_SPList[i] = bp->ChangeType(); m_SPList[i] = bp->ChangeType();
} }
m_SPList[i]->SetDelta(lsources[j].delta); m_SPList[i]->SetDelta(m_lsources[j].delta);
} }
} }
} }*/
///TODO add extra source points from the line souce to the octree ///TODO add extra source points from the line souce to the octree
} }
} }
......
...@@ -47,6 +47,57 @@ namespace Nektar ...@@ -47,6 +47,57 @@ namespace Nektar
namespace NekMeshUtils namespace NekMeshUtils
{ {
//struct to assist in the creation of linesources in the code
struct linesource
{
Array<OneD, NekDouble> x1, x2;
NekDouble R, delta;
linesource(Array<OneD, NekDouble> p1,
Array<OneD, NekDouble> p2,
NekDouble r,
NekDouble d)
: x1(p1), x2(p2), R(r), delta(d)
{
}
bool withinRange(Array<OneD, NekDouble> p)
{
Array<OneD, NekDouble> Le(3), Re(3), s(3);
for (int i = 0; i < 3; i++)
{
Le[i] = p[i] - x1[i];
Re[i] = p[i] - x2[i];
s[i] = x2[i] - x1[i];
}
Array<OneD, NekDouble> dev(3);
dev[0] = Le[1] * Re[2] - Re[1] * Le[2];
dev[1] = Le[2] * Re[0] - Re[2] * Le[0];
dev[2] = Le[0] * Re[1] - Re[0] * Le[1];
NekDouble dist =
sqrt(dev[0] * dev[0] + dev[1] * dev[1] + dev[2] * dev[2]) / Length();
NekDouble t = -1.0 * ((x1[0] - p[0]) * s[0] + (x1[1] - p[1]) * s[1] +
(x1[2] - p[2]) * s[2]) / Length() / Length();
if (dist < R && !(t > 1) && !(t < 0))
{
return true;
}
else
{
return false;
}
}
NekDouble Length()
{
return sqrt((x1[0] - x2[0]) * (x1[0] - x2[0]) +
(x1[1] - x2[1]) * (x1[1] - x2[1]) +
(x1[2] - x2[2]) * (x1[2] - x2[2]));
}
};
/** /**
* @brief class for octree * @brief class for octree
* *
...@@ -188,6 +239,7 @@ private: ...@@ -188,6 +239,7 @@ private:
MeshSharedPtr m_mesh; MeshSharedPtr m_mesh;
std::string m_refinement; std::string m_refinement;
std::vector<linesource> m_lsources;
}; };
typedef boost::shared_ptr<Octree> OctreeSharedPtr; typedef boost::shared_ptr<Octree> OctreeSharedPtr;
......
...@@ -933,7 +933,12 @@ bool FaceMesh::Validate() ...@@ -933,7 +933,12 @@ bool FaceMesh::Validate()
numValid++; numValid++;
} }
if (numValid != 3) NekDouble rmin = min(r[0],r[1]);
rmin = min(rmin,r[2]);
NekDouble rmax = max(r[0],r[1]);
rmax = min(rmax,r[2]);
if (numValid != 3 || rmax / rmin > 1.1)
{ {
Array<OneD, NekDouble> ainfo, binfo, cinfo; Array<OneD, NekDouble> ainfo, binfo, cinfo;
ainfo = m_connec[i][0]->GetCADSurfInfo(m_id); ainfo = m_connec[i][0]->GetCADSurfInfo(m_id);
......
...@@ -87,12 +87,6 @@ void SurfaceMesh::Process() ...@@ -87,12 +87,6 @@ void SurfaceMesh::Process()
if (m_mesh->m_verbose) if (m_mesh->m_verbose)
cout << endl << "\tFace meshing:" << endl << endl; cout << endl << "\tFace meshing:" << endl << endl;
int prefix = 100;
if(m_mesh->m_cad->GetNumSurf() > 1000)
{
prefix *= 10;
}
bool validError = false; bool validError = false;
for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++) for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
{ {
...@@ -103,7 +97,7 @@ void SurfaceMesh::Process() ...@@ -103,7 +97,7 @@ void SurfaceMesh::Process()
} }
m_facemeshes[i] = m_facemeshes[i] =
MemoryManager<FaceMesh>::