Commit dd07e528 authored by Michael Turner's avatar Michael Turner

tet meshing working

parent c41a3b47
......@@ -7,13 +7,13 @@ IF(NEKTAR_USE_MESH)
Octree/Octant.cpp
Octree/Octree.cpp
SurfaceMeshing/SurfaceMesh.cpp
# SurfaceMeshing/SurfaceMeshHOMesh.cpp
SurfaceMeshing/SurfaceMeshHOMesh.cpp
SurfaceMeshing/FaceMesh.cpp
# SurfaceMeshing/OptimiseFunctions.cpp
# SurfaceMeshing/OptimiseUtils.cpp
SurfaceMeshing/OptimiseFunctions.cpp
SurfaceMeshing/OptimiseUtils.cpp
ExtLibInterface/TriangleInterface.cpp
# ExtLibInterface/TetGenInterface.cpp
# TetMeshing/TetMesh.cpp
ExtLibInterface/TetGenInterface.cpp
TetMeshing/TetMesh.cpp
)
ENDIF()
......@@ -44,8 +44,8 @@ IF(NEKTAR_USE_MESH)
SurfaceMeshing/SurfaceMesh.h
SurfaceMeshing/FaceMesh.h
ExtLibInterface/TriangleInterface.h
# ExtLibInterface/TetGenInterface.h
# TetMeshing/TetMesh.h
ExtLibInterface/TetGenInterface.h
TetMeshing/TetMesh.h
)
ENDIF()
......
......@@ -39,16 +39,11 @@ using namespace std;
namespace Nektar{
namespace MeshUtils{
void TetGenInterface::InitialMesh(const std::vector<int> &nis,
const std::vector<int> &nsti,
std::map<int, MeshTriSharedPtr> &Tris,
std::map<int, MeshNodeSharedPtr> &Nodes)
void TetGenInterface::InitialMesh(map<int, NodeSharedPtr> nis,
vector<ElementSharedPtr> tri)
{
surface.initialize();
output.initialize();
additional.initialize();
nodemap.clear(); nodemapr.clear();
//build surface input
tetgenio::facet *f;
......@@ -58,45 +53,20 @@ void TetGenInterface::InitialMesh(const std::vector<int> &nis,
surface.numberofpoints = nis.size();
surface.pointlist = new REAL[surface.numberofpoints*3];
additional.firstnumber = 0;
additional.numberofpoints = nsti.size();
additional.pointlist = new REAL[additional.numberofpoints*3];
int pointc = 0;
for(int i = 0; i < nis.size(); i++)
{
nodemap[pointc] = nis[i];
nodemapr[nis[i]] = pointc;
Array<OneD, NekDouble> loc = Nodes[nis[i]]->GetLoc();
Array<OneD, NekDouble> loc = nis[i]->GetLoc();
surface.pointlist[i*3+0] = loc[0];
surface.pointlist[i*3+1] = loc[1];
surface.pointlist[i*3+2] = loc[2];
pointc++;
}
for(int i = 0; i < nsti.size(); i++)
{
nodemap[pointc] = nsti[i];
nodemapr[nsti[i]] = pointc;
Array<OneD, NekDouble> loc = Nodes[nsti[i]]->GetLoc();
additional.pointlist[i*3+0] = loc[0];
additional.pointlist[i*3+1] = loc[1];
additional.pointlist[i*3+2] = loc[2];
pointc++;
}
surface.numberoffacets = Tris.size();
surface.numberoffacets = tri.size();
surface.facetlist = new tetgenio::facet[surface.numberoffacets];
surface.facetmarkerlist = new int[surface.numberoffacets];
for(int i = 0; i < Tris.size(); i++)
for(int i = 0; i < tri.size(); i++)
{
f = &surface.facetlist[i];
f->numberofpolygons = 1;
......@@ -107,15 +77,15 @@ void TetGenInterface::InitialMesh(const std::vector<int> &nis,
p->numberofvertices = 3;
p->vertexlist = new int[p->numberofvertices];
Array<OneD, int> n = Tris[i]->GetN();
p->vertexlist[0] = nodemapr[n[0]];
p->vertexlist[1] = nodemapr[n[1]];
p->vertexlist[2] = nodemapr[n[2]];
vector<NodeSharedPtr> n = tri[i]->GetVertexList();
p->vertexlist[0] = n[0]->m_id;
p->vertexlist[1] = n[1]->m_id;
p->vertexlist[2] = n[2]->m_id;
surface.facetmarkerlist[i] = 0;
}
tetrahedralize("pYizqQ", &surface, &output, &additional);
tetrahedralize("pYzqQ", &surface, &output);
cout << output.numberofpoints << " " << output.numberoftetrahedra << endl;
}
......@@ -131,11 +101,7 @@ void TetGenInterface::GetNewPoints(int num, vector<Array<OneD, NekDouble> > &new
}
}
void TetGenInterface::RefineMesh(int num,
const std::vector<NekDouble> &ndel,
std::map<int, MeshTriSharedPtr> &Tris,
std::map<int, MeshNodeSharedPtr> &Nodes,
const std::vector<NekDouble> &newndel)
void TetGenInterface::RefineMesh(std::map<int, NekDouble> delta)
{
input = output;
......@@ -143,48 +109,29 @@ void TetGenInterface::RefineMesh(int num,
input.pointmtrlist = new REAL[input.numberofpoints];
for(int i = 0; i < num; i++)
{
input.pointmtrlist[i] = ndel[i];
}
int c = 0;
for(int i = num; i < input.numberofpoints; i++)
for(int i = 0; i < input.numberofpoints; i++)
{
input.pointmtrlist[i] = newndel[c];
c++;
input.pointmtrlist[i] = delta[i];
}
tetrahedralize("pYrmzq1.41/15QO2/7o/120", &input, &output);
}
tetrahedralize("pYrmzqQO2/7o/120", &input, &output);
cout << output.numberofpoints << " " << output.numberoftetrahedra << endl;
void TetGenInterface::AddNodes(int num, map<int, MeshNodeSharedPtr> &Nodes)
{
for(int i = num; i < output.numberofpoints; i++)
{
MeshNodeSharedPtr n = MemoryManager<MeshNode>::AllocateSharedPtr(
Nodes.size(), output.pointlist[i*3+0], output.pointlist[i*3+1],
output.pointlist[i*3+2]);
nodemap[i] = Nodes.size();
Nodes[Nodes.size()] = n;
}
}
void TetGenInterface::Extract(int &numtet,
Array<OneD, Array<OneD, int> > &tetconnect)
vector<Array<OneD, int> > TetGenInterface::Extract()
{
numtet = output.numberoftetrahedra;
tetconnect = Array<OneD, Array<OneD, int> >(numtet);
for(int i = 0; i < numtet; i++)
vector<Array<OneD, int> > tets;
for(int i = 0; i < output.numberoftetrahedra; i++)
{
Array<OneD, int> tet(4);
tet[0] = nodemap[output.tetrahedronlist[i*4+0]];
tet[1] = nodemap[output.tetrahedronlist[i*4+1]];
tet[2] = nodemap[output.tetrahedronlist[i*4+2]];
tet[3] = nodemap[output.tetrahedronlist[i*4+3]];
tetconnect[i] = tet;
tet[0] = output.tetrahedronlist[i*4+0];
tet[1] = output.tetrahedronlist[i*4+1];
tet[2] = output.tetrahedronlist[i*4+2];
tet[3] = output.tetrahedronlist[i*4+3];
tets.push_back(tet);
}
return tets;
}
void TetGenInterface::freetet()
......
......@@ -44,10 +44,13 @@
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <MeshUtils/MeshElements/MeshElements.h>
#include <MeshUtils/MeshElem.hpp>
namespace Nektar{
namespace MeshUtils{
namespace Nektar
{
namespace MeshUtils
{
class TetGenInterface
{
......@@ -64,10 +67,8 @@ class TetGenInterface
/**
* @brief assign parameters for meshing
*/
void InitialMesh(const std::vector<int> &nis,
const std::vector<int> &nsti,
std::map<int, MeshTriSharedPtr> &Tris,
std::map<int, MeshNodeSharedPtr> &Nodes);
void InitialMesh(std::map<int, NodeSharedPtr> nis,
std::vector<ElementSharedPtr> tri);
/**
* @brief gets the locations of the stiener points added by tetgen
......@@ -77,20 +78,9 @@ class TetGenInterface
/**
* @brief refines a previously made tetmesh with node delta information from the Octree
*/
void RefineMesh(int num,
const std::vector<NekDouble> &ndel,
std::map<int, MeshTriSharedPtr> &Tris,
std::map<int, MeshNodeSharedPtr> &Nodes,
const std::vector<NekDouble> &newndel);
/**
* @brief extract mesh information
*/
void Extract(int &numtet, Array<OneD, Array<OneD, int> > &tetconnect);
void RefineMesh(std::map<int, NekDouble> delta);
/**
* @brief from the list of stiener points creates new mesh node objects
*/
void AddNodes(int num, std::map<int, MeshNodeSharedPtr> &Nodes);
std::vector<Array<OneD, int> > Extract();
/**
* @brief clear previous mesh
......@@ -101,10 +91,6 @@ class TetGenInterface
///tetgen objects
tetgenio surface, output, input, additional;
/// map from meshconvert id to tetgen id
std::map<int, int> nodemap;
/// map in reverse
std::map<int, int> nodemapr;
};
typedef boost::shared_ptr<TetGenInterface> TetGenInterfaceSharedPtr;
......
......@@ -138,6 +138,7 @@ namespace MeshUtils
vector<pair<ElementSharedPtr, int> > m_elLink;
/// id of cad curve which edge lies on
unsigned int CADCurveID;
unsigned int CADSurfID;
private:
SpatialDomains::SegGeomSharedPtr m_geom;
......
......@@ -33,13 +33,15 @@
//
////////////////////////////////////////////////////////////////////////////////
#include <MeshUtils/SurfaceMeshing/SurfaceMeshing.h>
#include <MeshUtils/SurfaceMeshing/SurfaceMesh.h>
using namespace std;
namespace Nektar{
namespace MeshUtils{
namespace Nektar
{
namespace MeshUtils
{
Array<OneD, NekDouble> SurfaceMeshing::EdgeGrad(NekDouble ux, NekDouble vx,
Array<OneD, NekDouble> SurfaceMesh::EdgeGrad(NekDouble ux, NekDouble vx,
vector<Array<OneD,NekDouble> > bcs,
vector<NekDouble> weights,
int surf, bool &valid)
......@@ -88,7 +90,7 @@ Array<OneD, NekDouble> SurfaceMeshing::EdgeGrad(NekDouble ux, NekDouble vx,
return df;
}
NekDouble SurfaceMeshing::EdgeF(NekDouble ux, NekDouble vx,
NekDouble SurfaceMesh::EdgeF(NekDouble ux, NekDouble vx,
vector<Array<OneD,NekDouble> > bcs,
vector<NekDouble> weights,
int surf, bool &valid)
......@@ -115,7 +117,7 @@ NekDouble SurfaceMeshing::EdgeF(NekDouble ux, NekDouble vx,
}
Array<OneD, NekDouble> SurfaceMeshing::FaceGrad(NekDouble ux, NekDouble vx,
Array<OneD, NekDouble> SurfaceMesh::FaceGrad(NekDouble ux, NekDouble vx,
vector<Array<OneD,NekDouble> > bcs,
vector<NekDouble> weights,
int surf, bool &valid)
......@@ -172,7 +174,7 @@ Array<OneD, NekDouble> SurfaceMeshing::FaceGrad(NekDouble ux, NekDouble vx,
return df;
}
NekDouble SurfaceMeshing::FaceF(NekDouble ux, NekDouble vx,
NekDouble SurfaceMesh::FaceF(NekDouble ux, NekDouble vx,
vector<Array<OneD,NekDouble> > bcs,
vector<NekDouble> weights,
int surf, bool &valid)
......
......@@ -33,13 +33,15 @@
//
////////////////////////////////////////////////////////////////////////////////
#include <MeshUtils/SurfaceMeshing/SurfaceMeshing.h>
#include <MeshUtils/SurfaceMeshing/SurfaceMesh.h>
using namespace std;
namespace Nektar{
namespace MeshUtils{
namespace Nektar
{
namespace MeshUtils
{
NekDouble SurfaceMeshing::BrentOpti(NekDouble ax, NekDouble bx,
NekDouble SurfaceMesh::BrentOpti(NekDouble ax, NekDouble bx,
NekDouble cx, NekDouble &fx,
NekDouble tol, int surf,
Array<OneD, NekDouble> uvi,
......@@ -47,7 +49,7 @@ NekDouble SurfaceMeshing::BrentOpti(NekDouble ax, NekDouble bx,
Array<OneD, NekDouble> bounds,
std::vector<Array<OneD,NekDouble> > bcs,
std::vector<NekDouble> weights,
NekDouble (SurfaceMeshing::*GetF)(
NekDouble (SurfaceMesh::*GetF)(
NekDouble, NekDouble,
std::vector<Array<OneD,NekDouble> >,
std::vector<NekDouble>, int, bool &))
......@@ -139,7 +141,7 @@ NekDouble SurfaceMeshing::BrentOpti(NekDouble ax, NekDouble bx,
return xmin;
}
void SurfaceMeshing::Find1DBounds(NekDouble &a, NekDouble &b,
void SurfaceMesh::Find1DBounds(NekDouble &a, NekDouble &b,
Array<OneD, NekDouble> uvi,
Array<OneD, NekDouble> df,
Array<OneD, NekDouble> bounds)
......
......@@ -34,6 +34,7 @@
////////////////////////////////////////////////////////////////////////////////
#include <list>
#include <algorithm>
#include <MeshUtils/SurfaceMeshing/SurfaceMesh.h>
#include <LibUtilities/BasicUtils/Progressbar.hpp>
......@@ -51,99 +52,149 @@ namespace MeshUtils
* needed because OCC curves are not distorted, but should be included for
* completness
*/
void SurfaceMeshing::HOSurf()
void SurfaceMesh::HOSurf()
{
if(m_verbose)
if(m_mesh->m_verbose)
cout << endl << "\tHigh-Order Surface meshing" << endl;
map<int, MeshEdgeSharedPtr>::iterator eit;
EdgeSet::iterator eit;
int counter = 0;
LibUtilities::PointsKey ekey(m_order+1,
LibUtilities::PointsKey ekey(m_mesh->m_nummode,
LibUtilities::eGaussLobattoLegendre);
Array<OneD, NekDouble> gll;
LibUtilities::PointsManager()[ekey]->GetPoints(gll);
//loop over all edges in the surface
for(eit = Edges.begin(); eit != Edges.end(); eit++)
for(eit = m_mesh->m_edgeSet.begin(); eit != m_mesh->m_edgeSet.end(); eit++)
{
if(m_verbose)
if(m_mesh->m_verbose)
{
LibUtilities::PrintProgressbar(counter,Edges.size(),
LibUtilities::PrintProgressbar(counter,m_mesh->m_edgeSet.size(),
"\t\tEdges");
}
counter++;
MeshEdgeSharedPtr e = eit->second;
Array<OneD, int> n = e->GetN();
if((*eit)->m_n1->GetListCADCurve().size() > 0 && (*eit)->m_n2->GetListCADCurve().size() > 0)
{
vector<int> clist1, clist2, c;
clist1 = (*eit)->m_n1->GetListCADCurve();
clist2 = (*eit)->m_n2->GetListCADCurve();
sort(clist1.begin(), clist1.end());
sort(clist2.begin(), clist2.end());
set_intersection(clist1.begin(), clist1.end(),
clist2.begin(), clist2.end(),
back_inserter(c));
if(c.size() > 0)
{
(*eit)->CADCurveID = c[0];
}
else
{
vector<int> slist1, slist2, s;
slist1 = (*eit)->m_n1->GetListCADSurf();
slist2 = (*eit)->m_n2->GetListCADSurf();
sort(slist1.begin(), slist1.end());
sort(slist2.begin(), slist2.end());
set_intersection(slist1.begin(), slist1.end(),
slist2.begin(), slist2.end(),
back_inserter(s));
if(s.size()>0)
{
(*eit)->CADSurfID = s[0];
}
else
{
ASSERTL0(false,"failed to identify location of edge");
}
}
}
else
{
vector<int> slist1, slist2, s;
slist1 = (*eit)->m_n1->GetListCADSurf();
slist2 = (*eit)->m_n2->GetListCADSurf();
sort(slist1.begin(), slist1.end());
sort(slist2.begin(), slist2.end());
set_intersection(slist1.begin(), slist1.end(),
slist2.begin(), slist2.end(),
back_inserter(s));
if(s.size()>0)
{
(*eit)->CADSurfID = s[0];
}
else
{
ASSERTL0(false,"failed to identify location of edge");
}
}
//check if edge in on a surface or curve
if(e->GetCurve() != -1)
if((*eit)->CADCurveID != -1)
{
LibUtilities::CADCurveSharedPtr c = m_cad->GetCurve(
e->GetCurve());
int cid = (*eit)->CADCurveID;
CADCurveSharedPtr c = m_cad->GetCurve(cid);
//edge is on curve and needs hoing that way
NekDouble tb = Nodes[n[0]]->GetC(e->GetCurve());
NekDouble te = Nodes[n[1]]->GetC(e->GetCurve());
NekDouble tb = (*eit)->m_n1->GetCADCurve(cid);
NekDouble te = (*eit)->m_n2->GetCADCurve(cid);
Array<OneD, NekDouble> ti(m_order+1);
Array<OneD, NekDouble> ti(m_mesh->m_nummode);
for(int i = 0; i < m_order+1; i++)
for(int i = 0; i < m_mesh->m_nummode; i++)
{
ti[i] = tb*(1.0 - gll[i])/2.0 +
te*(1.0 + gll[i])/2.0;
}
vector<int> Surfs = c->GetAdjSurf();
ASSERTL0(Surfs.size() == 2, "Number of common surfs should be 2");
vector<CADSurfSharedPtr> s = c->GetAdjSurf();
vector<int> honodes(m_order-1);
ASSERTL0(s.size() == 2, "Number of common surfs should be 2");
LibUtilities::CADSurfSharedPtr s1,s2;
s1 = m_cad->GetSurf(Surfs[0]);
s2 = m_cad->GetSurf(Surfs[1]);
vector<NodeSharedPtr> honodes(m_mesh->m_nummode-2);
for(int i = 1; i < m_order +1 -1; i++)
for(int i = 1; i < m_mesh->m_nummode -1; i++)
{
Array<OneD, NekDouble> loc = c->P(ti[i]);
MeshNodeSharedPtr nn = MemoryManager<MeshNode>::
AllocateSharedPtr(Nodes.size(),loc[0],
loc[1],loc[2]);
nn->SetCurve(c->GetID(),ti[i]);
Array<OneD, NekDouble> uv = s1->locuv(loc);
nn->SetSurf(Surfs[0],uv);
uv = s2->locuv(loc);
nn->SetSurf(Surfs[1],uv);
honodes[i-1] = Nodes.size();
Nodes[Nodes.size()] = nn;
NodeSharedPtr nn = MemoryManager<Node>::
AllocateSharedPtr(-1,loc[0],loc[1],loc[2]);
nn->SetCADCurve(cid,ti[i]);
Array<OneD, NekDouble> uv = s[0]->locuv(loc);
nn->SetCADSurf(s[0]->GetId(),uv);
uv = s[1]->locuv(loc);
nn->SetCADSurf(s[1]->GetId(),uv);
honodes[i-1] = nn;
}
e->SetHONodes(honodes);
(*eit)->m_edgeNodes = honodes;
(*eit)->m_curveType = LibUtilities::eGaussLobattoLegendre;
}
else
{
//edge is on surface and needs 2d optimisation
LibUtilities::CADSurfSharedPtr s = m_cad->GetSurf(e->GetSurf());
int sid = (*eit)->CADSurfID;
CADSurfSharedPtr s = m_cad->GetSurf(sid);
Array<OneD, NekDouble> uvb,uve;
uvb = Nodes[n[0]]->GetS(e->GetSurf());
uve = Nodes[n[1]]->GetS(e->GetSurf());
uvb = (*eit)->m_n1->GetCADSurf(sid);
uve = (*eit)->m_n1->GetCADSurf(sid);
vector<int> honodes(m_order-1);
vector<NekDouble> gllint(m_order-1);
for(int i = 1; i < m_order+1 -1; i++)
vector<NodeSharedPtr> honodes(m_mesh->m_nummode-2);
vector<NekDouble> gllint(m_mesh->m_nummode-2);
for(int i = 1; i < m_mesh->m_nummode -1; i++)
{
Array<OneD, NekDouble> loc;
Array<OneD, NekDouble> uv(2);
uv[0] = uvb[0]*(1.0 - gll[i])/2.0 + uve[0]*(1.0 + gll[i])/2.0;
uv[1] = uvb[1]*(1.0 - gll[i])/2.0 + uve[1]*(1.0 + gll[i])/2.0;
loc = s->P(uv);
MeshNodeSharedPtr nn = MemoryManager<MeshNode>::
AllocateSharedPtr(Nodes.size(),loc[0],
loc[1],loc[2]);
nn->SetSurf(e->GetSurf(),uv);
honodes[i-1] = Nodes.size();
Nodes[Nodes.size()] = nn;
NodeSharedPtr nn = MemoryManager<Node>::
AllocateSharedPtr(-1,loc[0],loc[1],loc[2]);
nn->SetCADSurf(sid,uv);
honodes[i-1] = nn;
gllint[i-1] = gll[i];
}
......@@ -173,7 +224,7 @@ void SurfaceMeshing::HOSurf()
}
else
{
bcs.push_back(Nodes[honodes[i-1]]->GetS(e->GetSurf()));
bcs.push_back(honodes[i-1]->GetCADSurf(sid));
za = gllint[i-1];
}
if(i==honodes.size()-1)
......@@ -183,7 +234,7 @@ void SurfaceMeshing::HOSurf()
}
else
{
bcs.push_back(Nodes[honodes[i+1]]->GetS(e->GetSurf()));
bcs.push_back(honodes[i+1]->GetCADSurf(sid));
zb = gllint[i+1];
}
......@@ -193,11 +244,11 @@ void SurfaceMeshing::HOSurf()
weights[0] = 1.0/(zb - zm);
weights[1] = 1.0 /(zm - za);
uvi = Nodes[honodes[i]]->GetS(e->GetSurf());
uvi = honodes[i]->GetCADSurf(sid);