Commit 45dd4bd4 authored by Michael Turner's avatar Michael Turner
Browse files

fix to tetrahderon problem

parent 6cbcdb12
......@@ -133,6 +133,7 @@ vector<Array<OneD, int> > TetGenInterface::Extract()
tet[3] = output.tetrahedronlist[i * 4 + 3];
tets.push_back(tet);
}
return tets;
}
......
......@@ -42,6 +42,7 @@
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <NekMeshUtils/MeshElements/Node.h>
#include <NekMeshUtils/MeshElements/Mesh.h>
#define TETLIBRARY
#include <tetgen.h>
......
......@@ -57,11 +57,6 @@ public:
{
ElementSharedPtr e = boost::shared_ptr<Element>(
new Tetrahedron(pConf, pNodeList, pTagList));
std::vector<FaceSharedPtr> faces = e->GetFaceList();
for (int i = 0; i < faces.size(); ++i)
{
faces[i]->m_elLink.push_back(std::pair<ElementSharedPtr, int>(e, i));
}
return e;
}
/// Element type
......
......@@ -34,7 +34,6 @@
////////////////////////////////////////////////////////////////////////////////
#include <NekMeshUtils/TetMeshing/TetMesh.h>
#include <NekMeshUtils/ExtLibInterface/TetGenInterface.h>
using namespace std;
namespace Nektar
......@@ -47,24 +46,27 @@ void TetMesh::Mesh()
if (m_mesh->m_verbose)
cout << endl << endl << "Tetrahdral mesh generation" << endl;
TetGenInterfaceSharedPtr tetgen =
MemoryManager<TetGenInterface>::AllocateSharedPtr();
tetgen = MemoryManager<TetGenInterface>::AllocateSharedPtr();
map<int, NodeSharedPtr> TetgenIdToNode;
map<int, int> NodeIdToTetgenId;
map<int, NodeSharedPtr> IdToNode;
// at this point all nodes are in m_mesh->m_vertexset, but if there is a
// boundary layer, we dont want all of them, also, tetgen ids must be
// sequential so there is a map from tetgen id to real nodes
// build sequentially ordered maps of nodes that exist and there delta value
// in the octree
map<int, NekDouble> TetgenIdToDelta;
vector<Array<OneD, int> >
surfacetris; // surface mesh connectivity based on tetgenids
map<int, NekDouble> IdToDelta;
vector<Array<OneD, int> > surfacetris;
int cnt = 0;
if(!m_usePSurface)
{
NodeSet::iterator it;
for(it = m_mesh->m_vertexSet.begin(); it != m_mesh->m_vertexSet.end(); it++)
{
IdToNode[(*it)->m_id] = *it;
IdToDelta[(*it)->m_id] = m_octree->Query((*it)->GetLoc());
}
// build surface mesh and node map from all surface elements
for (int i = 0; i < m_mesh->m_element[2].size(); i++)
{
......@@ -76,27 +78,15 @@ void TetMesh::Mesh()
Array<OneD, int> tri(3);
for (int j = 0; j < n.size(); j++)
{
map<int, int>::iterator it;
it = NodeIdToTetgenId.find(n[j]->m_id);
if (it == NodeIdToTetgenId.end())
{
tri[j] = cnt;
NodeIdToTetgenId[n[j]->m_id] = cnt;
TetgenIdToNode[cnt] = n[j];
TetgenIdToDelta[cnt] = m_octree->Query(n[j]->GetLoc());
cnt++;
}
else
{
tri[j] = it->second;
}
tri[j] = n[j]->m_id;
}
surfacetris.push_back(tri);
}
}
else
{
vector<unsigned int> blsurfs = m_blmesh->GetBLSurfs();
ASSERTL0(false,"logic needs replacing will not work currently");
/*vector<unsigned int> blsurfs = m_blmesh->GetBLSurfs();
vector<ElementSharedPtr> Psurf = m_blmesh->GetPsuedoSurf();
//surface triangles will need to be checked against surftopriface to get the right face
//all surface elements are sequentially numbered so this should be easy to find in map
......@@ -166,18 +156,18 @@ void TetMesh::Mesh()
}
}
surfacetris.push_back(tri);
}
}*/
}
if (m_mesh->m_verbose)
{
cout << "\tInital Node Count: " << TetgenIdToNode.size() << endl;
cout << "\tInital Node Count: " << IdToNode.size() << endl;
}
tetgen->InitialMesh(TetgenIdToNode, surfacetris);
tetgen->InitialMesh(IdToNode, surfacetris);
vector<Array<OneD, NekDouble> > newp;
int ctbefore = TetgenIdToNode.size();
int ctbefore = IdToNode.size();
int newpb;
do
......@@ -188,33 +178,45 @@ void TetMesh::Mesh()
for (int i = 0; i < newp.size(); i++)
{
NekDouble d = m_octree->Query(newp[i]);
TetgenIdToDelta[ctbefore + i] = d;
IdToDelta[ctbefore + i] = d;
}
tetgen->RefineMesh(TetgenIdToDelta);
tetgen->RefineMesh(IdToDelta);
} while (newpb != newp.size());
// make new map of all nodes to build tets.
newp.clear();
tetgen->GetNewPoints(ctbefore, newp);
for (int i = 0; i < newp.size(); i++)
{
NodeSharedPtr n = boost::shared_ptr<Node>(
new Node(m_mesh->m_numNodes++, newp[i][0], newp[i][1], newp[i][2]));
TetgenIdToNode[ctbefore + i] = n;
new Node(ctbefore + i, newp[i][0], newp[i][1], newp[i][2]));
IdToNode[ctbefore + i] = n;
}
m_tetconnect = tetgen->Extract();
/*m_mesh->m_vertexSet.clear(); m_mesh->m_faceSet.clear(); m_mesh->m_element[2].clear();
m_mesh->m_edgeSet.clear(); m_mesh->m_element[3].clear();
newp.clear();
tetgen->GetNewPoints(0,newp);
vector<NodeSharedPtr> nodes;
for (int i = 0; i < newp.size(); i++)
{
nodes.push_back(boost::shared_ptr<Node>(
new Node(i, newp[i][0], newp[i][1], newp[i][2])));
}*/
// tetgen->freetet();
// create tets
for (int i = 0; i < m_tetconnect.size(); i++)
{
vector<NodeSharedPtr> n;
n.push_back(TetgenIdToNode[m_tetconnect[i][0]]);
n.push_back(TetgenIdToNode[m_tetconnect[i][1]]);
n.push_back(TetgenIdToNode[m_tetconnect[i][2]]);
n.push_back(TetgenIdToNode[m_tetconnect[i][3]]);
n.push_back(IdToNode[m_tetconnect[i][0]]);
n.push_back(IdToNode[m_tetconnect[i][1]]);
n.push_back(IdToNode[m_tetconnect[i][2]]);
n.push_back(IdToNode[m_tetconnect[i][3]]);
ElmtConfig conf(LibUtilities::eTetrahedron, 1, false, false);
vector<int> tags;
tags.push_back(0);
......
......@@ -47,6 +47,8 @@
#include <NekMeshUtils/SurfaceMeshing/SurfaceMesh.h>
#include <NekMeshUtils/BLMeshing/BLMesh.h>
#include <NekMeshUtils/ExtLibInterface/TetGenInterface.h>
namespace Nektar
{
namespace NekMeshUtils
......@@ -97,6 +99,8 @@ private:
/// conncetivity of the tets from the interface
std::vector<Array<OneD, int> > m_tetconnect;
TetGenInterfaceSharedPtr tetgen;
};
typedef boost::shared_ptr<TetMesh> TetMeshSharedPtr;
......
......@@ -326,6 +326,7 @@ void InputCAD::Process()
{
m_mesh->m_expDim = 3;
m_tet = MemoryManager<TetMesh>::AllocateSharedPtr(m_mesh, m_octree);
}
m_tet->Mesh();
......@@ -339,6 +340,24 @@ void InputCAD::Process()
ProcessElements();
ProcessComposites();
FaceSet::iterator fit;
count = 0;
for(fit = m_mesh->m_faceSet.begin(); fit != m_mesh->m_faceSet.end(); fit++)
{
if((*fit)->m_elLink.size() != 2)
{
count++;
}
}
if (count - m_mesh->m_element[2].size() > 0)
{
cerr << "Error: mesh contains unconnected faces and is not valid "
<< count - m_mesh->m_element[2].size()
<< endl;
abort();
}
//return;
m_surfacemesh->HOSurf();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment