Commit 4189f8fc by Michael Turner

Merge branch 'feature/CAD-update' into 'master'

feature/CAD-update

See merge request !822
parents 7f669950 2b00e54c
......@@ -34,6 +34,7 @@ v5.0.0
- Fix issue with reading CCM files due to definition of default arrays
rather than a vector (!797)
- Fix inverted triangles and small memory issue in surface meshing (!798)
- Update for the CAD system, more advance self-healing and analysis (!822)
- Additional curve types in GEO reader: BSpline, Circle, Ellipse (!800)
**FieldConvert**:
......
......@@ -24,19 +24,28 @@ if(NOT DEFINED OCE_DIR)
# Check for OSX needs to come first because UNIX evaluates to true on OSX
if(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
if(DEFINED MACPORTS_PREFIX)
find_package(OCE 0.15 QUIET HINTS ${MACPORTS_PREFIX}/Library/Frameworks)
find_package(OCE 0.17 QUIET HINTS ${MACPORTS_PREFIX}/Library/Frameworks)
elseif(DEFINED HOMEBREW_PREFIX)
find_package(OCE 0.15 QUIET HINTS ${HOMEBREW_PREFIX}/Cellar/oce/*)
find_package(OCE 0.17 QUIET HINTS ${HOMEBREW_PREFIX}/Cellar/oce/*)
endif()
elseif(UNIX)
set(OCE_DIR "/usr/local/share/cmake/")
endif()
endif()
find_package(OCE 0.15 QUIET)
find_package(OCE 0.17 QUIET)
if(OCE_FOUND)
message(STATUS "-- OpenCASCADE Community Edition has been found.")
set(OCC_INCLUDE_DIR ${OCE_INCLUDE_DIRS})
message(STATUS "OpenCASCADE Community Edition has been found.")
#check that the OCE package has CAF modules
FIND_LIBRARY(OCC_CAF_LIBRARY TKXCAF ${OCE_INCLUDE_DIRS}/../../lib )
if(OCC_CAF_LIBRARY)
set(OCC_INCLUDE_DIR ${OCE_INCLUDE_DIRS})
else()
message(STATUS "-- OCE does not have CAF libraries, will build from source.")
endif()
else(OCE_FOUND) #look for OpenCASCADE
FIND_PATH(OCC_INCLUDE_DIR Standard_Version.hxx
......@@ -45,7 +54,6 @@ else(OCE_FOUND) #look for OpenCASCADE
/usr/local/opt/opencascade/include
/opt/opencascade/include
/opt/opencascade/inc
/opt/local/include/oce
)
FIND_LIBRARY(OCC_LIBRARY TKernel
/usr/lib
......@@ -56,6 +64,7 @@ else(OCE_FOUND) #look for OpenCASCADE
)
if(OCC_LIBRARY)
message(STATUS "OpenCASCADE has been found.")
GET_FILENAME_COMPONENT(OCC_LIBRARY_DIR ${OCC_LIBRARY} PATH)
IF(NOT OCC_INCLUDE_DIR)
FIND_PATH(OCC_INCLUDE_DIR Standard_Version.hxx
......@@ -112,9 +121,11 @@ if(OCC_FOUND)
TKSTEPAttr
TKHLR
TKFeat
TKXCAF
TKXDESTEP
)
if(OCC_VERSION_STRING VERSION_LESS 6.7)
if(OCC_VERSION_STRING VERSION_LESS 6.8)
MESSAGE(SEND_ERROR "OCC version too low")
endif(OCC_VERSION_STRING VERSION_LESS 6.7)
endif(OCC_VERSION_STRING VERSION_LESS 6.8)
message(STATUS "-- Found OCE/OpenCASCADE with OCC version: ${OCC_VERSION_STRING}")
endif(OCC_FOUND)
......@@ -22,9 +22,34 @@ IF(NEKTAR_USE_MESHGEN)
IF (THIRDPARTY_BUILD_OCE)
INCLUDE(ExternalProject)
SET(OCC_LIBRARIES_TMP PTKernel TKernel TKMath TKBRep TKIGES TKSTEP TKSTEPAttr
TKSTEP209 TKSTEPBase TKShapeSchema TKGeomBase TKGeomAlgo TKG3d TKG2d
TKXSBase TKPShape TKTopAlgo TKShHealing)
SET(OCC_LIBRARIES_TMP
TKFillet
TKMesh
TKernel
TKG2d
TKG3d
TKMath
TKIGES
TKSTL
TKShHealing
TKXSBase
TKBool
TKBO
TKBRep
TKTopAlgo
TKGeomAlgo
TKGeomBase
TKOffset
TKPrim
TKSTEP
TKSTEPBase
TKSTEPAttr
TKHLR
TKFeat
TKXCAF
TKXDESTEP
)
FOREACH(OCC_LIB ${OCC_LIBRARIES_TMP})
LIST(APPEND OCC_LIBRARIES ${TPDIST}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}${OCC_LIB}${CMAKE_SHARED_LIBRARY_SUFFIX})
ENDFOREACH()
......@@ -51,8 +76,6 @@ IF(NEKTAR_USE_MESHGEN)
-DOCE_INSTALL_PREFIX:PATH=${TPDIST}
-DOCE_TESTING=OFF
-DOCE_VISUALISATION=OFF
-DOCE_DISABLE_X11=ON
-DOCE_OCAF=OFF
${TPSRC}/oce-0.17
)
......
......@@ -162,6 +162,8 @@ public:
* @brief locates a point in the parametric space
*/
virtual NekDouble loct(Array<OneD, NekDouble> xyz) = 0;
virtual NekDouble DistanceTo(Array<OneD, NekDouble> xyz) = 0;
CADOrientation::Orientation GetOrienationWRT(int surf)
{
......
......@@ -32,15 +32,14 @@
// Description: cad object methods.
//
////////////////////////////////////////////////////////////////////////////////
#include "CADSurf.h"
#include "CADCurve.h"
#include <boost/geometry.hpp>
#include <boost/geometry/algorithms/assign.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include "CADSurf.h"
#include "CADCurve.h"
namespace bg = boost::geometry;
typedef bg::model::d2::point_xy<double> point_xy;
......
......@@ -145,13 +145,17 @@ public:
* @brief takes a point from anywhere find the nearest surface point and its
* uv
*/
virtual void ProjectTo(Array<OneD, NekDouble> &tp,
virtual NekDouble ProjectTo(Array<OneD, NekDouble> &tp,
Array<OneD, NekDouble> &uv) = 0;
virtual Array<OneD, NekDouble> BoundingBox() = 0;
/**
* @brief returns curvature at point uv
*/
virtual NekDouble Curvature(Array<OneD, NekDouble> uv) = 0;
virtual bool IsPlanar() = 0;
/**
* @brief query reversed normal
......@@ -160,6 +164,16 @@ public:
{
return m_orientation;
}
void SetName(std::string i)
{
m_name = i;
}
std::string GetName()
{
return m_name;
}
protected:
/// List of bounding edges in loops with orientation.
......@@ -167,6 +181,8 @@ protected:
/// Function which tests the the value of uv used is within the surface
virtual void Test(Array<OneD, NekDouble> uv) = 0;
std::string m_name;
};
typedef std::shared_ptr<CADSurf> CADSurfSharedPtr;
......
......@@ -101,5 +101,12 @@ Array<OneD, NekDouble> CADSystem::GetPeriodicTranslationVector(int first,
return ret;
}
std::string CADSystem::GetSurfaceName(int i)
{
return m_surfs[i]->GetName();
}
}
}
......@@ -93,7 +93,7 @@ public:
m_2d = false;
}
~CADSystem()
virtual ~CADSystem()
{
}
......@@ -168,8 +168,8 @@ public:
*/
CADCurveSharedPtr GetCurve(int i)
{
std::map<int, CADCurveSharedPtr>::iterator search = m_curves.find(i);
ASSERTL0(search != m_curves.end(), "curve does not exist");
auto search = m_curves.find(i);
ASSERTL1(search != m_curves.end(), "curve does not exist");
return search->second;
}
......@@ -179,8 +179,16 @@ public:
*/
CADSurfSharedPtr GetSurf(int i)
{
std::map<int, CADSurfSharedPtr>::iterator search = m_surfs.find(i);
ASSERTL0(search != m_surfs.end(), "surface does not exist");
auto search = m_surfs.find(i);
ASSERTL1(search != m_surfs.end(), "surface does not exist");
return search->second;
}
CADVertSharedPtr GetVert(int i)
{
auto search = m_verts.find(i);
ASSERTL1(search != m_verts.end(), "vert does not exist");
return search->second;
}
......@@ -200,7 +208,9 @@ public:
{
return m_verts.size();
}
std::string GetSurfaceName(int i);
NEKMESHUTILS_EXPORT Array<OneD, NekDouble> GetPeriodicTranslationVector(
int first, int second);
......
......@@ -64,9 +64,7 @@ public:
m_type = CADType::eVert;
}
virtual ~CADVert()
{
}
virtual ~CADVert(){};
/**
* @brief Get x,y,z location of the vertex
......@@ -115,6 +113,8 @@ public:
return -1;
}
}
virtual NekDouble DistanceTo(Array<OneD, NekDouble> l) = 0;
protected:
/// mesh convert object of vert
......
......@@ -81,21 +81,31 @@ NekDouble CADCurveOCE::Length(NekDouble ti, NekDouble tf)
NekDouble CADCurveOCE::loct(Array<OneD, NekDouble> xyz)
{
NekDouble t = 0.0;
Array<OneD, NekDouble> b = GetBounds();
gp_Pnt loc(xyz[0]*1000.0, xyz[1]*1000.0, xyz[2]*1000.0);
ShapeAnalysis_Curve sac;
gp_Pnt p;
sac.Project(m_c,loc,1e-8 ,p,t);
WARNINGL2(p.Distance(loc) < 1e-3, "large locuv distance " +
boost::lexical_cast<string>(p.Distance(loc)/1000.0));
if(p.Distance(loc) > 1e-5)
{
cerr << "large loct distance" << endl;
}
return t;
}
NekDouble CADCurveOCE::DistanceTo(Array<OneD, NekDouble> xyz)
{
NekDouble t = 0.0;
gp_Pnt loc(xyz[0]*1000.0, xyz[1]*1000.0, xyz[2]*1000.0);
ShapeAnalysis_Curve sac;
gp_Pnt p;
sac.Project(m_c,loc,1e-8 ,p,t);
return p.Distance(loc)/1000.0;
}
Array<OneD, NekDouble> CADCurveOCE::P(NekDouble t)
{
Array<OneD, NekDouble> location(3);
......
......@@ -59,7 +59,7 @@ public:
{
}
virtual ~CADCurveOCE()
~CADCurveOCE()
{
}
......@@ -73,6 +73,7 @@ public:
virtual NekDouble Curvature(NekDouble t);
virtual Array<OneD, NekDouble> NormalWRT(NekDouble t, int surf);
virtual Array<OneD, NekDouble> N(NekDouble t);
virtual NekDouble DistanceTo(Array<OneD, NekDouble> xyz);
void Initialise(int i, TopoDS_Shape in)
{
......
......@@ -71,6 +71,10 @@ void CADSurfOCE::Initialise(int i, TopoDS_Shape in)
m_bounds[3]);
m_sas = new ShapeAnalysis_Surface(m_s);
m_sas->SetDomain(m_bounds[0], m_bounds[1], m_bounds[2], m_bounds[3]);
m_shape = in;
m_2Dclass = new BRepTopAdaptor_FClass2d(TopoDS::Face(m_shape), 1e-4);
}
Array<OneD, NekDouble> CADSurfOCE::GetBounds()
......@@ -78,6 +82,35 @@ Array<OneD, NekDouble> CADSurfOCE::GetBounds()
return m_bounds;
}
bool CADSurfOCE::IsPlanar()
{
if(m_occSurface.GetType() == GeomAbs_Plane)
{
return true;
}
return false;
}
Array<OneD, NekDouble> CADSurfOCE::BoundingBox()
{
BRepMesh_IncrementalMesh brmsh;
brmsh.SetShape(m_shape);
brmsh.SetDeflection(0.005);
brmsh.Perform();
Bnd_Box B;
BRepBndLib::Add(m_shape, B);
NekDouble e = sqrt(B.SquareExtent()) * 0.01;
e = min(e, 5e-3);
B.Enlarge(e);
Array<OneD, NekDouble> ret(6);
B.Get(ret[0],ret[1],ret[2],ret[3],ret[4],ret[5]);
return ret;
}
Array<OneD, NekDouble> CADSurfOCE::locuv(Array<OneD, NekDouble> p)
{
// has to transfer back to mm
......@@ -89,9 +122,9 @@ Array<OneD, NekDouble> CADSurfOCE::locuv(Array<OneD, NekDouble> p)
uvr[0] = p2.X();
uvr[1] = p2.Y();
WARNINGL2(m_sas->Value(p2).Distance(loc) < 1e-3, "large locuv distance " +
boost::lexical_cast<string>(
m_sas->Value(p2).Distance(loc)/1000.0) + " " +
gp_Pnt p3 = m_sas->Value(p2);
WARNINGL2(p3.Distance(loc) < 1e-3, "large locuv distance " +
boost::lexical_cast<string>(p3.Distance(loc)/1000.0) + " " +
boost::lexical_cast<string>(m_id));
// if the uv returned is slightly off the surface
......@@ -172,18 +205,47 @@ NekDouble CADSurfOCE::DistanceTo(Array<OneD, NekDouble> p)
gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
gp_Pnt2d p2 = m_sas->ValueOfUV(loc, Precision::Confusion());
TopAbs_State s = m_2Dclass->Perform(p2);
if(s == TopAbs_OUT)
{
BRepBuilderAPI_MakeVertex v(gp_Pnt(p[0],p[1],p[2]));
BRepExtrema_DistShapeShape dss(BRepTools::OuterWire(TopoDS::Face(m_shape)), v.Shape());
dss.Perform();
return dss.Value();
//return numeric_limits<double>::max();
}
gp_Pnt p3 = m_sas->Value(p2);
return p3.Distance(loc);
return p3.Distance(loc) / 1000.0;
}
void CADSurfOCE::ProjectTo(Array<OneD, NekDouble> &tp,
NekDouble CADSurfOCE::ProjectTo(Array<OneD, NekDouble> &tp,
Array<OneD, NekDouble> &uv)
{
gp_Pnt loc(tp[0] * 1000.0, tp[1] * 1000.0, tp[2] * 1000.0);
gp_Pnt2d p2 = m_sas->ValueOfUV(loc, Precision::Confusion());
TopAbs_State s = m_2Dclass->Perform(p2);
if(s == TopAbs_OUT)
{
BRepBuilderAPI_MakeVertex v(gp_Pnt(tp[0],tp[1],tp[2]));
BRepExtrema_DistShapeShape dss(BRepTools::OuterWire(TopoDS::Face(m_shape)), v.Shape());
dss.Perform();
gp_Pnt np = dss.PointOnShape1(1);
tp[0] = np.X();
tp[1] = np.Y();
tp[2] = np.Z();
gp_Pnt loc2(tp[0] * 1000.0, tp[1] * 1000.0, tp[2] * 1000.0);
gp_Pnt2d p22 = m_sas->ValueOfUV(loc2, Precision::Confusion());
uv[0] = p22.X();
uv[1] = p22.Y();
return dss.Value();
}
gp_Pnt p3 = m_sas->Value(p2);
......@@ -193,6 +255,8 @@ void CADSurfOCE::ProjectTo(Array<OneD, NekDouble> &tp,
uv[0] = p2.X();
uv[1] = p2.Y();
return p3.Distance(loc) / 1000.0;
}
Array<OneD, NekDouble> CADSurfOCE::P(Array<OneD, NekDouble> uv)
......
......@@ -59,7 +59,7 @@ public:
{
}
virtual ~CADSurfOCE()
~CADSurfOCE()
{
}
......@@ -72,9 +72,11 @@ public:
virtual Array<OneD, NekDouble> P (Array<OneD, NekDouble> uv);
virtual Array<OneD, NekDouble> locuv(Array<OneD, NekDouble> p);
virtual NekDouble DistanceTo(Array<OneD, NekDouble> p);
virtual void ProjectTo(Array<OneD, NekDouble> &tp,
virtual NekDouble ProjectTo(Array<OneD, NekDouble> &tp,
Array<OneD, NekDouble> &uv);
virtual NekDouble Curvature(Array<OneD, NekDouble> uv);
virtual Array<OneD, NekDouble> BoundingBox();
virtual bool IsPlanar();
private:
/// Function which tests the the value of uv used is within the surface
......@@ -87,6 +89,10 @@ private:
Array<OneD, NekDouble> m_bounds;
/// locuv object (stored because it gets faster with stored information)
ShapeAnalysis_Surface *m_sas;
/// original shape
TopoDS_Shape m_shape;
///
BRepTopAdaptor_FClass2d *m_2Dclass;
};
typedef std::shared_ptr<CADSurf> CADSurfSharedPtr;
......
......@@ -54,8 +54,30 @@ namespace NekMeshUtils
std::string CADSystemOCE::key = GetEngineFactory().RegisterCreatorFunction(
"oce", CADSystemOCE::create, "Uses OCE as cad engine");
void filterModShape(TopTools_DataMapOfShapeShape &modShape, TopoDS_Shape &S)
{
bool repeat = true;
while (repeat)
{
repeat = false;
if (modShape.IsBound(S))
{
repeat = true;
S = modShape.Find(S);
}
}
}
bool CADSystemOCE::LoadCAD()
{
Handle(XSControl_WorkSession) WS;
Handle(Interface_InterfaceModel) Model;
Handle(XSControl_TransferReader) TR;
Handle(Transfer_TransientProcess) TP;
Handle(XCAFDoc_ShapeTool) STool;
bool fromStep = false;
if (m_naca.size() == 0)
{
// not a naca profile behave normally
......@@ -69,9 +91,24 @@ bool CADSystemOCE::LoadCAD()
else
{
// Takes step file and makes OpenCascade shape
STEPControl_Reader reader;
reader = STEPControl_Reader();
reader.ReadFile(m_name.c_str());
STEPCAFControl_Reader readerCAF;
readerCAF.SetNameMode(true);
readerCAF.SetLayerMode(true);
readerCAF.SetColorMode(true);
Handle(TDocStd_Document) document =
new TDocStd_Document(Storage::Version());
readerCAF.ReadFile(m_name.c_str());
readerCAF.Transfer(document);
STEPControl_Reader reader = readerCAF.Reader();
WS = reader.WS();
Model = WS->Model();
TR = WS->TransferReader();
TP = TR->TransientProcess();
XCAFDoc_DocumentTool::ShapeTool(document->Main());
reader.NbRootsForTransfer();
reader.TransferRoots();
shape = reader.OneShape();
......@@ -79,6 +116,7 @@ bool CADSystemOCE::LoadCAD()
{
return false;
}
fromStep = true;
}
}
else
......@@ -88,6 +126,39 @@ bool CADSystemOCE::LoadCAD()
TopExp_Explorer explr;
TopTools_DataMapOfShapeShape modShape;
if(!m_2d)
{
BRepBuilderAPI_Sewing sew(1e-1);
for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
{
sew.Add(explr.Current());
}
sew.Perform();
for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
{
if (sew.IsModified(explr.Current()))
{
modShape.Bind(explr.Current(), sew.Modified(explr.Current()));
}
}
shape = sew.SewedShape();
int shell = 0;
for (explr.Init(shape, TopAbs_SHELL); explr.More(); explr.Next())
{
shell++;
}
ASSERTL0(shell == 1,
"Was not able to form a topological water tight shell");
}
// build map of verticies
for (explr.Init(shape, TopAbs_VERTEX); explr.More(); explr.Next())
{
......@@ -110,8 +181,8 @@ bool CADSystemOCE::LoadCAD()
{
continue;
}
BRepAdaptor_Curve curve = BRepAdaptor_Curve(TopoDS::Edge(e));
if (curve.GetType() != 7)
if (!BRep_Tool::Degenerated(TopoDS::Edge(e)))
{
int i = mapOfEdges.Add(e);
AddCurve(i, e);
......@@ -127,6 +198,59 @@ bool CADSystemOCE::LoadCAD()
AddSurf(i, f);
}
// attemps to extract patch names from STEP file
if(fromStep)
{
int nb = Model->NbEntities();
for (int i = 1; i <= nb; i++)
{
if (!Model->Value(i)->DynamicType()->SubType(
"StepRepr_RepresentationItem"))
continue;
Handle(StepRepr_RepresentationItem) enti =
Handle(StepRepr_RepresentationItem)::DownCast(Model->Value(i));
Handle(TCollection_HAsciiString) name = enti->Name();
if (name->IsEmpty())
continue;
Handle(Transfer_Binder) binder = TP->Find(Model->Value(i));
if (binder.IsNull() || !binder->HasResult())
continue;
TopoDS_Shape S = TransferBRep::ShapeResult(TP, binder);
if (S.IsNull())
continue;
if (S.ShapeType() == TopAbs_FACE)
{
string s(name->ToCString());
if (mapOfFaces.Contains(S))
{
int id = mapOfFaces.FindIndex(S);
m_surfs[id]->SetName(s);
}
else
{
filterModShape(modShape, S);
if (mapOfFaces.Contains(S))
{
int id = mapOfFaces.FindIndex(S);
m_surfs[id]->SetName(s);
}
else
{
ASSERTL0(false, "Name error");
}
}
}
}
}
// attempts to identify properties of the vertex on the degen edge
for (int i = 1; i <= mapOfFaces.Extent(); i++)
{
......@@ -156,11 +280,10 @@ bool CADSystemOCE::LoadCAD()
// This checks that all edges are bound by two surfaces, sanity check.
if (!m_2d)
{
map<int, CADCurveSharedPtr>::iterator it;
for (it = m_curves.begin(); it != m_curves.end(); it++)
for (auto &i : m_curves)
{
ASSERTL0(it->second->GetAdjSurf().size() == 2,
"curve is not joined to 2 surfaces");
ASSERTL0(i.second->GetAdjSurf().size() == 2,
"topolgy error found, surface not closed");
}
}
......
......@@ -59,13 +59,16 @@ public:
* @brief Default constructor.
*/
CADSystemOCE(std::string name) : CADSystem(name) {}
virtual ~CADSystemOCE()
{
}
~CADSystemOCE(){};
bool LoadCAD();
Array<OneD, NekDouble> GetBoundingBox();
TopoDS_Shape GetShape()
{
return shape;
}
private:
/// Function to add curve to CADSystem::m_verts.
......@@ -82,6 +85,8 @@ private:
TopTools_IndexedMapOfShape mapOfVerts, mapOfEdges, mapOfFaces;
};
typedef boost::shared_ptr<CADSystemOCE> CADSystemOCESharedPtr;
}
}
......
......@@ -62,11 +62,17 @@ public:
{
}
virtual ~CADVertOCE()
~CADVertOCE()
{
}
void Initialise(int i, TopoDS_Shape in);
NekDouble DistanceTo(Array<OneD, NekDouble> l)
{
gp_Pnt lp(l[0],l[1],l[2]);
return m_occVert.Distance(lp);
}
private:
/// OpenCascade object of the curve.
......
......@@ -38,49 +38,73 @@
/// This is a list of OpenCascade headers required for use with nektar
#include <STEPControl_Reader.hxx>
#include <TColStd_HSequenceOfTransient.hxx>
#include <TopoDS.hxx>
/// IO classes
#include <STEPCAFControl_Reader.hxx>
#include <StepRepr_RepresentationItem.hxx>
#include <TDocStd_Document.hxx>
#include <XSControl_WorkSession.hxx>
#include <XSControl_TransferReader.hxx>
#include <XCAFDoc_DocumentTool.hxx>
#include <Storage.hxx>
/// STL classes
#include <BRepMesh_IncrementalMesh.hxx>
/// Shape Analysis / exploration classes
#include <BRepTopAdaptor_FClass2d.hxx>
#include <ShapeAnalysis_Surface.hxx>
#include <BRepTools_WireExplorer.hxx>