Commit 666103d6 authored by Michael Turner's avatar Michael Turner

clean up mesh elements, this code depends on MR 697:

parent c83d04c9
SET(NEKMESHUTILS_SOURCES
MeshElements/Mesh.cpp
MeshElements/Element.cpp
MeshElements/Edge.cpp
MeshElements/Face.cpp
MeshElements/BooleanOperators.cpp
MeshElements/Point.cpp
MeshElements/Line.cpp
......@@ -90,8 +92,8 @@ IF(NEKTAR_USE_MESHGEN)
TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PRIVATE ${TETGEN_LIBRARY})
TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PUBLIC ${OCC_LIBRARIES})
SET(OCC_DEF LIN LININTEL HAVE_WOK_CONFIG_H HAVE_CONFIG_H CSFDB)
SET_PROPERTY(TARGET NekMeshUtils APPEND PROPERTY COMPILE_DEFINITIONS ${OCC_DEF})
#SET(OCC_DEF LIN LININTEL HAVE_WOK_CONFIG_H HAVE_CONFIG_H CSFDB)
#SET_PROPERTY(TARGET NekMeshUtils APPEND PROPERTY COMPILE_DEFINITIONS ${OCC_DEF})
ADD_DEPENDENCIES(NekMeshUtils opencascade-6.9 tetgen-1.5)
ENDIF()
......
////////////////////////////////////////////////////////////////////////////////
//
// File: Edge.h
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Mesh manipulation objects.
//
////////////////////////////////////////////////////////////////////////////////
#include <NekMeshUtils/MeshElements/Edge.h>
#include <NekMeshUtils/CADSystem/CADCurve.h>
#include <NekMeshUtils/CADSystem/CADSurf.h>
namespace Nektar
{
namespace NekMeshUtils
{
using namespace std;
string Edge::GetXmlCurveString()
{
std::vector<NodeSharedPtr> nodeList;
GetCurvedNodes(nodeList);
std::stringstream s;
std::string str;
// put them into a string and return
for (int k = 0; k < nodeList.size(); ++k)
{
s << std::scientific << std::setprecision(8) << " "
<< nodeList[k]->m_x << " " << nodeList[k]->m_y << " "
<< nodeList[k]->m_z << " ";
}
return s.str();
}
SpatialDomains::SegGeomSharedPtr Edge::GetGeom(int coordDim)
{
static boost::mutex io_mutex;
boost::mutex::scoped_lock lock(io_mutex);
// Create edge vertices.
SpatialDomains::PointGeomSharedPtr p[2];
SpatialDomains::SegGeomSharedPtr ret;
p[0] = m_n1->GetGeom(coordDim);
p[1] = m_n2->GetGeom(coordDim);
// Create a curve if high-order information exists.
if (m_edgeNodes.size() > 0)
{
SpatialDomains::CurveSharedPtr c =
MemoryManager<SpatialDomains::Curve>::AllocateSharedPtr(
m_id, m_curveType);
c->m_points.push_back(p[0]);
for (int i = 0; i < m_edgeNodes.size(); ++i)
{
c->m_points.push_back(m_edgeNodes[i]->GetGeom(coordDim));
}
c->m_points.push_back(p[1]);
ret = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
m_id, coordDim, p, c);
}
else
{
ret = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
m_id, coordDim, p);
}
return ret;
}
void Edge::MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType edgeType,
int coordDim,
int &id)
{
int nPoints = order + 1;
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
Array<OneD, NekDouble> edgePoints;
LibUtilities::PointsKey edgeKey(nPoints, edgeType);
LibUtilities::PointsManager()[edgeKey]->GetPoints(edgePoints);
Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
for (int i = 0; i < coordDim; ++i)
{
phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
}
m_edgeNodes.resize(nPoints - 2);
for (int i = 1; i < nPoints - 1; ++i)
{
Array<OneD, NekDouble> x(3, 0.0);
for (int j = 0; j < coordDim; ++j)
{
x[j] = xmap->PhysEvaluate(edgePoints + i, phys[j]);
}
m_edgeNodes[i-1] = boost::shared_ptr<Node>(
new Node(id++, x[0], x[1], x[2]));
}
m_curveType = edgeType;
#ifdef NEKTAR_USE_MESHGEN
if(onSurf || onCurve)
{
if(onCurve)
{
for(int i = 0; i < m_edgeNodes.size(); i++)
{
Array<OneD, NekDouble> loc(3);
loc[0] = m_edgeNodes[i]->m_x;
loc[1] = m_edgeNodes[i]->m_y;
loc[2] = m_edgeNodes[i]->m_z;
NekDouble t = CADCurve->loct(loc);
m_edgeNodes[i]->SetCADCurve(CADCurveId,CADCurve,t);
loc = CADCurve->P(t);
m_edgeNodes[i]->m_x = loc[0];
m_edgeNodes[i]->m_y = loc[1];
m_edgeNodes[i]->m_z = loc[2];
std::vector<CADSurfSharedPtr> s = CADCurve->GetAdjSurf();
for(int j = 0; j < s.size(); j++)
{
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
m_edgeNodes[i]->SetCADSurf(s[j]->GetId(),s[j],uv);
}
}
}
else
{
for(int i = 0; i < m_edgeNodes.size(); i++)
{
Array<OneD, NekDouble> loc(3);
loc[0] = m_edgeNodes[i]->m_x;
loc[1] = m_edgeNodes[i]->m_y;
loc[2] = m_edgeNodes[i]->m_z;
Array<OneD, NekDouble> uv = CADSurf->locuv(loc);
loc = CADSurf->P(uv);
m_edgeNodes[i]->m_x = loc[0];
m_edgeNodes[i]->m_y = loc[1];
m_edgeNodes[i]->m_z = loc[2];
m_edgeNodes[i]->SetCADSurf(CADSurfId,CADSurf,uv);
}
}
}
#endif
}
}
}
......@@ -42,7 +42,6 @@
#include <SpatialDomains/SegGeom.h>
#include <NekMeshUtils/NekMeshUtilsDeclspec.h>
#include <NekMeshUtils/MeshElements/Element.h>
#include <NekMeshUtils/MeshElements/Node.h>
namespace Nektar
......@@ -72,6 +71,7 @@ public:
{
#ifdef NEKTAR_USE_MESHGEN
onCurve = false;
onSurf = false;
#endif
}
......@@ -81,6 +81,7 @@ public:
{
#ifdef NEKTAR_USE_MESHGEN
onCurve = false;
onSurf = false;
#endif
}
......@@ -114,100 +115,17 @@ public:
/// Creates a Nektar++ string listing the coordinates of all the
/// nodes.
NEKMESHUTILS_EXPORT std::string GetXmlCurveString() const
{
std::vector<NodeSharedPtr> nodeList;
GetCurvedNodes(nodeList);
std::stringstream s;
std::string str;
// put them into a string and return
for (int k = 0; k < nodeList.size(); ++k)
{
s << std::scientific << std::setprecision(8) << " "
<< nodeList[k]->m_x << " " << nodeList[k]->m_y << " "
<< nodeList[k]->m_z << " ";
}
return s.str();
}
NEKMESHUTILS_EXPORT std::string GetXmlCurveString();
/// Generate a SpatialDomains::SegGeom object for this edge.
NEKMESHUTILS_EXPORT SpatialDomains::SegGeomSharedPtr GetGeom(int coordDim)
{
// Create edge vertices.
SpatialDomains::PointGeomSharedPtr p[2];
SpatialDomains::SegGeomSharedPtr ret;
p[0] = m_n1->GetGeom(coordDim);
p[1] = m_n2->GetGeom(coordDim);
// Create a curve if high-order information exists.
if (m_edgeNodes.size() > 0)
{
SpatialDomains::CurveSharedPtr c =
MemoryManager<SpatialDomains::Curve>::AllocateSharedPtr(
m_id, m_curveType);
c->m_points.push_back(p[0]);
for (int i = 0; i < m_edgeNodes.size(); ++i)
{
c->m_points.push_back(m_edgeNodes[i]->GetGeom(coordDim));
}
c->m_points.push_back(p[1]);
ret = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
m_id, coordDim, p, c);
}
else
{
ret = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
m_id, coordDim, p);
}
return ret;
}
NEKMESHUTILS_EXPORT SpatialDomains::SegGeomSharedPtr GetGeom(int coordDim);
/// Make this edge an order @p order edge. @see Element::MakeOrder.
void MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType edgeType,
int coordDim,
int &id)
{
int nPoints = order + 1;
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
Array<OneD, NekDouble> edgePoints;
LibUtilities::PointsKey edgeKey(nPoints, edgeType);
LibUtilities::PointsManager()[edgeKey]->GetPoints(edgePoints);
Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
for (int i = 0; i < coordDim; ++i)
{
phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
}
m_edgeNodes.resize(nPoints - 2);
for (int i = 1; i < nPoints - 1; ++i)
{
Array<OneD, NekDouble> x(3, 0.0);
for (int j = 0; j < coordDim; ++j)
{
x[j] = xmap->PhysEvaluate(edgePoints + i, phys[j]);
}
m_edgeNodes[i-1] = boost::shared_ptr<Node>(
new Node(id++, x[0], x[1], x[2]));
}
m_curveType = edgeType;
}
int &id);
/// ID of edge.
unsigned int m_id;
......@@ -224,9 +142,11 @@ public:
#ifdef NEKTAR_USE_MESHGEN
bool onCurve;
/// id of cad curve which edge lies on
int CADCurveId;
CADCurveSharedPtr CADCurve;
bool onSurf;
int CADSurfId;
CADSurfSharedPtr CADSurf;
#endif
private:
......
......@@ -49,6 +49,142 @@ namespace Nektar
{
namespace NekMeshUtils
{
/**
* @brief A lightweight struct for dealing with high-order triangle
* alignment.
*
* The logic underlying these routines is taken from the original Nektar
* code.
*/
template <typename T> struct HOTriangle
{
HOTriangle(std::vector<int> pVertId, std::vector<T> pSurfVerts)
: vertId(pVertId), surfVerts(pSurfVerts)
{
}
HOTriangle(std::vector<int> pVertId) : vertId(pVertId)
{
}
/// The triangle vertex IDs
std::vector<int> vertId;
/// The triangle surface vertices -- templated so that this can
/// either be nodes or IDs.
std::vector<T> surfVerts;
/**
* @brief Rotates the triangle of data points inside #surfVerts
* counter-clockwise nrot times.
*
* @param nrot Number of times to rotate triangle.
*/
void Rotate(int nrot)
{
int n, i, j, cnt;
int np = ((int)sqrt(8.0 * surfVerts.size() + 1.0) - 1) / 2;
std::vector<T> tmp(np * np);
for (n = 0; n < nrot; ++n)
{
for (cnt = i = 0; i < np; ++i)
{
for (j = 0; j < np - i; ++j, cnt++)
{
tmp[i * np + j] = surfVerts[cnt];
}
}
for (cnt = i = 0; i < np; ++i)
{
for (j = 0; j < np - i; ++j, cnt++)
{
surfVerts[cnt] = tmp[(np - 1 - i - j) * np + i];
}
}
}
}
/**
* @brief Reflect data points inside #surfVerts.
*
* This applies a mapping essentially doing the following
* reordering:
*
* 9 9
* 7 8 -> 8 7
* 4 5 6 6 5 4
* 0 1 2 3 3 2 1 0
*/
void Reflect()
{
int i, j, cnt;
int np = ((int)sqrt(8.0 * surfVerts.size() + 1.0) - 1) / 2;
std::vector<T> tmp(np * np);
for (cnt = i = 0; i < np; ++i)
{
for (j = 0; j < np - i; ++j, cnt++)
{
tmp[i * np + np - i - 1 - j] = surfVerts[cnt];
}
}
for (cnt = i = 0; i < np; ++i)
{
for (j = 0; j < np - i; ++j, cnt++)
{
surfVerts[cnt] = tmp[i * np + j];
}
}
}
/**
* @brief Align this surface to a given vertex ID.
*/
void Align(std::vector<int> vertId)
{
if (vertId[0] == this->vertId[0])
{
if (vertId[1] == this->vertId[1] || vertId[1] == this->vertId[2])
{
if (vertId[1] == this->vertId[2])
{
Rotate(1);
Reflect();
}
}
}
else if (vertId[0] == this->vertId[1])
{
if (vertId[1] == this->vertId[0] || vertId[1] == this->vertId[2])
{
if (vertId[1] == this->vertId[0])
{
Reflect();
}
else
{
Rotate(2);
}
}
}
else if (vertId[0] == this->vertId[2])
{
if (vertId[1] == this->vertId[0] || vertId[1] == this->vertId[1])
{
if (vertId[1] == this->vertId[1])
{
Rotate(2);
Reflect();
}
else
{
Rotate(1);
}
}
}
}
};
/**
* @brief Basic information about an element.
*
......@@ -83,19 +219,18 @@ struct ElmtConfig
/// Element type (e.g. triangle, quad, etc).
LibUtilities::ShapeType m_e;
/// Denotes whether the element contains face nodes. For 2D
/// elements, if this is true then the element contains interior
/// nodes.
/// Denotes whether the element contains face nodes. For 2D elements, if
/// this is true then the element contains interior nodes.
bool m_faceNodes;
/// Denotes whether the element contains volume (i.e. interior)
/// nodes. These are not supported by either the mesh converter or
/// Nektar++ but are included for completeness and are required
/// for some output modules (e.g. Gmsh).
/// Denotes whether the element contains volume (i.e. interior) nodes. These
/// are not supported by either the mesh converter or Nektar++ but are
/// included for completeness and are required for some output modules
/// (e.g. Gmsh).
bool m_volumeNodes;
/// Order of the element.
unsigned int m_order;
/// Denotes whether the element needs to be re-orientated for a
/// spectral element framework.
/// Denotes whether the element needs to be re-orientated for a spectral
/// element framework.
bool m_reorient;
/// Distribution of points in edges.
LibUtilities::PointsType m_edgeCurveType;
......@@ -138,6 +273,11 @@ public:
{
return m_conf;
}
///returns the shapetype
NEKMESHUTILS_EXPORT LibUtilities::ShapeType GetShapeType() const
{
return m_conf.m_e;
}
/// Returns the tag which defines the element shape.
NEKMESHUTILS_EXPORT std::string GetTag() const
{
......@@ -213,6 +353,15 @@ public:
}
else
{
for (int i = 0; i < m_face.size(); ++i)
{
n += m_face[i]->GetNodeCount();
}
for (int i = 0; i < m_edge.size(); ++i)
{
n -= m_edge[i]->GetNodeCount();
}
n += m_vertex.size();
std::cerr << "Not supported." << std::endl;
exit(1);
}
......@@ -325,9 +474,13 @@ public:
return s.str();
}
NEKMESHUTILS_EXPORT void GetCurvedNodes(
NEKMESHUTILS_EXPORT virtual void GetCurvedNodes(
std::vector<NodeSharedPtr> &nodeList) const
{
ASSERTL0(false,
"This function should be implemented at a shape level.");
/*std::cerr << "WARNING: Unsupported curvature for a " << m_vertex.size()
<< "-vertex element is not yet implemented." << std::endl;
// Node orderings are different for different elements.
// Triangle
if (m_vertex.size() == 2)
......@@ -417,7 +570,7 @@ public:
std::cerr << "GetXmlCurveString for a " << m_vertex.size()
<< "-vertex element is not yet implemented."
<< std::endl;
}
}*/
}
/// Generates a string listing the coordinates of all nodes
......
////////////////////////////////////////////////////////////////////////////////
//
// File: Face.h
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER