Commit ff0ae48e authored by Michael Turner's avatar Michael Turner

most of linear meshing working, high-order a bit off with a surf finding error

parent dd07e528
......@@ -173,6 +173,32 @@ 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
......
......@@ -43,6 +43,8 @@
#include <MeshUtils/CADSystem/OpenCascade.h>
#include <MeshUtils/MeshElements/MeshElements.h>
namespace Nektar
{
namespace MeshUtils
......@@ -72,6 +74,10 @@ public:
m_ID = i;
m_occVert = BRep_Tool::Pnt(TopoDS::Vertex(in));
m_node = boost::shared_ptr<Node>(
new Node(0, m_occVert.X(), m_occVert.Y(), m_occVert.Z()));
degen = false;
}
Array<OneD, NekDouble> GetLoc()
......@@ -83,11 +89,41 @@ public:
int GetId(){return m_ID;}
NodeSharedPtr GetNode(){return m_node;}
void SetDegen(int s, NekDouble u, NekDouble v)
{
degen = true;
degensurf = s;
Array<OneD, NekDouble> uv(2);
uv[0] = u;
uv[1] = v;
m_node->SetCADSurf(s,uv);
}
int IsDegen()
{
if(degen)
{
return degensurf;
}
else
{
return -1;
}
}
private:
/// ID of the vert.
int m_ID;
/// OpenCascade object of the curve.
gp_Pnt m_occVert;
/// mesh convert object of vert
NodeSharedPtr m_node;
///degen marker
bool degen;
//degen surface
int degensurf;
};
typedef boost::shared_ptr<CADVert> CADVertSharedPtr;
......
......@@ -73,45 +73,38 @@ void TriangleInterface::Mesh(bool Quiet, bool Quality)
for(int i = 0; i < m_boundingloops.size(); i++)
{
for(int j = 0; j < m_boundingloops[i].size(); j++)
for(int j = 0; j < m_boundingloops[i].size(); j++, pointc++)
{
nodemap[pointc] = m_boundingloops[i][j]->m_id;
nodemapr[m_boundingloops[i][j]->m_id] = pointc;
nodemap[pointc] = m_boundingloops[i][j];
Array<OneD, NekDouble> uv = m_boundingloops[i][j]->GetCADSurf(sid);
in.pointlist[pointc*2+0] = uv[0]*m_str;
in.pointlist[pointc*2+1] = uv[1];
pointc++;
}
}
for(int i = 0; i < m_stienerpoints.size(); i++)
for(int i = 0; i < m_stienerpoints.size(); i++, pointc++)
{
nodemap[pointc] = m_stienerpoints[i]->m_id;
nodemapr[m_stienerpoints[i]->m_id] = pointc;
nodemap[pointc] = m_stienerpoints[i];
Array<OneD, NekDouble> uv = m_stienerpoints[i]->GetCADSurf(sid);
in.pointlist[pointc*2+0] = uv[0]*m_str;
in.pointlist[pointc*2+1] = uv[1];
pointc++;
}
in.numberofsegments = numSeg;
in.segmentlist = (int *) malloc(in.numberofsegments*2*sizeof(int));
pointc=0;
for(int i = 0; i < m_boundingloops.size(); i++)
for(int i = 0; i < m_boundingloops.size(); i++, pointc++)
{
for(int j = 0; j < m_boundingloops[i].size()-1; j++)
int first = pointc;
for(int j = 0; j < m_boundingloops[i].size()-1; j++, pointc++)
{
in.segmentlist[pointc*2+0] = nodemapr[m_boundingloops[i][j]->m_id];
in.segmentlist[pointc*2+1] = nodemapr[m_boundingloops[i][j+1]->m_id];
pointc++;
in.segmentlist[pointc*2+0] = pointc;
in.segmentlist[pointc*2+1] = pointc+1;
}
in.segmentlist[pointc*2+0] = nodemapr[m_boundingloops[i].back()->m_id];
in.segmentlist[pointc*2+1] = nodemapr[m_boundingloops[i][0]->m_id];
pointc++;
in.segmentlist[pointc*2+0] = pointc;
in.segmentlist[pointc*2+1] = first;
}
in.numberofregions = 0;
......@@ -142,7 +135,7 @@ void TriangleInterface::Mesh(bool Quiet, bool Quality)
}
else if(Quiet && !Quality)
{
triangulate("pzenQYY", &in, &out, NULL);
triangulate("pzenYYQ", &in, &out, NULL);
}
else if(!Quiet && Quality)
{
......@@ -153,11 +146,6 @@ void TriangleInterface::Mesh(bool Quiet, bool Quality)
triangulate("pzenYY", &in, &out, NULL);
}
for(int i = 0; i < out.numberofpoints; i++)
{
out.pointlist[i*2+0] = out.pointlist[i*2+0]/m_str;
}
//verify the mesh a bit
if(out.numberofpoints-out.numberofedges+out.numberoftriangles != 2 - m_centers.size())
{
......@@ -167,12 +155,12 @@ void TriangleInterface::Mesh(bool Quiet, bool Quality)
}
void TriangleInterface::Extract(std::vector<Array<OneD, int> > &Connec)
void TriangleInterface::Extract(std::vector<std::vector<NodeSharedPtr> > &Connec)
{
Connec.clear();
for(int i = 0; i < out.numberoftriangles; i++)
{
Array<OneD, int> tri(3);
vector<NodeSharedPtr> tri(3);
tri[0] = nodemap[out.trianglelist[i*3+0]];
tri[1] = nodemap[out.trianglelist[i*3+1]];
tri[2] = nodemap[out.trianglelist[i*3+2]];
......
......@@ -77,8 +77,7 @@ public:
* @brief assign meshing paramters
*/
void Assign(std::vector<std::vector<NodeSharedPtr> > &boundingloops,
std::vector<Array<OneD, NekDouble> > &centers,
std::vector<NodeSharedPtr> stiner, int i,
std::vector<Array<OneD, NekDouble> > &centers, int i,
NekDouble str = 1.0)
{
if(meshloaded)
......@@ -87,7 +86,6 @@ public:
}
m_boundingloops = boundingloops;
m_centers = centers;
m_stienerpoints = stiner;
m_str = str;
sid = i;
}
......@@ -114,7 +112,7 @@ public:
/**
* @brief extract mesh
*/
void Extract(std::vector<Array<OneD, int> > &Connec);
void Extract(std::vector<std::vector<NodeSharedPtr> > &Connec);
private:
......@@ -130,9 +128,7 @@ private:
/// coordinates of the centers of the loops
std::vector<Array<OneD, NekDouble> > m_centers;
/// map from meshconvert id to triangle id
std::map<int, int> nodemap;
/// reverse map
std::map<int, int> nodemapr;
std::map<int, NodeSharedPtr> nodemap;
/// id of the surface
int sid;
/// has mesh been run before
......
......@@ -137,8 +137,8 @@ namespace MeshUtils
/// Element(s) which are linked to this edge.
vector<pair<ElementSharedPtr, int> > m_elLink;
/// id of cad curve which edge lies on
unsigned int CADCurveID;
unsigned int CADSurfID;
int CADCurveID;
int CADSurfID;
private:
SpatialDomains::SegGeomSharedPtr m_geom;
......
......@@ -511,8 +511,6 @@ namespace MeshUtils
}
};
bool operator==(ElementSharedPtr const &e1, ElementSharedPtr const &e2);
}
}
......
......@@ -87,8 +87,6 @@ namespace MeshUtils
unsigned int m_spaceDim;
/// a order tag to aid output, a bit of a hack
unsigned int m_nummode;
///list of nodes used in meshing routines before its a full mesh
std::vector<NodeSharedPtr> m_meshnode;
/// List of mesh nodes.
std::vector<NodeSharedPtr> m_node;
/// Set of element vertices.
......
......@@ -150,18 +150,6 @@ namespace Nektar
return true;
}
bool operator==(ElementSharedPtr const &e1, ElementSharedPtr const &e2)
{
if(e1->GetVertexCount() != e2->GetVertexCount())
return false;
if(e1->GetId() != e2->GetId())
return false;
return true;
}
/**
* @brief Test equality of two conditions - i.e. compare types, fields
* and values but _not_ composite ids.
......
......@@ -61,9 +61,9 @@ namespace MeshUtils
Node(int pId, NekDouble pX, NekDouble pY, NekDouble pZ)
: m_id(pId), m_x(pX), m_y(pY), m_z(pZ), m_geom() {}
/// Copy an existing node.
Node(const Node& pSrc)
: m_id(pSrc.m_id), m_x(pSrc.m_x), m_y(pSrc.m_y),
m_z(pSrc.m_z), m_geom() {}
//Node(const Node& pSrc)
// : m_id(pSrc.m_id), m_x(pSrc.m_x), m_y(pSrc.m_y),
// m_z(pSrc.m_z), m_geom() {}
Node() : m_id(0), m_x(0.0), m_y(0.0), m_z(0.0), m_geom() {}
~Node() {}
......@@ -212,36 +212,27 @@ namespace MeshUtils
return list;
}
NekDouble Distance(NodeSharedPtr p)
NekDouble Distance(NodeSharedPtr &p)
{
return sqrt((m_x-p->m_x)*(m_x-p->m_x)+(m_y-p->m_y)*(m_y-p->m_y)
+(m_z-p->m_z)*(m_z-p->m_z));
}
NekDouble Angle(NodeSharedPtr a, NodeSharedPtr b)
NekDouble Angle(NodeSharedPtr &a, NodeSharedPtr &b)
{
Array<OneD,NekDouble> la = a->GetLoc();
Array<OneD,NekDouble> lb = b->GetLoc();
Array<OneD,NekDouble> va(3),vb(3);
va[0] = la[0] - m_x;
va[1] = la[1] - m_y;
va[2] = la[2] - m_z;
vb[0] = lb[0] - m_x;
vb[1] = lb[1] - m_y;
vb[2] = lb[2] - m_z;
va[0] = a->m_x - m_x;
va[1] = a->m_y - m_y;
va[2] = a->m_z - m_z;
vb[0] = b->m_x - m_x;
vb[1] = b->m_y - m_y;
vb[2] = b->m_z - m_z;
NekDouble ca = va[0]*vb[0] + va[1]*vb[1] + va[2]*vb[2];
ca /= sqrt(va[0]*va[0] + va[1]*va[1] + va[2]*va[2]);
ca /= sqrt(vb[0]*vb[0] + vb[1]*vb[1] + vb[2]*vb[2]);
if(ca < 0)
{
return 3.142 + acos(ca);
}
else
{
return acos(ca);
}
return acos(ca);
}
void Move(Array<OneD, NekDouble> l, int s, Array<OneD, NekDouble> uv)
......
......@@ -106,30 +106,64 @@ void CurveMesh::Mesh()
Array<OneD, NekDouble> loc;
vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();
vector<CADSurfSharedPtr> s = m_cadcurve->GetAdjSurf();
ASSERTL0(s.size() == 2, "invalid curve");
int id = verts[0]->GetId() - 1;
NodeSharedPtr n = verts[0]->GetNode();
t = m_bounds[0];
m_mesh->m_meshnode[id]->SetCADCurve(m_id,t);
m_meshpoints.push_back(m_mesh->m_meshnode[id]);
n->SetCADCurve(m_id,t);
loc = n->GetLoc();
for(int j = 0; j < 2; j++)
{
if(verts[0]->IsDegen() == s[j]->GetId()) //if the degen has been set for this node the node already knows its corrected location
continue;
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
n->SetCADSurf(s[j]->GetId(), uv);
}
m_meshpoints.push_back(n);
for(int i = 1; i < meshsvalue.size()-1; i++)
{
t = m_cadcurve->tAtArcLength(meshsvalue[i]);
loc = m_cadcurve->P(t);
NodeSharedPtr n2 = boost::shared_ptr<Node>(
new Node(m_mesh->m_meshnode.size(),loc[0],loc[1],loc[2]));
new Node(0,loc[0],loc[1],loc[2]));
n2->SetCADCurve(m_id,t);
m_mesh->m_meshnode.push_back(n2);
for(int j = 0; j < 2; j++)
{
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
n2->SetCADSurf(s[j]->GetId(), uv);
}
m_meshpoints.push_back(n2);
}
id = verts[1]->GetId() - 1;
n = verts[1]->GetNode();
t = m_bounds[1];
m_mesh->m_meshnode[id]->SetCADCurve(m_id,t);
m_meshpoints.push_back(m_mesh->m_meshnode[id]);
n->SetCADCurve(m_id,t);
loc = n->GetLoc();
for(int j = 0; j < 2; j++)
{
if(verts[1]->IsDegen() == s[j]->GetId()) //if the degen has been set for this node the node already knows its corrected location
continue;
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
n->SetCADSurf(s[j]->GetId(), uv);
}
m_meshpoints.push_back(n);
ASSERTL0(Ne+1 == m_meshpoints.size(),"incorrect number of points in curve mesh");
for(int i = 0; i < m_meshpoints.size(); i++)
{
Array<OneD, NekDouble> loc = m_meshpoints[i]->GetLoc();
for(int j = 0; j < 2; j++)
{
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
m_meshpoints[i]->SetCADSurf(s[j]->GetId(), uv);
}
}
}
......
......@@ -105,6 +105,8 @@ class CurveMesh
std::map<int, MeshNodeSharedPtr> &Nodes,
std::map<int, MeshEdgeSharedPtr> &Edges);
NekDouble GetLength(){return m_curvelength;}
private:
/**
......
......@@ -33,6 +33,7 @@
//
////////////////////////////////////////////////////////////////////////////////
#include <limits>
#include <MeshUtils/SurfaceMeshing/FaceMesh.h>
#include <MeshUtils/ExtLibInterface/TriangleInterface.h>
......@@ -80,20 +81,11 @@ void FaceMesh::Mesh()
centers.push_back(m_edgeloops[i].center);
}
pplanemesh->Assign(orderedLoops, centers, m_stienerpoints, m_id, asr/pasr);
pplanemesh->Assign(orderedLoops, centers, m_id, asr/pasr);
pplanemesh->Mesh();
vector<Array<OneD, int> > intconnec;
pplanemesh->Extract(intconnec);
for(int i = 0; i < intconnec.size(); i++)
{
vector<NodeSharedPtr> tri(3);
tri[0] = m_mesh->m_meshnode[intconnec[i][0]];
tri[1] = m_mesh->m_meshnode[intconnec[i][1]];
tri[2] = m_mesh->m_meshnode[intconnec[i][2]];
m_connec.push_back(tri);
}
pplanemesh->Extract(m_connec);
bool repeat = true;
int meshcounter = 1;
......@@ -108,17 +100,25 @@ void FaceMesh::Mesh()
m_connec.clear();
pplanemesh->AssignStiener(m_stienerpoints);
pplanemesh->Mesh();
pplanemesh->Extract(intconnec);
for(int i = 0; i < intconnec.size(); i++)
pplanemesh->Extract(m_connec);
meshcounter++;
}
NekDouble area = numeric_limits<double>::max();
//idiot check the elements
for(int i = 0; i < m_connec.size(); i++)
{
Array<OneD, NekDouble> a = m_connec[i][0]->GetCADSurf(m_id);
Array<OneD, NekDouble> b = m_connec[i][1]->GetCADSurf(m_id);
Array<OneD, NekDouble> c = m_connec[i][2]->GetCADSurf(m_id);
area = min(area, 0.5*(-b[0]*a[1] + c[0]*a[1] + a[0]*b[1] - c[0]*b[1] - a[0]*c[1] + b[0]*c[1]));
if(area <= 0)
{
vector<NodeSharedPtr> tri(3);
for(int j = 0; j < 3; j++)
{
tri[j] = m_mesh->m_meshnode[intconnec[i][j]];
}
m_connec.push_back(tri);
cout << area << endl;
exit(-1);
}
meshcounter++;
}
//make new elements and add to list from list of nodes and connectivity from triangle
......@@ -270,8 +270,7 @@ void FaceMesh::AddNewPoint(Array<OneD, NekDouble> uv)
Array<OneD, NekDouble> np = m_cadsurf->P(uv);
NekDouble npDelta = m_octree->Query(np);
NodeSharedPtr n = boost::shared_ptr<Node>(
new Node(m_mesh->m_meshnode.size(),np[0],np[1],np[2]));
NodeSharedPtr n = boost::shared_ptr<Node>(new Node(0,np[0],np[1],np[2]));
bool add = true;
......@@ -306,7 +305,6 @@ void FaceMesh::AddNewPoint(Array<OneD, NekDouble> uv)
if(add)
{
n->SetCADSurf(m_id,uv);
m_mesh->m_meshnode.push_back(n);
m_stienerpoints.push_back(n);
}
}
......
......@@ -79,7 +79,7 @@ Array<OneD, NekDouble> SurfaceMesh::EdgeGrad(NekDouble ux, NekDouble vx,
NekDouble dfmag = sqrt(df[0]*df[0] + df[1]*df[1]);
df[0] = df[0]/dfmag; df[1] = df[1]/dfmag;
if(dfmag < 1E-30)
if(dfmag < 1E-15)
{
valid = false;
}
......
......@@ -34,6 +34,7 @@
////////////////////////////////////////////////////////////////////////////////
#include <MeshUtils/SurfaceMeshing/SurfaceMesh.h>
#include <limits>
using namespace std;
namespace Nektar
......@@ -146,104 +147,109 @@ void SurfaceMesh::Find1DBounds(NekDouble &a, NekDouble &b,
Array<OneD, NekDouble> df,
Array<OneD, NekDouble> bounds)
{
bool aset = false; bool bset = false;
NekDouble K, L, inter;
a = -1.0*numeric_limits<double>::max();
b = numeric_limits<double>::max();
bool aset = false, bset = false;
NekDouble K;
//want a to be negative, against the gradient, but properly bounded!!
//check edges of bounding box one by one for intersect, because paralell
//lines some cases can be ingnored
//line 1 left edge;
if(!(fabs(df[0]) < 1E-30)) //wouldnt exist on this edge
if(!(fabs(df[0]) < 1E-15)) //wouldnt exist on this edge
{
K = (bounds[0] - uvi[0]) / df[0];
L = df[1] * K + uvi[1] - bounds[2];
inter = bounds[2] + L;
if(!(inter < bounds[2] || inter > bounds[3]))
if(K<0)
{
//hit
if(K < 0)
//talking about a
if(fabs(K) < fabs(a))
{
ASSERTL0(!aset, "parameter previously set");
a = K;
aset = true;
}
else
}
else
{
if(fabs(K) < fabs(b))
{
ASSERTL0(!bset, "parameter previously set")
b = K;
bset = true;
}
}
}
//line 2 right edge;
if(!(fabs(df[0]) < 1E-30)) //wouldnt exist on this edge
{
//line 2 right edge;
K = (bounds[1] - uvi[0]) / df[0];
L = df[1] * K + uvi[1] - bounds[2];
inter = bounds[2] + L;
if(!(inter < bounds[2] || inter > bounds[3]))
if(K<0)
{
//hit
if(K < 0)