Commit 8e19333c authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'feature/geo-reader' into 'master'

feature/geo-reader

See merge request !731
parents 7702eef3 0fd8d822
...@@ -3,7 +3,9 @@ Changelog ...@@ -3,7 +3,9 @@ Changelog
v4.5.0 v4.5.0
------ ------
**NekMesh**: **NekMesh**:
- Add feature to read basic 2D geo files as CAD (!731)
- Add periodic boundary condition meshing in 2D (!733) - Add periodic boundary condition meshing in 2D (!733)
- Adjust boundary layer thickness in corners in 2D (!739) - Adjust boundary layer thickness in corners in 2D (!739)
- Add non-O BL meshing in 2D (!757) - Add non-O BL meshing in 2D (!757)
...@@ -152,6 +154,7 @@ v4.4.0 ...@@ -152,6 +154,7 @@ v4.4.0
(!712) (!712)
- 2D to 3D mesh extrusion module (!715) - 2D to 3D mesh extrusion module (!715)
- Add new two-dimensional mesher from NACA code or step file (!720) - Add new two-dimensional mesher from NACA code or step file (!720)
- Add basic gmsh cad (.geo) reader to the meshing system (!731)
- Fix inverted boundary layer in 2D (!736) - Fix inverted boundary layer in 2D (!736)
- More sensible element sizing with boundary layers in 2D (!736) - More sensible element sizing with boundary layers in 2D (!736)
- Change variable names in mcf file to make more sense (!736) - Change variable names in mcf file to make more sense (!736)
......
...@@ -842,22 +842,22 @@ example which generates a 2D NACA wing. ...@@ -842,22 +842,22 @@ example which generates a 2D NACA wing.
<MESHING> <MESHING>
<INFORMATION> <INFORMATION>
<I PROPERTY="CADFile" VALUE="6412" /> <I PROPERTY="CADFile" VALUE="6412" />
<I PROPERTY="MeshType" VALUE="2D" /> <I PROPERTY="MeshType" VALUE="2D" />
</INFORMATION> </INFORMATION>
<PARAMETERS> <PARAMETERS>
<P PARAM="MinDelta" VALUE="0.01" /> <P PARAM="MinDelta" VALUE="0.01" />
<P PARAM="MaxDelta" VALUE="1.0" /> <P PARAM="MaxDelta" VALUE="1.0" />
<P PARAM="EPS" VALUE="0.1" /> <P PARAM="EPS" VALUE="0.1" />
<P PARAM="Order" VALUE="4" /> <P PARAM="Order" VALUE="4" />
<!-- 2D Domain !--> <!-- 2D Domain !-->
<P PARAM="Xmin" VALUE="-1.0" /> <P PARAM="Xmin" VALUE="-1.0" />
<P PARAM="Ymin" VALUE="-2.0" /> <P PARAM="Ymin" VALUE="-2.0" />
<P PARAM="Xmax" VALUE="3.0" /> <P PARAM="Xmax" VALUE="3.0" />
<P PARAM="Ymax" VALUE="2.0" /> <P PARAM="Ymax" VALUE="2.0" />
<P PARAM="AOA" VALUE="15.0" /> <P PARAM="AOA" VALUE="15.0" />
</PARAMETERS> </PARAMETERS>
</MESHING> </MESHING>
...@@ -867,10 +867,11 @@ In all cases the mesh generator needs two pieces of information and four ...@@ -867,10 +867,11 @@ In all cases the mesh generator needs two pieces of information and four
parameters. It firstly needs to know the CAD file with which to work. In the parameters. It firstly needs to know the CAD file with which to work. In the
example above this is listed as a 4 digit number, this is because the mesh example above this is listed as a 4 digit number, this is because the mesh
generator is equiped with a NACA wing generator. In all other cases this generator is equiped with a NACA wing generator. In all other cases this
parameter would be a STEP file. Secondly, what type of mesh to make, the options parameter would be the name of a CAD file (in either STEP or GEO format).
are \inltt{EULER} and \inltt{BL} for 3D meshes and \inltt{2D} and \inltt{2DBL} Secondly, what type of mesh to make, the options
are \inltt{EULER} and \inltt{BndLayer} for 3D meshes and \inltt{2D} and \inltt{2DBndLayer}
for 2D meshes. In the case of \inltt{EULER} the mesh will be made with only for 2D meshes. In the case of \inltt{EULER} the mesh will be made with only
tetrahedra. For \inltt{BL} the mesh generator will attempt to insert a single tetrahedra. For \inltt{BndLayer} the mesh generator will attempt to insert a single
macro prism layer onto the geometry surface. This option requires additional macro prism layer onto the geometry surface. This option requires additional
parameters. This is similar for the 2D scenarios. The automatic mesh parameters. This is similar for the 2D scenarios. The automatic mesh
specification system requires three parameters to build the specification of a specification system requires three parameters to build the specification of a
...@@ -889,35 +890,35 @@ An example is shown. ...@@ -889,35 +890,35 @@ An example is shown.
<MESHING> <MESHING>
<INFORMATION> <INFORMATION>
<I PROPERTY="CADFile" VALUE="6412" /> <I PROPERTY="CADFile" VALUE="6412" />
<I PROPERTY="MeshType" VALUE="2DBL" /> <I PROPERTY="MeshType" VALUE="2DBndLayer" />
</INFORMATION> </INFORMATION>
<PARAMETERS> <PARAMETERS>
<P PARAM="MinDelta" VALUE="0.01" /> <P PARAM="MinDelta" VALUE="0.01" />
<P PARAM="MaxDelta" VALUE="1.0" /> <P PARAM="MaxDelta" VALUE="1.0" />
<P PARAM="EPS" VALUE="0.1" /> <P PARAM="EPS" VALUE="0.1" />
<P PARAM="Order" VALUE="4" /> <P PARAM="Order" VALUE="4" />
<!-- Boundary layer !--> <!-- Boundary layer !-->
<P PARAM="BLSurfs" VALUE="5,6" /> <P PARAM="BndLayerSurfaces" VALUE="5,6" />
<P PARAM="BLThick" VALUE="0.03" /> <P PARAM="BndLayerThickness" VALUE="0.03" />
<P PARAM="BLLayers" VALUE="4" /> <P PARAM="BndLayerLayers" VALUE="4" />
<P PARAM="BLProg" VALUE="2.0" /> <P PARAM="BndLayerProgression" VALUE="2.0" />
<!-- 2D Domain !--> <!-- 2D Domain !-->
<P PARAM="Xmin" VALUE="-1.0" /> <P PARAM="Xmin" VALUE="-1.0" />
<P PARAM="Ymin" VALUE="-2.0" /> <P PARAM="Ymin" VALUE="-2.0" />
<P PARAM="Xmax" VALUE="3.0" /> <P PARAM="Xmax" VALUE="3.0" />
<P PARAM="Ymax" VALUE="2.0" /> <P PARAM="Ymax" VALUE="2.0" />
<P PARAM="AOA" VALUE="15.0" /> <P PARAM="AOA" VALUE="15.0" />
</PARAMETERS> </PARAMETERS>
</MESHING> </MESHING>
</NEKTAR> </NEKTAR>
\end{lstlisting} \end{lstlisting}
A list of the CAD surfaces which will have a prism generated on is described by A list of the CAD surfaces which will have a prism generated on is described by
\inltt{BLSurfs} and the thickness of the boundary to aim for is \inltt{BLThick}. \inltt{BndLayerSurfaces} and the thickness of the boundary to aim for is \inltt{BndLayerThickness}.
% %
The mesh generator has been created with a range of error messages to aid in The mesh generator has been created with a range of error messages to aid in
debugging. If you encounter an error and the mesh generator fails, run \nm with debugging. If you encounter an error and the mesh generator fails, run \nm with
...@@ -925,6 +926,32 @@ the verbose \inltt{-v} flag and send the stdout with the .mcf and .step files ...@@ -925,6 +926,32 @@ the verbose \inltt{-v} flag and send the stdout with the .mcf and .step files
to \inltt{m.turner14@imperial.ac.uk}. Without the feedback this functionality to \inltt{m.turner14@imperial.ac.uk}. Without the feedback this functionality
cannot improve. cannot improve.
\subsubsection{GEO format}
Recent developments have been made to facilitate the generation of meshes from
simple 2D geometries. The GEO file format, used by Gmsh, is a popular option
that allows the user to script geometrical and meshing operations without the
need of a GUI. A simplified reader has been implemented in NekMesh for 2D geometries.
Although very basic this reader may be extended in the future to cover a wider
range of geometrical features.
For a full description of the GEO format the user should refer to Gmsh's
documentation. The following commands are currently supported:
\begin{itemize}
\item \inltt{//} (comments)
\item \inltt{Point}
\item \inltt{Line}
\item \inltt{Spline}
\item \inltt{Line Loop}
\item \inltt{Plane Surface}
\end{itemize}
At the present time, NekMesh does not support the full scripting capabilities of the
GEO format. The used GEO files should be a straightforward succession of entity
creations (see list above). This should however allow for the creation of quite
a wide range of 2D geometries by transformation of arbitrary curves into generic
splines.
%%% Local Variables: %%% Local Variables:
%%% mode: latex %%% mode: latex
......
...@@ -64,7 +64,7 @@ namespace Nektar ...@@ -64,7 +64,7 @@ namespace Nektar
{ {
SymbolFunctor symbolFunctor(&symbol); SymbolFunctor symbolFunctor(&symbol);
ValueFunctor valueFunctor(&value); ValueFunctor valueFunctor(&value);
return parse(str, return parse(str,
// Begin grammar // Begin grammar
( (
...@@ -93,7 +93,7 @@ namespace Nektar ...@@ -93,7 +93,7 @@ namespace Nektar
space_p).full; space_p).full;
} }
static bool GenerateOrderedVector(const char *const str, std::vector<unsigned int> &vec) static bool GenerateOrderedVector(const char *const str, std::vector<unsigned int> &vec)
{ {
// Functors used to parse the sequence. // Functors used to parse the sequence.
...@@ -114,7 +114,7 @@ namespace Nektar ...@@ -114,7 +114,7 @@ namespace Nektar
{ {
// Functors used to parse the sequence. // Functors used to parse the sequence.
fctor4 functor4(&vec); fctor4 functor4(&vec);
return parse(str, return parse(str,
// Begin grammar // Begin grammar
( (
...@@ -124,12 +124,12 @@ namespace Nektar ...@@ -124,12 +124,12 @@ namespace Nektar
// End grammar // End grammar
space_p).full; space_p).full;
} }
static bool GenerateUnOrderedVector(const char *const str, std::vector<NekDouble> &vec) static bool GenerateUnOrderedVector(const char *const str, std::vector<NekDouble> &vec)
{ {
// Functors used to parse the sequence. // Functors used to parse the sequence.
fctor5 functor5(&vec); fctor5 functor5(&vec);
return parse(str, return parse(str,
// Begin grammar // Begin grammar
( (
...@@ -139,7 +139,22 @@ namespace Nektar ...@@ -139,7 +139,22 @@ namespace Nektar
// End grammar // End grammar
space_p).full; space_p).full;
} }
static bool GenerateUnOrderedVector(const char *const str, std::vector<unsigned int> &vec)
{
// Functors used to parse the sequence.
fctor6 functor6(&vec);
return parse(str,
// Begin grammar
(
uint_p[functor6] >> *(',' >> uint_p[functor6])
)
,
// End grammar
space_p).full;
}
static bool GenerateOrderedStringVector(const char *const str, std::vector<std::string> &vec) static bool GenerateOrderedStringVector(const char *const str, std::vector<std::string> &vec)
{ {
// Functors used to parse the sequence. // Functors used to parse the sequence.
...@@ -226,7 +241,7 @@ namespace Nektar ...@@ -226,7 +241,7 @@ namespace Nektar
m_value(value) m_value(value)
{ {
} }
void operator()(NekDouble val) const void operator()(NekDouble val) const
{ {
*m_value = val; *m_value = val;
...@@ -245,7 +260,7 @@ namespace Nektar ...@@ -245,7 +260,7 @@ namespace Nektar
void operator()(unsigned int n) const void operator()(unsigned int n) const
{ {
#ifdef NOTREQUIRED //SJS: I do not think we need this check #ifdef NOTREQUIRED //SJS: I do not think we need this check
if (!m_vector->empty()) if (!m_vector->empty())
{ {
unsigned int prevElem = m_vector->back(); unsigned int prevElem = m_vector->back();
...@@ -304,7 +319,7 @@ namespace Nektar ...@@ -304,7 +319,7 @@ namespace Nektar
std::vector<std::string> *m_vector; std::vector<std::string> *m_vector;
}; };
// Probably should template fctor1 if that is possible? // Probably should template fctor1 if that is possible?
struct fctor4 struct fctor4
{ {
fctor4(std::vector<NekDouble> *vec): fctor4(std::vector<NekDouble> *vec):
...@@ -333,24 +348,41 @@ namespace Nektar ...@@ -333,24 +348,41 @@ namespace Nektar
std::vector<NekDouble> *m_vector; std::vector<NekDouble> *m_vector;
fctor4(); fctor4();
}; };
struct fctor5 struct fctor5
{ {
fctor5(std::vector<NekDouble> *vec): fctor5(std::vector<NekDouble> *vec):
m_vector(vec) m_vector(vec)
{ {
} }
void operator()(NekDouble n) const void operator()(NekDouble n) const
{ {
m_vector->push_back(n); m_vector->push_back(n);
} }
private: private:
std::vector<NekDouble> *m_vector; std::vector<NekDouble> *m_vector;
fctor5(); fctor5();
}; };
struct fctor6
{
fctor6(std::vector<unsigned int> *vec):
m_vector(vec)
{
}
void operator()(unsigned int n) const
{
m_vector->push_back(n);
}
private:
std::vector<unsigned int> *m_vector;
fctor6();
};
}; };
} }
......
...@@ -40,6 +40,9 @@ ...@@ -40,6 +40,9 @@
#include <NekMeshUtils/CADSystem/OCE/CADSystemOCE.h> #include <NekMeshUtils/CADSystem/OCE/CADSystemOCE.h>
#include <NekMeshUtils/CADSystem/OCE/CADVertOCE.h> #include <NekMeshUtils/CADSystem/OCE/CADVertOCE.h>
#include <boost/algorithm/string.hpp>
#include <boost/filesystem.hpp>
using namespace std; using namespace std;
namespace Nektar namespace Nektar
...@@ -54,17 +57,27 @@ bool CADSystemOCE::LoadCAD() ...@@ -54,17 +57,27 @@ bool CADSystemOCE::LoadCAD()
{ {
if (m_naca.size() == 0) if (m_naca.size() == 0)
{ {
// not a naca profile behave normally //not a naca profile behave normally
// Takes step file and makes OpenCascade shape //could be a geo
STEPControl_Reader reader; string ext = boost::filesystem::extension(m_name);
reader = STEPControl_Reader();
reader.ReadFile(m_name.c_str()); if (boost::iequals(ext, ".geo"))
reader.NbRootsForTransfer(); {
reader.TransferRoots(); shape = BuildGeo(m_name);
shape = reader.OneShape(); }
if (shape.IsNull()) else
{ {
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 else
...@@ -388,8 +401,166 @@ TopoDS_Shape CADSystemOCE::BuildNACA(string naca) ...@@ -388,8 +401,166 @@ TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
ShapeFix_Face sf(face.Face()); ShapeFix_Face sf(face.Face());
sf.FixOrientation(); sf.FixOrientation();
return sf.Face();
}
TopoDS_Shape CADSystemOCE::BuildGeo(string geo)
{
ifstream f;
f.open(geo.c_str());
map<int, string> points;
map<int, string> lines;
map<int, string> splines;
map<int, string> loops;
map<int, string> surfs;
string fline;
string flinetmp;
while (!f.eof())
{
getline(f, fline);
if (fline.size() == 0)
{
continue;
}
if (boost::starts_with(fline, "//"))
{
continue;
}
if (!boost::contains(fline, ";"))
{
flinetmp += fline;
continue;
}
fline = flinetmp + fline;
flinetmp.clear();
vector<string> tmp1, tmp2;
boost::split(tmp1, fline, boost::is_any_of("="));
boost::split(tmp2, tmp1[0], boost::is_any_of("("));
string type = tmp2[0];
boost::erase_all(tmp2[1], ")");
boost::erase_all(tmp2[1], " ");
int id = boost::lexical_cast<int>(tmp2[1]);
boost::erase_all(tmp1[1], " ");
boost::erase_all(tmp1[1], "{");
boost::erase_all(tmp1[1], "}");
boost::erase_all(tmp1[1], ";");
string var = tmp1[1];
if (boost::iequals(type, "Point"))
{
points[id] = var;
}
else if (boost::iequals(type, "Line"))
{
lines[id] = var;
}
else if (boost::iequals(type, "Spline"))
{
splines[id] = var;
}
else if (boost::iequals(type, "Line Loop"))
{
//line loops sometimes have negative entries for gmsh
//orientaton purposes
//we dont care so remove it
boost::erase_all(var, "-");
loops[id] = var;
}
else if (boost::iequals(type, "Plane Surface"))
{
surfs[id] = var;
}
else
{
cout << "not sure " << type << endl;
}
}
map<int, string>::iterator it;
// build points
map<int, gp_Pnt> cPoints;
for (it = points.begin(); it != points.end(); it++)
{
vector<NekDouble> data;
ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
cPoints[it->first] =
gp_Pnt(data[0] * 1000.0, data[1] * 1000.0, data[2] * 1000.0);
}
// build edges
map<int, TopoDS_Edge> cEdges;
for (it = lines.begin(); it != lines.end(); it++)
{
vector<unsigned int> data;
ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
BRepBuilderAPI_MakeEdge em(cPoints[data[0]], cPoints[data[1]]);
cEdges[it->first] = em.Edge();
}
for (it = splines.begin(); it != splines.end(); it++)
{
vector<unsigned int> data;
ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
TColgp_Array1OfPnt pointArray(0, data.size() - 1);
for (int i = 0; i < data.size(); i++)
{
pointArray.SetValue(i, cPoints[data[i]]);
}
GeomAPI_PointsToBSpline spline(pointArray);
Handle(Geom_BSplineCurve) curve = spline.Curve();
BRepBuilderAPI_MakeEdge em(curve);
cEdges[it->first] = em.Edge();
}
// build wires
map<int, TopoDS_Wire> cWires;
for (it = loops.begin(); it != loops.end(); it++)
{
vector<unsigned int> data;
ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
BRepBuilderAPI_MakeWire wm;
for (int i = 0; i < data.size(); i++)
{
wm.Add(cEdges[data[i]]);
}
cWires[it->first] = wm.Wire();
}
// make surface, at this point assuming its 2D (therefore only 1)
// also going to assume that the first loop in the list is the outer domain
ASSERTL0(surfs.size() == 1, "more than 1 surf");
it = surfs.begin();
vector<unsigned int> data;
ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
BRepBuilderAPI_MakeFace face(cWires[data[0]], true);
for (int i = 1; i < data.size(); i++)
{
face.Add(cWires[data[i]]);
}
ASSERTL0(face.Error() == BRepBuilderAPI_FaceDone, "build geo failed");
ShapeFix_Face sf(face.Face());
sf.FixOrientation();
return sf.Face(); return sf.Face();
} }
} }
} }
...@@ -76,6 +76,7 @@ private: ...@@ -76,6 +76,7 @@ private: