Commit 2b18237b authored by Dave Moxey's avatar Dave Moxey
Browse files

Various updates to MeshConvert.

- Updating BLPoints class to add a static variable to supply spacing option.
- Adding automatic layer width determination into ProcessBL along with more
  configuration options.
- Added 3D spherigon support to ProcessSpherigon, plus support for reading
  vertex normals from InputNek and InputPly.
- Fixed reorientation issue with Nektar meshes where prisms read after tets.
- Added documentation for spherigons/configuration option classes.


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4028 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent a7024c88
......@@ -46,25 +46,23 @@ namespace Nektar
{
namespace LibUtilities
{
// Default value.
double BLPoints::delta_star = 1.0;
void BLPoints::CalculatePoints()
{
// Allocate the storage for points
// Allocate the storage for points.
PointsBaseType::CalculatePoints();
unsigned int npts = m_pointsKey.GetNumPoints();
double rr = pow(10.0,log10(4000.0)/(npts-2.0));
//cerr << "power coefficiant : " << rr << endl;
for(unsigned int i=0;i<npts;++i)
unsigned int npts = m_pointsKey.GetNumPoints();
// Derived power coefficient.
double rn = powf(2.0/BLPoints::delta_star,1.0/(npts-2.0));
m_points[0][0] = -1.0;
for (unsigned int i = 1; i < npts; ++i)
{
if ( i == 0 )
{
m_points[0][i] = -1.0;
//cerr << m_points[0][i] +1.0 << endl;
}
else
{
m_points[0][i] = -1 + (1.0/2000.0)*pow(rr,(double(i-1)));
//cerr << m_points[0][i] +1.0 << endl;
}
m_points[0][i] = -1.0 + delta_star*pow(rn,(double)i-1.0);
}
}
......
......@@ -49,67 +49,91 @@ namespace Nektar
{
class BLPoints: public Points<NekDouble>
{
public:
typedef Points<NekDouble> PointsBaseType;
virtual ~BLPoints()
{
}
public:
typedef Points<NekDouble> PointsBaseType;
virtual ~BLPoints()
{
}
LIB_UTILITIES_EXPORT static boost::shared_ptr< PointsBaseType > Create(const PointsKey &key);
LIB_UTILITIES_EXPORT boost::shared_ptr< NekMatrix<NekDouble> > CreateMatrix(const PointsKey &pkey);
LIB_UTILITIES_EXPORT static boost::shared_ptr< PointsBaseType >
Create(const PointsKey &key);
LIB_UTILITIES_EXPORT boost::shared_ptr< NekMatrix<NekDouble> >
CreateMatrix(const PointsKey &pkey);
LIB_UTILITIES_EXPORT const MatrixSharedPtrType GetI(const PointsKey &pkey);
LIB_UTILITIES_EXPORT const MatrixSharedPtrType GetI(const Array<OneD, const NekDouble>& x);
LIB_UTILITIES_EXPORT const MatrixSharedPtrType GetI(unsigned int numpoints, const Array<OneD, const NekDouble>& x);
LIB_UTILITIES_EXPORT const MatrixSharedPtrType
GetI(const PointsKey &pkey);
LIB_UTILITIES_EXPORT const MatrixSharedPtrType
GetI(const Array<OneD, const NekDouble>& x);
LIB_UTILITIES_EXPORT const MatrixSharedPtrType
GetI(unsigned int numpoints,
const Array<OneD, const NekDouble>& x);
BLPoints(const PointsKey &key):PointsBaseType(key)
{
m_InterpManager.RegisterCreator(PointsKey(0, eGaussGaussLegendre),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussRadauMLegendre),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussRadauPLegendre),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussLobattoLegendre),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussGaussChebyshev),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussRadauMChebyshev),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussRadauPChebyshev),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussLobattoChebyshev),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha0Beta1),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha0Beta2),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha1Beta0),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eGaussRadauMAlpha2Beta0),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, ePolyEvenlySpaced),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eFourierEvenlySpaced),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(PointsKey(0, eBoundaryLayerPoints),
boost::bind(&BLPoints::CreateMatrix, this, _1));
}
BLPoints(const PointsKey &key):PointsBaseType(key)
{
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussGaussLegendre),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussRadauMLegendre),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussRadauPLegendre),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussLobattoLegendre),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussGaussChebyshev),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussRadauMChebyshev),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussRadauPChebyshev),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussLobattoChebyshev),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussRadauMAlpha0Beta1),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussRadauMAlpha0Beta2),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussRadauMAlpha1Beta0),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eGaussRadauMAlpha2Beta0),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, ePolyEvenlySpaced),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eFourierEvenlySpaced),
boost::bind(&BLPoints::CreateMatrix, this, _1));
m_InterpManager.RegisterCreator(
PointsKey(0, eBoundaryLayerPoints),
boost::bind(&BLPoints::CreateMatrix, this, _1));
}
private:
/// Default constructor should not be called except by Create method.
BLPoints();
static double delta_star;
private:
/// Default constructor should not be called except by Create
/// method.
BLPoints();
/// Copy constructor should not be called.
BLPoints(const BLPoints &points);
/// Copy constructor should not be called.
BLPoints(const BLPoints &points);
void CalculatePoints();
void CalculateWeights();
void CalculateDerivMatrix();
void CalculateInterpMatrix(
unsigned int npts,
const Array<OneD, const NekDouble>& xpoints,
Array<OneD, NekDouble>& interp);
void CalculatePoints();
void CalculateWeights();
void CalculateDerivMatrix();
void CalculateInterpMatrix(
unsigned int npts,
const Array<OneD, const NekDouble>& xpoints,
Array<OneD, NekDouble>& interp);
}; // class BLPoints
} // end of namespace
} // end of namespace
......
......@@ -72,6 +72,12 @@ namespace Nektar
typedef boost::shared_ptr<HOSurf> HOSurfSharedPtr;
enum NekCurve
{
eFile,
eRecon
};
/**
* Hash class for high-order surfaces.
*/
......@@ -122,7 +128,7 @@ namespace Nektar
/**
* Maps a curve tag to a filename containing surface information.
*/
std::map<string, string> curveTags;
std::map<string, pair<NekCurve, string> > curveTags;
/**
* Maps a curve tag to the high-order surface data for that tag.
......
......@@ -62,12 +62,6 @@ namespace Nektar
/**
* Gmsh file contains a list of nodes and their coordinates, along with
* a list of elements and those nodes which define them. We read in and
* store the list of nodes in #m_node and store the list of elements in
* #m_element. Each new element is supplied with a list of entries from
* #m_node which defines the element. Finally some mesh statistics are
* printed.
*
* @param pFilename Filename of Gmsh file to read.
*/
......@@ -82,9 +76,12 @@ namespace Nektar
int nEntities = 0;
int nElements = 0;
int nBoundaryElements = 0;
int nProperties = 0;
ElementType elType = eTriangle;
map<string, int> propMap;
cout << "Start reading InputPly..." << endl;
while (!mshFile.eof())
{
getline(mshFile, line);
......@@ -104,17 +101,29 @@ namespace Nektar
}
continue;
}
else if (word == "property")
{
s >> word >> word;
propMap[word] = nProperties++;
}
else if (word == "end_header")
{
// Read nodes
int id = 0;
vector<double> data(nProperties);
for (int i = 0; i < nVertices; ++i)
{
getline(mshFile, line);
stringstream st(line);
double x = 0, y = 0, z = 0;
st >> x >> y >> z;
for (int j = 0; j < nProperties; ++j)
{
st >> data[j];
}
double x = data[propMap["x"]];
double y = data[propMap["y"]];
double z = data[propMap["z"]];
if ((y * y) > 0.000001 && m->spaceDim != 3)
{
m->spaceDim = 2;
......@@ -123,13 +132,20 @@ namespace Nektar
{
m->spaceDim = 3;
}
id -= 1; // counter starts at 0
m->node.push_back(boost::shared_ptr<Node>(new Node(id, x, y, z)));
id++;
m->node.push_back(
boost::shared_ptr<Node>(new Node(i, x, y, z)));
// Read vertex normals.
if (propMap.count("nx") > 0)
{
double nx = data[propMap["nx"]];
double ny = data[propMap["ny"]];
double nz = data[propMap["nz"]];
m->vertexNormals[i] = Node(0, nx, ny, nz);
}
}
// Read elements
int zeroDid = 0, oneDid = 0, twoDid = 0, threeDid = 0;
for (int i = 0; i < nEntities; ++i)
{
getline(mshFile, line);
......@@ -139,8 +155,7 @@ namespace Nektar
// Create element tags
vector<int> tags;
tags.push_back(0); // composite
tags.push_back(eTriangle); // element type
// Read element node list
st >> id;
vector<NodeSharedPtr> nodeList;
......@@ -150,14 +165,15 @@ namespace Nektar
st >> node;
nodeList.push_back(m->node[node]);
}
// Create element
ElmtConfig conf(elType,1,false,false);
ElementSharedPtr E = GetElementFactory().
CreateInstance(elType,conf,nodeList,tags);
// Determine mesh expansion dimension
if (E->GetDim() > m->expDim) {
if (E->GetDim() > m->expDim)
{
m->expDim = E->GetDim();
}
m->element[E->GetDim()].push_back(E);
......
......@@ -222,13 +222,13 @@ int main(int argc, char* argv[])
if (tmp2.size() == 1)
{
mod->RegisterConfig(tmp2[0], "true");
mod->RegisterConfig(tmp2[0], "1");
}
else if (tmp2.size() == 2)
{
mod->RegisterConfig(tmp2[0], tmp2[1]);
}
if (tmp2.size() != 2)
else
{
cerr << "ERROR: Invalid module configuration: format is "
<< "either :arg or :arg=val" << endl;
......
......@@ -1127,7 +1127,6 @@ namespace Nektar
// point. Place second highest global vertex at base degenerate
// point.
sort(vertex.begin(), vertex.end());
//reverse(vertex.begin(), vertex.end());
// Calculate a.(b x c) to determine tet volume; if negative,
// reverse order of non-degenerate points to correctly orientate
......@@ -1153,7 +1152,7 @@ namespace Nektar
{
swap(vertex[0], vertex[1]);
}
TetOrientSet::iterator it;
// Search for the face in the original set of face nodes. Then use
......@@ -1247,13 +1246,13 @@ namespace Nektar
vector<NodeSharedPtr> faceVertices;
vector<EdgeSharedPtr> faceEdges;
vector<NodeSharedPtr> faceNodes;
for (int k = 0; k < 4; ++k)
int nEdge = 3 - (j % 2 - 1);
for (int k = 0; k < nEdge; ++k)
{
if (face_ids[j][k] == -1)
break;
faceVertices.push_back(vertex[face_ids[j][k]]);
NodeSharedPtr a = vertex[face_ids[j][k]];
NodeSharedPtr b = vertex[face_ids[j][(k+1)%(j%2==0?4:3)]];
NodeSharedPtr b = vertex[face_ids[j][(k+1) % nEdge]];
for (unsigned int i = 0; i < edge.size(); ++i)
{
if ((edge[i]->n1 == a && edge[i]->n2 == b) ||
......@@ -1279,7 +1278,7 @@ namespace Nektar
face.push_back(FaceSharedPtr(
new Face(faceVertices, faceNodes, faceEdges, m_conf.faceCurveType)));
}
// Re-order edge array to be consistent with Nektar++ ordering.
vector<EdgeSharedPtr> tmp(9);
tmp[0] = edge[face_edges[0][0]];
......@@ -1459,12 +1458,12 @@ namespace Nektar
* - 2 if the prism is rotated counter-clockwise.
*
* This is necessary for some input modules (e.g. #InputNek) which add
* high-order information to faces.
* high-order information or bounary conditions to faces.
*/
void Prism::OrientPrism()
{
int lid[6], gid[6];
// Re-order vertices.
for (int i = 0; i < 6; ++i)
{
......
......@@ -44,6 +44,7 @@
#include <boost/shared_ptr.hpp>
#include <boost/unordered_set.hpp>
#include <boost/unordered_map.hpp>
#include <LibUtilities/Foundations/Foundations.hpp>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
......@@ -96,7 +97,7 @@ namespace Nektar
Node(const Node& pSrc)
: id(pSrc.id), x(pSrc.x), y(pSrc.y),
z(pSrc.z), m_geom() {}
Node() {}
Node() : id(0), x(0.0), y(0.0), z(0.0), m_geom() {}
~Node() {}
/// Define node ordering based on ID.
......@@ -130,6 +131,11 @@ namespace Nektar
return Node(id, alpha*x, alpha*y, alpha*z);
}
Node operator/(const double &alpha) const
{
return Node(id, x/alpha, y/alpha, z/alpha);
}
void operator+=(const Node &pSrc)
{
x += pSrc.x;
......@@ -843,6 +849,8 @@ namespace Nektar
*/
class Composite {
public:
Composite() {}
/// Generate the list of IDs of elements within this composite.
std::string GetXmlString();
......@@ -918,6 +926,11 @@ namespace Nektar
ConditionMap condition;
/// List of fields names.
std::vector<std::string> fields;
/// Map of vertex normals.
boost::unordered_map<int, Node> vertexNormals;
/// Set of all pairs of element ID and face number on which to apply
/// spherigon surface smoothing.
set<pair<int,int> > spherigonFaces;
/// Returns the total number of elements in the mesh with
/// dimension expDim.
unsigned int GetNumElements();
......
......@@ -71,6 +71,9 @@ namespace Nektar
config["outfile"] = ConfigOption(false, "", "Output filename.");
}
/**
* @brief Open a file for input.
*/
void InputModule::OpenStream()
{
string fname = config["infile"].as<string>();
......@@ -82,6 +85,9 @@ namespace Nektar
}
}
/**
* @brief Open a file for output.
*/
void OutputModule::OpenStream()
{
string fname = config["outfile"].as<string>();
......@@ -94,6 +100,9 @@ namespace Nektar
}
/**
* @brief Create a unique set of mesh vertices from elements stored in
* Mesh::element.
*
* Each element is processed in turn and the vertices extracted and
* inserted into #m_vertexSet, which at the end of the routine
* contains all unique vertices in the mesh.
......@@ -123,8 +132,9 @@ namespace Nektar
}
/**
* This routine only proceeds if the expansion dimension is 2 or 3.
*
* @brief Create a unique set of mesh edges from elements stored in
* Mesh::element.
*
* All elements are first scanned and a list of unique, enumerated
* edges produced in #m_edgeSet. Since each element generated its
* edges independently, we must now ensure that each element only uses
......@@ -133,6 +143,8 @@ namespace Nektar
* 1-D boundary elements which correspond to an edge in
* #m_edgeSet. For such elements, we set its edgeLink to reference the
* corresponding edge in #m_edgeSet.
*
* This routine only proceeds if the expansion dimension is 2 or 3.
*/
void Module::ProcessEdges()
{
......@@ -191,8 +203,9 @@ namespace Nektar
/**
* This routine only proceeds if the expansion dimension is 3.
*
* @brief Create a unique set of mesh faces from elements stored in
* Mesh::element.
*
* All elements are scanned and a unique list of enumerated faces is
* produced in #m_faceSet. Since elements created their own faces
* independently, we examine each element only uses face objects from
......@@ -201,6 +214,8 @@ namespace Nektar
* elements for 2-D boundary faces which correspond to faces in
* #m_faceSet. For such elements, we set its faceLink to reference the
* corresponding face in #m_faceSet.
*
* This routine only proceeds if the expansion dimension is 3.
*/
void Module::ProcessFaces()
{
......@@ -260,8 +275,9 @@ namespace Nektar
}
}
/**
* @brief Enumerate elements stored in Mesh::element.
*
* For all elements of equal dimension to the mesh dimension, we
* enumerate sequentially. All other elements in the list should be of
* lower dimension and have ID set by a corresponding edgeLink or
......@@ -276,12 +292,14 @@ namespace Nektar
}
}
/**
* Each element is assigned to a composite ID by Gmsh. First we scan
* the element list and generate a list of composite IDs. We then
* generate the composite objects and populate them with a second scan
* through the element list.
* @brief Generate a list of composites (groups of elements) from tag
* IDs stored in mesh vertices/edges/faces/elements.
*
* Each element is assigned to a composite ID by an input module. First
* we scan the element list and generate a list of composite IDs. We
* then generate the composite objects and populate them with a second
* scan through the element list.
*/
void Module::ProcessComposites()
{
......@@ -386,6 +404,9 @@ namespace Nektar
}
}
/**
* @brief Print a brief summary of information.
*/
void InputModule::PrintSummary()
{
// Compute the number of full-dimensional elements and boundary
......
......@@ -67,14 +67,28 @@ namespace Nektar
"Output"
};
/**
* @brief Represents a command-line configuration option.
*/
struct ConfigOption
{
/**
* @brief Construct a new configuration option.
*
* @param isBool True if the option is boolean type.
* @param defValue Default value of the option.
* @param desc Description of the option.
*/
ConfigOption(bool isBool, string defValue, string desc) :
isBool(isBool), beenSet(false), defValue(defValue),
desc(desc), value() {}
ConfigOption() :
isBool(false), beenSet(false), defValue(), desc(), value() {}
/**
* @brief Re-interpret the value stored in #value as some type using
* boost::lexical_cast.
*/
template<typename T>
T as()
{
......@@ -89,10 +103,17 @@ namespace Nektar
}
}
/// True if the configuration option is a boolean (thus does not
/// need additional arguments).
bool isBool;
/// True if the configuration option has been set at command
/// line. If false, the default value will be put into #value.
bool beenSet;
/// The value of the configuration option.
string value;
/// Default value of the configuration option.
string defValue;