Commit e17f2d81 authored by cadfix's avatar cadfix

working version

parent 6af249c0
......@@ -75,7 +75,7 @@ void Generator2D::Process()
// Check that cad is 2D
Array<OneD, NekDouble> bndBox = m_mesh->m_cad->GetBoundingBox();
ASSERTL0(fabs(bndBox[5] - bndBox[4]) < 1.0e-7, "CAD isn't 2D");
if (m_mesh->m_verbose)
{
cout << endl << "2D meshing" << endl;
......@@ -385,15 +385,14 @@ void Generator2D::MakeBL(int faceid)
// constructed by computing a normal but found on the adjacent curve
if (it->second.size() == 1)
{
vector<pair<int, CADCurveSharedPtr> > curves =
it->first->GetCADCurves();
vector<CADCurveSharedPtr> curves = it->first->GetCADCurves();
vector<EdgeSharedPtr> edges =
m_curvemeshes[curves[0].first]->GetMeshEdges();
m_curvemeshes[curves[0]->GetId()]->GetMeshEdges();
vector<EdgeSharedPtr>::iterator ie =
find(edges.begin(), edges.end(), it->second[0]);
int rightCurve =
(ie == edges.end()) ? curves[0].first : curves[1].first;
(ie == edges.end()) ? curves[0]->GetId() : curves[1]->GetId();
vector<NodeSharedPtr> nodes =
m_curvemeshes[rightCurve]->GetMeshPoints();
......@@ -429,7 +428,8 @@ void Generator2D::MakeBL(int faceid)
NodeSharedPtr nn = std::shared_ptr<Node>(
new Node(m_mesh->m_numNodes++, n[0], n[1], 0.0));
CADSurfSharedPtr s = m_mesh->m_cad->GetSurf(faceid);
Array<OneD, NekDouble> uv = s->locuv(n);
Array<OneD, NekDouble> uv;
s->locuv(n, uv);
nn->SetCADSurf(s, uv);
nodeNormals[it->first] = nn;
}
......
......@@ -160,11 +160,11 @@ public:
}
/*
* @brief locates a point in the parametric space
* @brief locates a point in the parametric space. returns the
* distance to the point and passes t by reference and updates it
*/
virtual NekDouble loct(Array<OneD, NekDouble> xyz) = 0;
virtual NekDouble DistanceTo(Array<OneD, NekDouble> xyz) = 0;
virtual NekDouble loct(Array<OneD, NekDouble> xyz,
NekDouble &t) = 0;
CADOrientation::Orientation GetOrienationWRT(int surf)
{
......@@ -182,7 +182,6 @@ public:
virtual Array<OneD, NekDouble> NormalWRT(NekDouble t, int surf) = 0;
virtual Array<OneD, NekDouble> N(NekDouble t) = 0;
virtual NekDouble DistanceTo(Array<OneD, NekDouble> xyz) = 0;
protected:
/// Length of edge
......
......@@ -70,7 +70,8 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
{
NekDouble t = bnds[0] + dt * k;
Array<OneD, NekDouble> l = ein[i]->edges[j]->P(t);
Array<OneD, NekDouble> uv = surf->locuv(l);
Array<OneD, NekDouble> uv(2);
surf->locuv(l, uv);
loop.push_back(uv);
}
}
......@@ -80,7 +81,8 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
{
NekDouble t = bnds[0] + dt * k;
Array<OneD, NekDouble> l = ein[i]->edges[j]->P(t);
Array<OneD, NekDouble> uv = surf->locuv(l);
Array<OneD, NekDouble> uv(2);
surf->locuv(l, uv);
loop.push_back(uv);
}
}
......
......@@ -127,34 +127,20 @@ public:
/**
* @brief Performs a reverse look up to find u,v and x,y,z.
* if xyz is off the surface it will return the closest point
*
* @param p Array of xyz location
* @return The parametric location of xyz on this surface
*/
virtual Array<OneD, NekDouble> locuv(Array<OneD, NekDouble> p) = 0;
virtual NekDouble locuv(Array<OneD, NekDouble> p, Array<OneD, NekDouble> &uv) = 0;
/**
* @brief does unconstrained locuv to project point from anywhere
* and calculate the distance between the orthonormal projection to the
* surface
* and the point
*/
virtual NekDouble DistanceTo(Array<OneD, NekDouble> p) = 0;
/**
* @brief takes a point from anywhere find the nearest surface point and its
* uv
*/
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;
/**
......@@ -164,12 +150,12 @@ public:
{
return m_orientation;
}
void SetName(std::string i)
{
m_name = i;
}
std::string GetName()
{
return m_name;
......@@ -181,7 +167,7 @@ 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;
};
......
......@@ -169,17 +169,6 @@ public:
return m_curves.size();
}
/**
* @brief Gets a curve from the map.
*/
CADVertSharedPtr GetVert(int i)
{
std::map<int, CADVertSharedPtr>::iterator search = m_verts.find(i);
ASSERTL0(search != m_verts.end(), "vert does not exist");
return search->second;
}
/**
* @brief Gets a curve from the map.
*/
......@@ -201,7 +190,10 @@ public:
return search->second;
}
/**
* @brief Gets a vert from the map.
*/
CADVertSharedPtr GetVert(int i)
{
auto search = m_verts.find(i);
......@@ -225,9 +217,9 @@ public:
{
return m_verts.size();
}
std::string GetSurfaceName(int i);
NEKMESHUTILS_EXPORT Array<OneD, NekDouble> GetPeriodicTranslationVector(
int first, int second);
......
......@@ -49,11 +49,11 @@ namespace NekMeshUtils
//forward decleration
class Node;
typedef boost::shared_ptr<Node> NodeSharedPtr;
typedef std::shared_ptr<Node> NodeSharedPtr;
class CADSurf;
typedef boost::shared_ptr<CADSurf> CADSurfSharedPtr;
typedef std::shared_ptr<CADSurf> CADSurfSharedPtr;
class CADCurve;
typedef boost::shared_ptr<CADCurve> CADCurveSharedPtr;
typedef std::shared_ptr<CADCurve> CADCurveSharedPtr;
/**
* @brief base class for CAD verticies.
......@@ -106,10 +106,8 @@ public:
return -1;
}
}
virtual NekDouble DistanceTo(Array<OneD, NekDouble> l) = 0;
NekDouble DistanceTo(Array<OneD, NekDouble> xyz);
virtual NekDouble DistanceTo(Array<OneD, NekDouble> l) = 0;
void AddAdjCurve(CADCurveSharedPtr c)
{
......
......@@ -81,7 +81,7 @@ NekDouble CADCurveCFI::tAtArcLength(NekDouble s)
return t - dt;
}
NekDouble CADCurveCFI::loct(Array<OneD, NekDouble> xyz)
NekDouble CADCurveCFI::loct(Array<OneD, NekDouble> xyz, NekDouble &t)
{
cfi::Position p;
p.x = xyz[0] / m_scal;
......@@ -91,18 +91,7 @@ NekDouble CADCurveCFI::loct(Array<OneD, NekDouble> xyz)
boost::optional<cfi::Projected<double> > pj =
m_cfiEdge->calcTFromXYZ(p, -1);
return pj.value().parameters;
}
NekDouble CADCurveCFI::DistanceTo(Array<OneD, NekDouble> xyz)
{
cfi::Position p;
p.x = xyz[0] / m_scal;
p.y = xyz[1] / m_scal;
p.z = xyz[2] / m_scal;
boost::optional<cfi::Projected<double> > pj =
m_cfiEdge->calcTFromXYZ(p, -1);
t = pj.value().parameters;
return pj.value().distance * m_scal;
}
......
......@@ -82,8 +82,7 @@ public:
}
virtual NekDouble tAtArcLength(NekDouble s);
virtual Array<OneD, NekDouble> GetMinMax();
virtual NekDouble loct(Array<OneD, NekDouble> xyz);
virtual NekDouble DistanceTo(Array<OneD, NekDouble> xyz);
virtual NekDouble loct(Array<OneD, NekDouble> xyz, NekDouble &t);
void Initialise(int i, cfi::Line *in, NekDouble s);
......
......@@ -65,7 +65,7 @@ Array<OneD, NekDouble> CADSurfCFI::GetBounds()
return b;
}
Array<OneD, NekDouble> CADSurfCFI::locuv(Array<OneD, NekDouble> p)
NekDouble CADSurfCFI::locuv(Array<OneD, NekDouble> p, Array<OneD, NekDouble> &uv)
{
cfi::Position px;
px.x = p[0] / m_scal;
......@@ -74,7 +74,6 @@ Array<OneD, NekDouble> CADSurfCFI::locuv(Array<OneD, NekDouble> p)
boost::optional<cfi::UVPosition> r = m_cfiSurface->calcUVFromXYZ(px);
Array<OneD, NekDouble> uv(2);
uv[0] = r.value().u;
uv[1] = r.value().v;
......@@ -83,7 +82,8 @@ Array<OneD, NekDouble> CADSurfCFI::locuv(Array<OneD, NekDouble> p)
NekDouble dist =
sqrt((p[0] - p2[0]) * (p[0] - p2[0]) + (p[1] - p2[1]) * (p[1] - p2[1]) +
(p[2] - p2[2]) * (p[2] - p2[2]));
return uv;
return dist * m_scal;
}
NekDouble CADSurfCFI::Curvature(Array<OneD, NekDouble> uv)
......@@ -94,42 +94,6 @@ NekDouble CADSurfCFI::Curvature(Array<OneD, NekDouble> uv)
return mxp.maxCurv.curvature * m_scal;
}
NekDouble CADSurfCFI::DistanceTo(Array<OneD, NekDouble> p)
{
cfi::Position px;
px.x = p[0] / m_scal;
px.y = p[1] / m_scal;
px.z = p[2] / m_scal;
boost::optional<cfi::UVPosition> r = m_cfiSurface->calcUVFromXYZ(px);
Array<OneD, NekDouble> uv(2);
uv[0] = r.value().u;
uv[1] = r.value().v;
Array<OneD, NekDouble> p2 = P(uv);
NekDouble dist =
sqrt((p[0] - p2[0]) * (p[0] - p2[0]) + (p[1] - p2[1]) * (p[1] - p2[1]) +
(p[2] - p2[2]) * (p[2] - p2[2]));
return dist * m_scal;
}
void CADSurfCFI::ProjectTo(Array<OneD, NekDouble> &tp,
Array<OneD, NekDouble> &uv)
{
cfi::Position px;
px.x = tp[0] / m_scal;
px.y = tp[1] / m_scal;
px.z = tp[2] / m_scal;
boost::optional<cfi::UVPosition> r = m_cfiSurface->calcUVFromXYZ(px);
uv[0] = r.value().u;
uv[1] = r.value().v;
}
Array<OneD, NekDouble> CADSurfCFI::P(Array<OneD, NekDouble> uv)
{
cfi::UVPosition uvp(uv[0], uv[1]);
......
......@@ -74,13 +74,21 @@ public:
Array<OneD, NekDouble> P(Array<OneD, NekDouble> uv);
Array<OneD, NekDouble> locuv(Array<OneD, NekDouble> p);
NekDouble locuv(Array<OneD, NekDouble> p, Array<OneD, NekDouble> &uv);
NekDouble DistanceTo(Array<OneD, NekDouble> p);
NekDouble Curvature(Array<OneD, NekDouble> uv);
void ProjectTo(Array<OneD, NekDouble> &tp, Array<OneD, NekDouble> &uv);
Array<OneD, NekDouble> BoundingBox()
{
ASSERTL0(false, "Not implemented in CFI");
return Array<OneD, NekDouble>();
}
NekDouble Curvature(Array<OneD, NekDouble> uv);
bool IsPlanar()
{
ASSERTL0(false, "Not implemented in CFI");
return false;
}
private:
/// Function which tests the the value of uv used is within the surface
......
......@@ -102,7 +102,7 @@ private:
NekDouble m_scal;
};
typedef boost::shared_ptr<CADSystemCFI> CADSystemCFISharedPtr;
typedef std::shared_ptr<CADSystemCFI> CADSystemCFISharedPtr;
}
}
......
......@@ -56,7 +56,7 @@ void CADVertCFI::Initialise(int i, cfi::Point* in, NekDouble s)
cfi::Position pos = m_cfipoint->getGeometry();
m_node = boost::shared_ptr<Node>(
m_node = std::shared_ptr<Node>(
new Node(i - 1, pos.x*m_scal, pos.y*m_scal, pos.z*m_scal));
degen = false;
}
......
......@@ -67,6 +67,11 @@ public:
void Initialise(int i, cfi::Point* in, NekDouble s);
NekDouble DistanceTo(Array<OneD, NekDouble> l)
{
ASSERTL0(false, "Not implemented in CFI");
return 0;
}
private:
/// cfi object
......
......@@ -78,31 +78,16 @@ NekDouble CADCurveOCE::Length(NekDouble ti, NekDouble tf)
return System.Mass() / 1000.0;
}
NekDouble CADCurveOCE::loct(Array<OneD, NekDouble> xyz)
NekDouble CADCurveOCE::loct(Array<OneD, NekDouble> xyz, NekDouble &t)
{
NekDouble t = 0.0;
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);
WARNINGL2(p.Distance(loc) < 1e-3, "large locuv distance " +
boost::lexical_cast<string>(p.Distance(loc)/1000.0));
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;
}
......@@ -145,7 +130,8 @@ Array<OneD, NekDouble> CADCurveOCE::NormalWRT(NekDouble t, int surf)
ASSERTL0(m_adjSurfs.size() == 1, "This will only work in 2D for one surface at the moment");
surface = m_adjSurfs[0];
Array<OneD, NekDouble> uv = surface.first->locuv(p);
Array<OneD, NekDouble> uv;
surface.first->locuv(p, uv);
Array<OneD, NekDouble> d1 = surface.first->D1(uv);
NekDouble t1 = t - 1e-8;
......@@ -156,8 +142,10 @@ Array<OneD, NekDouble> CADCurveOCE::NormalWRT(NekDouble t, int surf)
swap(t1, t2);
}
Array<OneD, NekDouble> uv1 = surface.first->locuv(P(t1));
Array<OneD, NekDouble> uv2 = surface.first->locuv(P(t2));
Array<OneD, NekDouble> uv1;
surface.first->locuv(P(t1), uv1);
Array<OneD, NekDouble> uv2;
surface.first->locuv(P(t2), uv2);
NekDouble du = uv2[1] - uv1[1];
NekDouble dv = -1.0*(uv2[0] - uv1[0]);
......
......@@ -69,11 +69,10 @@ public:
virtual Array<OneD, NekDouble> D2(NekDouble t);
virtual NekDouble tAtArcLength(NekDouble s);
virtual Array<OneD, NekDouble> GetMinMax();
virtual NekDouble loct(Array<OneD, NekDouble> xyz);
virtual NekDouble loct(Array<OneD, NekDouble> xyz, NekDouble &t);
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,9 +71,9 @@ 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);
}
......@@ -88,19 +88,19 @@ bool CADSurfOCE::IsPlanar()
{
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;
......@@ -111,50 +111,34 @@ Array<OneD, NekDouble> CADSurfOCE::BoundingBox()
return ret;
}
Array<OneD, NekDouble> CADSurfOCE::locuv(Array<OneD, NekDouble> p)
NekDouble CADSurfOCE::locuv(Array<OneD, NekDouble> p, Array<OneD, NekDouble>& uv)
{
// has to transfer back to mm
gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
Array<OneD, NekDouble> uvr(2);
uv = Array<OneD, NekDouble>(2);
gp_Pnt2d p2 = m_sas->ValueOfUV(loc, Precision::Confusion());
uvr[0] = p2.X();
uvr[1] = p2.Y();
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
//(which ShapeAnalysis_Surface can do sometimes)
if (uvr[0] < m_bounds[0] || uvr[0] > m_bounds[1] || uvr[1] < m_bounds[2] ||
uvr[1] > m_bounds[3])
TopAbs_State s = m_2Dclass->Perform(p2);
if(s == TopAbs_OUT)
{
if (uvr[0] < m_bounds[0])
{
uvr[0] = m_bounds[0];
}
else if (uvr[0] > m_bounds[1])
{
uvr[0] = m_bounds[1];
}
else if (uvr[1] < m_bounds[2])
{
uvr[1] = m_bounds[2];
}
else if (uvr[1] > m_bounds[3])
{
uvr[1] = m_bounds[3];
}
else
{
ASSERTL0(false, "Cannot correct locuv");
}
BRepBuilderAPI_MakeVertex v(loc);
BRepExtrema_DistShapeShape dss(BRepTools::OuterWire(TopoDS::Face(m_shape)), v.Shape());
dss.Perform();
gp_Pnt np = dss.PointOnShape1(1);
gp_Pnt2d p22 = m_sas->ValueOfUV(np, Precision::Confusion());
uv[0] = p22.X();
uv[1] = p22.Y();
}
else
{
uv[0] = p2.X();
uv[1] = p2.Y();
}
return uvr;
gp_Pnt p3 = m_sas->Value(p2);
return p3.Distance(loc) / 1000.0;
}
NekDouble CADSurfOCE::Curvature(Array<OneD, NekDouble> uv)
......@@ -200,65 +184,6 @@ NekDouble CADSurfOCE::Curvature(Array<OneD, NekDouble> uv)
return kv[0] > kv[1] ? kv[0] : kv[1];
}
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);