Commit f558526f authored by Michael Turner's avatar Michael Turner
Browse files

basic quad adding for sufaces, not quite working

parent 22fef1b5
......@@ -203,6 +203,29 @@ NekDouble CADSurf::DistanceTo(Array<OneD, NekDouble> p)
return p3.Distance(loc);
}
void CADSurf::ProjectTo(Array<OneD, NekDouble> &tp, Array<OneD, NekDouble> &uv)
{
gp_Pnt loc(tp[0] * 1000.0, tp[1] * 1000.0, tp[2] * 1000.0);
// alternative locuv methods
ShapeAnalysis_Surface sas(m_s);
sas.SetDomain(m_occSurface.FirstUParameter(),
m_occSurface.LastUParameter(),
m_occSurface.FirstVParameter(),
m_occSurface.LastVParameter());
gp_Pnt2d p2 = sas.ValueOfUV(loc,1e-7);
gp_Pnt p3 = sas.Value(p2);
tp[0] = p3.X()/1000.0;
tp[1] = p3.Y()/1000.0;
tp[2] = p3.Z()/1000.0;
uv[0] = p2.X();
uv[1] = p2.X();
}
bool CADSurf::IsPlane()
{
if(m_occSurface.GetType() == GeomAbs_Plane)
......
......@@ -180,6 +180,8 @@ public:
*/
NekDouble DistanceTo(Array<OneD, NekDouble> p);
void ProjectTo(Array<OneD, NekDouble> &tp, Array<OneD, NekDouble> &uv);
/**
* @brief returns curvature at point uv
*/
......
......@@ -207,33 +207,6 @@ bool CADSystem::LoadCAD()
}
}
//attempts to identify properties of the vertex on the degen edge
for(int i = 1; i <= mapOfFaces.Extent(); i++)
{
TopoDS_Shape face= mapOfFaces.FindKey(i);
TopTools_IndexedMapOfShape localEdges;
TopExp::MapShapes(face, TopAbs_EDGE, localEdges);
for(int j = 1; j <= localEdges.Extent(); j++)
{
TopoDS_Shape edge = localEdges.FindKey(j);
if(BRep_Tool::Degenerated(TopoDS::Edge(edge)))
{
cout << "degen edge on face " << i << endl;
gp_Pnt2d p1,p2;
BRep_Tool::UVPoints(TopoDS::Edge(edge),TopoDS::Face(face),p1,p2);
m_verts[mapOfVerts.FindIndex(TopExp::FirstVertex(
TopoDS::Edge(edge), Standard_True))]->SetDegen(
i, (p1.X() + p2.X())/2.0, (p1.Y() + p2.Y())/2.0);
}
}
}
map<int, vector<int> > adjsurfmap; //from id of curve to list of ids of surfs
// Adds edges to our type and map
......@@ -358,6 +331,31 @@ bool CADSystem::LoadCAD()
AddSurf(i, face, edgeloops);
}
//attempts to identify properties of the vertex on the degen edge
for(int i = 1; i <= mapOfFaces.Extent(); i++)
{
TopoDS_Shape face= mapOfFaces.FindKey(i);
TopTools_IndexedMapOfShape localEdges;
TopExp::MapShapes(face, TopAbs_EDGE, localEdges);
for(int j = 1; j <= localEdges.Extent(); j++)
{
TopoDS_Shape edge = localEdges.FindKey(j);
if(BRep_Tool::Degenerated(TopoDS::Edge(edge)))
{
cout << "degen edge on face " << i << endl;
gp_Pnt2d p1,p2;
BRep_Tool::UVPoints(TopoDS::Edge(edge),TopoDS::Face(face),p1,p2);
m_verts[mapOfVerts.FindIndex(TopExp::FirstVertex(
TopoDS::Edge(edge), Standard_True))]->SetDegen(
i, m_surfs[i], (p1.X() + p2.X())/2.0, (p1.Y() + p2.Y())/2.0);
}
}
}
// This checks that all edges are bound by two surfaces, sanity check.
for(map<int,vector<int> >::iterator it = adjsurfmap.begin();
it != adjsurfmap.end(); it++)
......
......@@ -165,7 +165,7 @@ public:
void SmallFeatureAnalysis(NekDouble min);
private:
/// Private function to add curve to CADSystem::m_verts.
void AddVert(int i, TopoDS_Shape in);
/// Private function to add curve to CADSystem::m_curves.
......
......@@ -103,18 +103,18 @@ public:
/**
* @brief if the vertex is degenerate manually set uv for that surface
*/
void SetDegen(int s, NekDouble u, NekDouble v)
void SetDegen(int s, CADSurfSharedPtr su, NekDouble u, NekDouble v)
{
degen = true;
degensurf = s;
Array<OneD, NekDouble> uv(2);
uv[0] = u;
uv[1] = v;
m_node->SetCADSurf(s,uv);
m_node->SetCADSurf(s, su, uv);
}
/**
* @brief query is degenerate
* @brief query is degenerate
*/
int IsDegen()
{
......
......@@ -91,6 +91,10 @@ IF(NEKTAR_USE_OCC)
ADD_DEPENDENCIES(NekMeshUtils opencascade-6.8)
ENDIF(NEKTAR_USE_OCC)
IF (NEKTAR_USE_MESH)
SET_TARGET_PROPERTIES(NekMeshUtils PROPERTIES COMPILE_DEFINITIONS "MESHGEN")
ENDIF()
INSTALL(DIRECTORY ./
DESTINATION ${NEKTAR_INCLUDE_DIR}/NekMeshUtils
COMPONENT dev
......
......@@ -77,7 +77,7 @@ void TriangleInterface::Mesh(bool Quiet, bool Quality)
{
nodemap[pointc] = m_boundingloops[i][j];
Array<OneD, NekDouble> uv = m_boundingloops[i][j]->GetCADSurf(sid);
Array<OneD, NekDouble> uv = m_boundingloops[i][j]->GetCADSurfInfo(sid);
in.pointlist[pointc*2+0] = uv[0]*m_str;
in.pointlist[pointc*2+1] = uv[1];
}
......@@ -87,7 +87,7 @@ void TriangleInterface::Mesh(bool Quiet, bool Quality)
{
nodemap[pointc] = m_stienerpoints[i];
Array<OneD, NekDouble> uv = m_stienerpoints[i]->GetCADSurf(sid);
Array<OneD, NekDouble> uv = m_stienerpoints[i]->GetCADSurfInfo(sid);
in.pointlist[pointc*2+0] = uv[0]*m_str;
in.pointlist[pointc*2+1] = uv[1];
}
......
......@@ -53,16 +53,27 @@ namespace NekMeshUtils
std::vector<NodeSharedPtr> pEdgeNodes,
LibUtilities::PointsType pCurveType)
: m_n1(pVertex1), m_n2(pVertex2), m_edgeNodes(pEdgeNodes),
m_curveType(pCurveType), CADCurveID(-1), m_geom() {}
m_curveType(pCurveType), m_geom()
{
#ifdef MESHGEN
onCurve = false;
#endif
}
/// Creates a new linear edge.
Edge(NodeSharedPtr pVertex1, NodeSharedPtr pVertex2)
: m_n1(pVertex1), m_n2(pVertex2), m_edgeNodes(),
m_curveType(), CADCurveID(-1), m_geom() {}
m_curveType(), m_geom()
{
#ifdef MESHGEN
onCurve = false;
#endif
}
/// Copies an existing edge.
Edge(const Edge& pSrc)
: m_n1(pSrc.m_n1), m_n2(pSrc.m_n2), m_edgeNodes(pSrc.m_edgeNodes),
m_curveType(pSrc.m_curveType), CADCurveID(pSrc.CADCurveID), m_geom(pSrc.m_geom) {}
~Edge() {}
m_curveType(pSrc.m_curveType), m_geom(pSrc.m_geom) {}
~Edge()
{}
/// Returns the total number of nodes defining the edge.
unsigned int GetNodeCount() const
......@@ -136,9 +147,13 @@ namespace NekMeshUtils
LibUtilities::PointsType m_curveType;
/// Element(s) which are linked to this edge.
vector<pair<ElementSharedPtr, int> > m_elLink;
#ifdef MESHGEN
bool onCurve;
/// id of cad curve which edge lies on
int CADCurveID;
std::vector<int> CADSurfID;
int CADCurveId;
CADCurveSharedPtr CADCurve;
#endif
private:
SpatialDomains::SegGeomSharedPtr m_geom;
......
......@@ -124,14 +124,6 @@ namespace NekMeshUtils
if (m_edgeLink.get() != 0) return m_edgeLink->m_id;
return m_id;
}
/// returns id of the Tri which this elemnt was created from in
/// mesh utils, only applys to surface elements
int GetCADSurf() const {
return m_meshtriID;
}
void SetCADSurf(int i) {
m_meshtriID = i;
}
/// Returns the expansion dimension of the element.
unsigned int GetDim() const {
return m_dim;
......@@ -446,11 +438,13 @@ namespace NekMeshUtils
}
}
#ifdef MESHGEN
int CADSurfId;
#endif
protected:
/// ID of the element.
unsigned int m_id;
/// ID of mesh utilse Tri
int m_meshtriID;
/// Dimension of the element.
unsigned int m_dim;
/// Contains configuration of the element.
......
......@@ -49,6 +49,19 @@
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/ShapeType.hpp>
#ifdef MESHGEN //forard declare cad objects for mesh elements
namespace Nektar
{
namespace NekMeshUtils
{
class CADCurve;
class CADSurf;
typedef boost::shared_ptr<CADCurve> CADCurveSharedPtr;
typedef boost::shared_ptr<CADSurf> CADSurfSharedPtr;
}
}
#endif
#include <SpatialDomains/SegGeom.h>
#include <SpatialDomains/TriGeom.h>
#include <SpatialDomains/QuadGeom.h>
......
......@@ -153,16 +153,6 @@ namespace NekMeshUtils
m_z*pSrc.m_x-m_x*pSrc.m_z, m_x*pSrc.m_y-m_y*pSrc.m_x);
}
void SetCADCurve(int i, NekDouble t)
{
CADCurve[i] = t;
}
void SetCADSurf(int i, Array<OneD, NekDouble> uv)
{
CADSurf[i] = uv;
}
Array<OneD, NekDouble> GetLoc()
{
Array<OneD, NekDouble> out(3);
......@@ -170,46 +160,14 @@ namespace NekMeshUtils
return out;
}
NekDouble GetCADCurve(int i)
{
std::map<int, NekDouble>::iterator search =
CADCurve.find(i);
ASSERTL0(search->first == i, "node not on this curve");
return search->second;
}
Array<OneD, NekDouble> GetCADSurf(int i)
{
//I dont know why I ahev to do this to get it to work
//this really needs bound checking
std::map<int, Array<OneD, NekDouble> >::iterator search =
CADSurf.find(i);
ASSERTL0(search->first == i,"surface not found");
return search->second;
}
std::vector<int> GetListCADCurve()
/// Generate a %SpatialDomains::PointGeom for this node.
SpatialDomains::PointGeomSharedPtr GetGeom(int coordDim)
{
std::vector<int> list;
std::map<int, NekDouble >::iterator c;
for(c = CADCurve.begin(); c != CADCurve.end(); c++)
{
list.push_back(c->first);
}
return list;
}
SpatialDomains::PointGeomSharedPtr ret =
MemoryManager<SpatialDomains::PointGeom>
::AllocateSharedPtr(coordDim,m_id,m_x,m_y,m_z);
std::vector<int> GetListCADSurf()
{
std::vector<int> list;
std::map<int, Array<OneD, NekDouble> >::iterator s;
for(s = CADSurf.begin(); s != CADSurf.end(); s++)
{
list.push_back(s->first);
}
return list;
return ret;
}
NekDouble Distance(NodeSharedPtr &p)
......@@ -251,22 +209,78 @@ namespace NekMeshUtils
return an;
}
void Move(Array<OneD, NekDouble> l, int s, Array<OneD, NekDouble> uv)
#ifdef MESHGEN //fucntions for cad information
void SetCADCurve(int i, CADCurveSharedPtr c, NekDouble t)
{
m_x = l[0]; m_y = l[1]; m_z = l[2];
CADSurf[s] = uv;
CADCurveList[i] = std::pair<CADCurveSharedPtr, NekDouble>(c,t);
}
/// Generate a %SpatialDomains::PointGeom for this node.
SpatialDomains::PointGeomSharedPtr GetGeom(int coordDim)
void SetCADSurf(int i, CADSurfSharedPtr s, Array<OneD, NekDouble> uv)
{
SpatialDomains::PointGeomSharedPtr ret =
MemoryManager<SpatialDomains::PointGeom>
::AllocateSharedPtr(coordDim,m_id,m_x,m_y,m_z);
CADSurfList[i] = std::pair<CADSurfSharedPtr, Array<OneD, NekDouble> >(s,uv);
}
return ret;
CADCurveSharedPtr GetCADCurve(int i)
{
std::map<int, std::pair<CADCurveSharedPtr, NekDouble> >::iterator search =
CADCurveList.find(i);
ASSERTL0(search->first == i, "node not on this curve");
return search->second.first;
}
CADSurfSharedPtr GetCADSurf(int i)
{
std::map<int, std::pair<CADSurfSharedPtr, Array<OneD, NekDouble> > >::iterator search =
CADSurfList.find(i);
ASSERTL0(search->first == i,"surface not found");
return search->second.first;
}
NekDouble GetCADCurveInfo(int i)
{
std::map<int, std::pair<CADCurveSharedPtr, NekDouble> >::iterator search =
CADCurveList.find(i);
ASSERTL0(search->first == i, "node not on this curve");
return search->second.second;
}
Array<OneD, NekDouble> GetCADSurfInfo(int i)
{
std::map<int, std::pair<CADSurfSharedPtr, Array<OneD, NekDouble> > >::iterator search =
CADSurfList.find(i);
ASSERTL0(search->first == i,"surface not found");
return search->second.second;
}
std::vector<std::pair<int, CADSurfSharedPtr> > GetCADSurfs()
{
std::vector<std::pair<int, CADSurfSharedPtr> > lst;
std::map<int, std::pair<CADSurfSharedPtr, Array<OneD, NekDouble> > >::iterator s;
for(s = CADSurfList.begin(); s != CADSurfList.end(); s++)
{
lst.push_back(std::pair<int, CADSurfSharedPtr>(s->first,s->second.first));
}
return lst;
}
int GetNumCadCurve() {return CADCurveList.size();}
int GetNumCADSurf() {return CADSurfList.size();}
void Move(Array<OneD, NekDouble> l, int s, Array<OneD, NekDouble> uv)
{
m_x = l[0]; m_y = l[1]; m_z = l[2];
CADSurfSharedPtr su = CADSurfList[s].first;
CADSurfList[s] = std::pair<CADSurfSharedPtr, Array<OneD, NekDouble> >(su,uv);
}
#endif
/// ID of node.
int m_id;
/// X-coordinate.
......@@ -275,10 +289,14 @@ namespace NekMeshUtils
NekDouble m_y;
/// Z-coordinate.
NekDouble m_z;
#ifdef MESHGEN //tag to tell the meshelemnets to include cad information
///list of cadcurves the node lies on
std::map<int, NekDouble> CADCurve;
std::map<int, std::pair<CADCurveSharedPtr, NekDouble> > CADCurveList;
///list of cadsurfs the node lies on
std::map<int, Array<OneD, NekDouble> > CADSurf;
std::map<int, std::pair<CADSurfSharedPtr, Array<OneD, NekDouble> > > CADSurfList;
#endif
private:
SpatialDomains::PointGeomSharedPtr m_geom;
......
......@@ -1203,7 +1203,6 @@ void Octree::AssignNeigbours(OctantSharedPtr const &o)
}
}
o->SetNeigbourList(nlist);
ASSERTL0(nlist.size()>6,"not enough neigbours");
ASSERTL0(eqhit == 1,"found itself more than once");
}
......
......@@ -111,7 +111,7 @@ void CurveMesh::Mesh()
NodeSharedPtr n = verts[0]->GetNode();
t = m_bounds[0];
n->SetCADCurve(m_id,t);
n->SetCADCurve(m_id, m_cadcurve,t);
loc = n->GetLoc();
for(int j = 0; j < 2; j++)
{
......@@ -119,7 +119,7 @@ void CurveMesh::Mesh()
continue;
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
n->SetCADSurf(s[j]->GetId(), uv);
n->SetCADSurf(s[j]->GetId(), s[j], uv);
}
m_meshpoints.push_back(n);
......@@ -130,18 +130,18 @@ 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,t);
n2->SetCADCurve(m_id, m_cadcurve, t);
for(int j = 0; j < 2; j++)
{
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
n2->SetCADSurf(s[j]->GetId(), uv);
n2->SetCADSurf(s[j]->GetId(), s[j], uv);
}
m_meshpoints.push_back(n2);
}
n = verts[1]->GetNode();
t = m_bounds[1];
n->SetCADCurve(m_id,t);
n->SetCADCurve(m_id, m_cadcurve, t);
loc = n->GetLoc();
for(int j = 0; j < 2; j++)
{
......@@ -149,7 +149,7 @@ void CurveMesh::Mesh()
continue;
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
n->SetCADSurf(s[j]->GetId(), uv);
n->SetCADSurf(s[j]->GetId(), s[j], uv);
}
m_meshpoints.push_back(n);
......@@ -161,7 +161,7 @@ void CurveMesh::Mesh()
for(int j = 0; j < 2; j++)
{
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
m_meshpoints[i]->SetCADSurf(s[j]->GetId(), uv);
m_meshpoints[i]->SetCADSurf(s[j]->GetId(), s[j], uv);
}
}
......@@ -176,8 +176,8 @@ void CurveMesh::Mesh()
for(int j = 0; j < 2; j++)
{
Array<OneD, NekDouble> uv1, uv2;
uv1 = m_meshpoints[i]->GetCADSurf(s[j]->GetId());
uv2 = m_meshpoints[i+1]->GetCADSurf(s[j]->GetId());
uv1 = m_meshpoints[i]->GetCADSurfInfo(s[j]->GetId());
uv2 = m_meshpoints[i+1]->GetCADSurfInfo(s[j]->GetId());
Array<OneD, NekDouble> N1, N2;
N1 = s[j]->N(uv1);
N2 = s[j]->N(uv2);
......@@ -192,17 +192,17 @@ void CurveMesh::Mesh()
{
ct++;
NekDouble t1, t2;
t1 = m_meshpoints[i]->GetCADCurve(m_id);
t2 = m_meshpoints[i+1]->GetCADCurve(m_id);
t1 = m_meshpoints[i]->GetCADCurveInfo(m_id);
t2 = m_meshpoints[i+1]->GetCADCurveInfo(m_id);
NekDouble tn = (t1 + t2)/2.0;
Array<OneD, NekDouble> loc = m_cadcurve->P(tn);
NodeSharedPtr nn = boost::shared_ptr<Node>(new Node(m_mesh->m_numNodes++,
loc[0],loc[1],loc[2]));
nn->SetCADCurve(m_id, tn);
nn->SetCADCurve(m_id, m_cadcurve, tn);
for(int j = 0; j < 2; j++)
{
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
nn->SetCADSurf(s[j]->GetId(), uv);
nn->SetCADSurf(s[j]->GetId(), s[j], uv);
}
m_meshpoints.insert(m_meshpoints.begin() + i+1, nn);
break;
......@@ -215,7 +215,9 @@ void CurveMesh::Mesh()
{
EdgeSharedPtr e = boost::shared_ptr<Edge>(new Edge(m_meshpoints[i],
m_meshpoints[i+1]));
e->CADCurveID = m_id;
e->CADCurveId = m_id;
e->CADCurve = m_cadcurve;
e->onCurve = true;
m_mesh->m_edgeSet.insert(e);
}
......
......@@ -53,7 +53,10 @@ void FaceMesh::Mesh()
OrientateCurves();
if(m_makebl)
{
MakeBL();
cout << "done" << endl;
}
int numPoints = 0;
for(int i = 0; i < orderedLoops.size(); i++)
......@@ -109,7 +112,7 @@ void FaceMesh::Mesh()
BuildLocalMesh();
OptimiseLocalMesh();
//OptimiseLocalMesh();
//clear local element links
EdgeSet::iterator eit;
......@@ -146,7 +149,38 @@ void FaceMesh::MakeBL()
for(int j = 0; j < orderedLoops[i].size(); j++)
{
//for each of the nodes make a new node which exists off the surface
vector<pair<int, CADSurfSharedPtr> > surfs = orderedLoops[i][j]->GetCADSurfs();
ASSERTL0(surfs.size() > 1, "point does not have enough surfs to make quad");
Array<OneD, NekDouble> AN(3,0.0);
//make a averaged normal ignoring surface mid
for(int s = 0; s < surfs.size(); s++)
{
if(surfs[s].first == m_id)
continue; <