Commit 382be47b authored by Michael Turner's avatar Michael Turner

merge conflict

parents bb7852af e025603b
......@@ -75,10 +75,11 @@ v4.4.0
- Add flag to `insertsurface` process for non-conforming geometries (!700)
- Bug fix to get two meshgen regression tests working (!700)
- Remove libANN in deference to boost::geometry (!703)
- 2D to 3D mesh extrusion module (!715)
- Add a mesh extract option to the linearise module to visualise the result
(!712)
- Refactor library to use NekMesh modules for CAD generation (!704)
- Add a mesh extract option to the linearise module to visualise the result
(!712)
- 2D to 3D mesh extrusion module (!715)
- Add new two-dimensional mesher from NACA code or step file (!720)
**FieldConvert:**
- Move all modules to a new library, FieldUtils, to support post-processing
......
This diff is collapsed.
////////////////////////////////////////////////////////////////////////////////
//
// File: SurfaceMeshing.h
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: class containing all surfacemeshing routines and classes.
//
////////////////////////////////////////////////////////////////////////////////
#ifndef NEKMESHUTILS_2D_2D
#define NEKMESHUTILS_2D_2D
#include <NekMeshUtils/Module/Module.h>
#include <NekMeshUtils/SurfaceMeshing/CurveMesh.h>
#include <NekMeshUtils/SurfaceMeshing/FaceMesh.h>
namespace Nektar
{
namespace NekMeshUtils
{
/**
* @brief class containing all surface meshing routines methods and classes
*/
class Generator2D : public ProcessModule
{
public:
/// Creates an instance of this class
static boost::shared_ptr<Module> create(MeshSharedPtr m)
{
return MemoryManager<Generator2D>::AllocateSharedPtr(m);
}
static ModuleKey className;
Generator2D(MeshSharedPtr m);
virtual ~Generator2D();
virtual void Process();
private:
void MakeBLPrep();
void MakeBL(int faceid, std::vector<EdgeLoop> e);
void Report();
/// map of individual surface meshes from parametric surfaces
std::map<int, FaceMeshSharedPtr> m_facemeshes;
/// map of individual curve meshes of the curves in the domain
std::map<int, CurveMeshSharedPtr> m_curvemeshes;
std::vector<unsigned int> m_blCurves;
NekDouble m_thickness;
std::map<NodeSharedPtr, std::vector<EdgeSharedPtr> > m_nodesToEdge;
};
}
}
#endif
......@@ -97,6 +97,8 @@ public:
*/
virtual Array<OneD, NekDouble> D2(NekDouble t) = 0;
virtual NekDouble Curvature(NekDouble t) = 0;
/**
* @brief Calculates the parametric coordinate and arclength location
* defined by \p s.
......
......@@ -86,6 +86,7 @@ public:
*/
CADSystem(std::string name) : m_name(name)
{
m_2d = false;
}
~CADSystem()
......@@ -100,6 +101,16 @@ public:
return m_name;
}
void Set2D()
{
m_2d = true;
}
bool Is2D()
{
return m_2d;
}
/**
* @brief Initialises CAD and makes surface, curve and vertex maps.
*
......@@ -190,6 +201,8 @@ protected:
std::map<int, CADSurfSharedPtr> m_surfs;
/// Map of vertices
std::map<int, CADVertSharedPtr> m_verts;
bool m_2d;
};
typedef boost::shared_ptr<CADSystem> CADSystemSharedPtr;
......
......@@ -128,6 +128,14 @@ Array<OneD, NekDouble> CADCurveOCE::D2(NekDouble t)
return out;
}
NekDouble CADCurveOCE::Curvature(NekDouble t)
{
GeomLProp_CLProps d(m_c,2,1e-8);
d.SetParameter(t);
return d.Curvature() * 1000.0;
}
Array<OneD, NekDouble> CADCurveOCE::Bounds()
{
Array<OneD, NekDouble> t(2);
......
......@@ -70,6 +70,7 @@ public:
virtual NekDouble tAtArcLength(NekDouble s);
virtual Array<OneD, NekDouble> GetMinMax();
virtual NekDouble loct(Array<OneD, NekDouble> xyz);
virtual NekDouble Curvature(NekDouble t);
void Initialise(int i, TopoDS_Shape in)
{
......
......@@ -51,18 +51,25 @@ std::string CADSystemOCE::key = GetEngineFactory().RegisterCreatorFunction(
bool CADSystemOCE::LoadCAD()
{
cout << "trying " << m_name << endl;
// Takes step file and makes OpenCascade shape
STEPControl_Reader reader;
reader = STEPControl_Reader();
reader.ReadFile(m_name.c_str());
reader.NbRootsForTransfer();
reader.TransferRoots();
shape = reader.OneShape();
if (shape.IsNull())
if (m_name.find('.') != std::string::npos)
{
return false;
// Takes step file and makes OpenCascade shape
STEPControl_Reader reader;
reader = STEPControl_Reader();
reader.ReadFile(m_name.c_str());
reader.NbRootsForTransfer();
reader.TransferRoots();
shape = reader.OneShape();
if (shape.IsNull())
{
return false;
}
}
else
{
cout << m_name << " is not a STEP file, assuming it is "
<< "a 4 digit NACA code" << endl;
shape = BuildNACA(m_name);
}
// faces and verts can be extracted straight from shape
......@@ -105,8 +112,8 @@ bool CADSystemOCE::LoadCAD()
}
}
map<int, vector<int> > adjsurfmap; // from id of curve to list of ids of
// surfs
// from id of curve to list of ids of surfs
map<int, vector<int> > adjsurfmap;
// Adds edges to our type and map
for (int i = 1; i <= mapOfEdges.Extent(); i++)
......@@ -267,7 +274,11 @@ bool CADSystemOCE::LoadCAD()
for (map<int, vector<int> >::iterator it = adjsurfmap.begin();
it != adjsurfmap.end(); it++)
{
ASSERTL0(it->second.size() == 2, "no three curve surfaces");
if(!m_2d)
{
ASSERTL0(it->second.size() == 2, "no three curve surfaces");
}
vector<CADSurfSharedPtr> sfs;
for (int i = 0; i < it->second.size(); i++)
{
......@@ -347,5 +358,111 @@ Array<OneD, NekDouble> CADSystemOCE::GetBoundingBox()
return bound;
}
TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
{
ASSERTL0(naca.length() == 4, "not a 4 digit code");
int n = boost::lexical_cast<int>(naca);
NekDouble T = (n%100) / 100.0;
n/=100;
NekDouble P = (n%10)/10.0;
n/=10;
NekDouble M = (n%10)/100.0;
int np = 25;
Array<OneD, NekDouble> xc(np);
NekDouble dtheta = M_PI/(np-1);
for(int i = 0; i < np; i++)
{
xc[i] = (1.0 - cos(i*dtheta)) / 2.0;
}
Array<OneD, NekDouble> yc(np), dyc(np);
for(int i = 0; i < np; i++)
{
if(xc[i] < P)
{
yc[i] = M / P / P * (2.0 * P * xc[i] - xc[i] * xc[i]);
dyc[i] = 2.0 * M / P / P * (P - xc[i]);
}
else
{
yc[i] = M / (1.0 - P) / (1.0 - P) * (
1.0 - 2.0 * P + 2.0 * P * xc[i] - xc[i] * xc[i]);
dyc[i] = 2.0 * M / (1.0 - P) / (1.0 - P) * (P - xc[i]);
}
}
Array<OneD, NekDouble> yt(np);
for(int i = 0; i < np; i++)
{
yt[i] = T / 0.2 * ( 0.2969 * sqrt(xc[i])
-0.1260 * xc[i]
-0.3516 * xc[i] * xc[i]
+0.2843 * xc[i] * xc[i] * xc[i]
-0.1015 * xc[i] * xc[i] * xc[i] * xc[i]);
}
Array<OneD, NekDouble> x(2*np-1), y(2*np-1);
int l = 0;
for(int i = np-1; i >= 0; i--, l++)
{
NekDouble theta = atan(dyc[i]);
x[l] = xc[i] - yt[i] * sin(theta);
y[l] = yc[i] + yt[i] * cos(theta);
}
for(int i = 1; i < np; i++)
{
NekDouble theta = atan(dyc[i]);
x[i+np-1] = xc[i] + yt[i] * sin(theta);
y[i+np-1] = yc[i] - yt[i] * cos(theta);
}
TColgp_Array1OfPnt pointArray(0,2*np-2);
for(int i = 0; i < 2*np-1; i++)
{
pointArray.SetValue(i,gp_Pnt(x[i]*1000.0,y[i]*1000.0,0.0));
}
GeomAPI_PointsToBSpline spline(pointArray);
Handle(Geom_BSplineCurve) curve = spline.Curve();
BRepBuilderAPI_MakeEdge areoEdgeBuilder(curve);
TopoDS_Edge aeroEdge = areoEdgeBuilder.Edge();
BRepBuilderAPI_MakeEdge aeroTEBuilder(
gp_Pnt(x[0]*1000.0, y[0]*1000.0, 0.0),
gp_Pnt(x[2*np-2]*1000.0, y[2*np-2]*1000.0, 0.0));
TopoDS_Edge TeEdge = aeroTEBuilder.Edge();
BRepBuilderAPI_MakeWire aeroWireBuilder(aeroEdge, TeEdge);
TopoDS_Wire aeroWire = aeroWireBuilder.Wire();
BRepBuilderAPI_MakeEdge domInlBuilder(gp_Pnt(-2000.0,-2000.0,0.0),
gp_Pnt(-2000.0,2000.0,0.0));
TopoDS_Edge inlEdge = domInlBuilder.Edge();
BRepBuilderAPI_MakeEdge domTopBuilder(gp_Pnt(-2000.0,2000.0,0.0),
gp_Pnt(5000.0,2000.0,0.0));
TopoDS_Edge topEdge = domTopBuilder.Edge();
BRepBuilderAPI_MakeEdge domOutBuilder(gp_Pnt(5000.0,2000.0,0.0),
gp_Pnt(5000.0,-2000.0,0.0));
TopoDS_Edge outEdge = domOutBuilder.Edge();
BRepBuilderAPI_MakeEdge domBotBuilder(gp_Pnt(5000.0,-2000.0,0.0),
gp_Pnt(-2000.0,-2000.0,0.0));
TopoDS_Edge botEdge = domBotBuilder.Edge();
BRepBuilderAPI_MakeWire domWireBuilder(inlEdge, topEdge, outEdge, botEdge);
TopoDS_Wire domWire = domWireBuilder.Wire();
BRepBuilderAPI_MakeFace face(domWire, true);
face.Add(aeroWire);
return face.Face();
}
}
}
......@@ -72,6 +72,8 @@ private:
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);
TopoDS_Shape BuildNACA(std::string naca);
/// OCC master object
TopoDS_Shape shape;
};
......
......@@ -57,6 +57,7 @@
#include <GProp_GProps.hxx>
#include <BRepGProp.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <GeomLProp_CLProps.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <GCPnts_AbscissaPoint.hxx>
#include <TopAbs_State.hxx>
......@@ -70,4 +71,12 @@
#include <ShapeAnalysis_Curve.hxx>
#include <Standard_Macro.hxx>
#include <GeomAPI_PointsToBSpline.hxx>
#include <Geom_BSplineCurve.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <STEPControl_Writer.hxx>
#endif
......@@ -51,6 +51,8 @@ ProcessLoadCAD::ProcessLoadCAD(MeshSharedPtr m) : ProcessModule(m)
{
m_config["filename"] =
ConfigOption(false, "", "Generate prisms on these surfs");
m_config["2D"] =
ConfigOption(true, "", "allow 2d loading");
}
ProcessLoadCAD::~ProcessLoadCAD()
......@@ -67,6 +69,12 @@ void ProcessLoadCAD::Process()
}
m_mesh->m_cad = GetEngineFactory().CreateInstance("oce",name);
if(m_config["2D"].beenSet)
{
m_mesh->m_cad->Set2D();
}
ASSERTL0(m_mesh->m_cad->LoadCAD(), "Failed to load CAD");
if (m_mesh->m_verbose)
......
......@@ -38,6 +38,7 @@ IF(NEKTAR_USE_MESHGEN)
CADSystem/OCE/CADVertOCE.cpp
CADSystem/OCE/CADCurveOCE.cpp
CADSystem/OCE/CADSurfOCE.cpp
2DGenerator/2DGenerator.cpp
)
ENDIF()
......@@ -92,6 +93,7 @@ IF(NEKTAR_USE_MESHGEN)
Optimisation/BGFS-B.h
Optimisation/OptimiseObj.h
Triangle/Triangle.h
2DGenerator/2DGenerator.h
)
ENDIF()
......
......@@ -35,10 +35,14 @@
#include "Octree.h"
#include <NekMeshUtils/CADSystem/CADSurf.h>
#include <NekMeshUtils/CADSystem/CADCurve.h>
#include <NekMeshUtils/Module/Module.h>
#include <LibUtilities/BasicUtils/ParseUtils.hpp>
#include <LibUtilities/BasicUtils/Progressbar.hpp>
#include <boost/algorithm/string.hpp>
using namespace std;
namespace Nektar
{
......@@ -868,7 +872,53 @@ struct linesource
void Octree::CompileSourcePointList()
{
//first sample surfaces
for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
{
CADCurveSharedPtr curve = m_mesh->m_cad->GetCurve(i);
Array<OneD, NekDouble> bds = curve->Bounds();
int samples = 100;
NekDouble dt = (bds[1] - bds[0]) / (samples + 1);
for (int j = 1; j < samples -1; j++) //dont want first and last point
{
NekDouble t = bds[0] + dt * j;
NekDouble C = curve->Curvature(t);
Array<OneD, NekDouble> loc = curve->P(t);
vector<CADSurfSharedPtr> ss = curve->GetAdjSurf();
Array<OneD, NekDouble> uv = ss[0]->locuv(loc);
if (C != 0.0)
{
NekDouble del = 2.0 * (1.0 / C) * sqrt(m_eps * (2.0 - m_eps));
if (del > m_maxDelta)
{
del = m_maxDelta;
}
if (del < m_minDelta)
{
del = m_minDelta;
}
CPointSharedPtr newCPoint =
MemoryManager<CPoint>::AllocateSharedPtr(ss[0]->GetId(), uv,
loc, del);
m_SPList.push_back(newCPoint);
}
else
{
BPointSharedPtr newBPoint =
MemoryManager<BPoint>::AllocateSharedPtr(ss[0]->GetId(), uv,
loc);
m_SPList.push_back(newBPoint);
}
}
}
for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
{
if(m_mesh->m_verbose)
......@@ -1000,35 +1050,34 @@ void Octree::CompileSourcePointList()
cout << endl;
}
if (m_udsfileset)
if (m_refinement.size() > 0)
{
if(m_mesh->m_verbose)
{
cout << "\t\tModifying based on uds files" << endl;
cout << "\t\tModifying based on refinement lines" << endl;
}
// now deal with the user defined spacing
vector<linesource> lsources;
fstream fle;
fle.open(m_udsfile.c_str());
vector<string> lines;
string fileline;
boost::split(lines, m_refinement, boost::is_any_of(":"));
while (!fle.eof())
for(int i = 0; i < lines.size(); i++)
{
getline(fle, fileline);
stringstream s(fileline);
string word;
s >> word;
if (word == "L")
{
Array<OneD, NekDouble> x1(3), x2(3);
NekDouble r, d;
s >> x1[0] >> x1[1] >> x1[2] >> x2[0] >> x2[1] >> x2[2]
>> r >> d;
lsources.push_back(linesource(x1, x2, r, d));
}
vector<NekDouble> data;
ParseUtils::GenerateUnOrderedVector(lines[i].c_str(), data);
Array<OneD, NekDouble> x1(3), x2(3);
x1[0] = data[0];
x1[1] = data[1];
x1[2] = data[2];
x2[0] = data[3];
x2[1] = data[4];
x2[2] = data[5];
lsources.push_back(linesource(x1, x2, data[6], data[7]));
}
fle.close();
// this takes any existing sourcepoints within the influence range
// and modifies them
......
......@@ -59,7 +59,6 @@ public:
Octree(MeshSharedPtr m) : m_mesh(m)
{
m_udsfileset = false;
}
/**
......@@ -113,10 +112,9 @@ public:
*
* @param nm name of the user defined spacing file
*/
void UDS(std::string nm)
void Refinement(std::string nm)
{
m_udsfile = nm;
m_udsfileset = true;
m_refinement = nm;
}
private:
......@@ -188,10 +186,8 @@ private:
int m_numoct;
/// Mesh object
MeshSharedPtr m_mesh;
/// user defined spacing has been set
bool m_udsfileset;
/// name of the user defined spacing file
std::string m_udsfile;
std::string m_refinement;
};
typedef boost::shared_ptr<Octree> OctreeSharedPtr;
......
......@@ -55,8 +55,8 @@ ProcessLoadOctree::ProcessLoadOctree(MeshSharedPtr m) : ProcessModule(m)
ConfigOption(false, "0", "mindelta.");
m_config["eps"] =
ConfigOption(false, "0", "mindelta.");
m_config["udsfile"] =
ConfigOption(false, "0", "mindelta.");
m_config["refinement"] =
ConfigOption(false, "", "mindelta.");
m_config["writeoctree"] =
ConfigOption(true, "0", "dump octree as xml mesh");
}
......@@ -87,9 +87,9 @@ void ProcessLoadOctree::Process()
m_mesh->m_octree->SetParameters(minDelta, maxDelta, eps);
if(m_config["udsfile"].beenSet)
if(m_config["refinement"].beenSet)
{
m_mesh->m_octree->UDS(m_config["udsfile"].as<string>());
m_mesh->m_octree->Refinement(m_config["refinement"].as<string>());
}
m_mesh->m_octree->Process();
......
......@@ -170,6 +170,7 @@ private:
int sid;
/// uv coord on surf
Array<OneD, NekDouble> m_uv;
NekDouble m_ti;
/// delta parameter
NekDouble m_delta;
};
......@@ -226,6 +227,7 @@ private:
int sid;
/// uv coord on surf
Array<OneD, NekDouble> m_uv;
NekDouble m_ti;
};
typedef boost::shared_ptr<BPoint> BPointSharedPtr;
......
......@@ -108,13 +108,13 @@ void CurveMesh::Mesh()
vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();
vector<CADSurfSharedPtr> s = m_cadcurve->GetAdjSurf();
ASSERTL0(s.size() == 2, "invalid curve");
ASSERTL0(s.size() == (m_mesh->m_cad->Is2D()) ? 1 : 2, "invalid curve");
NodeSharedPtr n = verts[0]->GetNode();
t = m_bounds[0];
n->SetCADCurve(m_id, m_cadcurve, t);
loc = n->GetLoc();
for (int j = 0; j < 2; j++)
for (int j = 0; j < s.size(); j++)
{
if (verts[0]->IsDegen() == s[j]->GetId()) // if the degen has been set
// for this node the node
......@@ -134,7 +134,7 @@ void CurveMesh::Mesh()
NodeSharedPtr n2 = boost::shared_ptr<Node>(
new Node(m_mesh->m_numNodes++, loc[0], loc[1], loc[2]));
n2->SetCADCurve(m_id, m_cadcurve, t);
for (int j = 0; j < 2; j++)
for (int j = 0; j < s.size(); j++)
{
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
n2->SetCADSurf(s[j]->GetId(), s[j], uv);
......@@ -146,7 +146,7 @@ void CurveMesh::Mesh()
t = m_bounds[1];
n->SetCADCurve(m_id, m_cadcurve, t);
loc = n->GetLoc();
for (int j = 0; j < 2; j++)
for (int j = 0; j < s.size(); j++)
{
if (verts[1]->IsDegen() == s[j]->GetId()) // if the degen has been set
// for this node the node
......
......@@ -1195,6 +1195,18 @@ void FaceMesh::OrientateCurves()
vector<NodeSharedPtr> tmp = orderedLoops[0];
reverse(tmp.begin(), tmp.end());
orderedLoops[0] = tmp;
//need to flip edgeo
for(int i = 0; i < m_edgeloops[0].edgeo.size(); i++)
{
if(m_edgeloops[0].edgeo[i] == 0)
{
m_edgeloops[0].edgeo[i] = 1;
}
else