#include <ProcessSpherigon.h>

Inheritance diagram for Nektar::Utilities::ProcessSpherigon:

Collaboration diagram for Nektar::Utilities::ProcessSpherigon:

Public Member Functions
	ProcessSpherigon (NekMeshUtils::MeshSharedPtr m)
	Default constructor. More...

virtual	~ProcessSpherigon ()
	Destructor. More...

virtual void	Process ()
	Write mesh to output file. More...

Public Member Functions inherited from Nektar::NekMeshUtils::ProcessModule
NEKMESHUTILS_EXPORT	ProcessModule (MeshSharedPtr p_m)

Public Member Functions inherited from Nektar::NekMeshUtils::Module
NEKMESHUTILS_EXPORT	Module (MeshSharedPtr p_m)

NEKMESHUTILS_EXPORT void	RegisterConfig (std::string key, std::string value)
	Register a configuration option with a module. More...

NEKMESHUTILS_EXPORT void	PrintConfig ()
	Print out all configuration options for a module. More...

NEKMESHUTILS_EXPORT void	SetDefaults ()
	Sets default configuration options for those which have not been set. More...

NEKMESHUTILS_EXPORT MeshSharedPtr	GetMesh ()

virtual NEKMESHUTILS_EXPORT void	ProcessVertices ()
	Extract element vertices. More...

virtual NEKMESHUTILS_EXPORT void	ProcessEdges (bool ReprocessEdges=true)
	Extract element edges. More...

virtual NEKMESHUTILS_EXPORT void	ProcessFaces (bool ReprocessFaces=true)
	Extract element faces. More...

virtual NEKMESHUTILS_EXPORT void	ProcessElements ()
	Generate element IDs. More...

virtual NEKMESHUTILS_EXPORT void	ProcessComposites ()
	Generate composites. More...

virtual NEKMESHUTILS_EXPORT void	ClearElementLinks ()

Static Public Member Functions
static boost::shared_ptr< Module >	create (NekMeshUtils::MeshSharedPtr m)
	Creates an instance of this class. More...

Static Public Attributes
static NekMeshUtils::ModuleKey	className

Protected Member Functions
void	GenerateNormals (std::vector< NekMeshUtils::ElementSharedPtr > &el, NekMeshUtils::MeshSharedPtr &mesh)
	Generate a set of approximate vertex normals to a surface represented by line segments in 2D and a hybrid triangular/quadrilateral mesh in 3D. More...

NekDouble	CrossProdMag (NekMeshUtils::Node &a, NekMeshUtils::Node &b)
	Calculate the magnitude of the cross product $\vec{a}\times\vec{b}$ . More...

void	UnitCrossProd (NekMeshUtils::Node &a, NekMeshUtils::Node &b, NekMeshUtils::Node &c)

NekDouble	Blend (NekDouble r)
	Calculate the blending function for spherigon implementation. More...

void	SuperBlend (std::vector< NekDouble > &r, std::vector< NekMeshUtils::Node > &Q, NekMeshUtils::Node &P, std::vector< NekDouble > &blend)
	Calculate the blending function for spherigon implementation. More...

void	FindNormalFromPlyFile (NekMeshUtils::MeshSharedPtr &plymesh, std::map< int, NekMeshUtils::NodeSharedPtr > &surfverts)

Protected Member Functions inherited from Nektar::NekMeshUtils::Module
NEKMESHUTILS_EXPORT void	ReorderPrisms (PerMap &perFaces)
	Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...

NEKMESHUTILS_EXPORT void	PrismLines (int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)

Additional Inherited Members
Protected Attributes inherited from Nektar::NekMeshUtils::Module
MeshSharedPtr	m_mesh
	Mesh object. More...

std::map< std::string, ConfigOption >	m_config
	List of configuration values. More...

Detailed Description

This class implements the spherigon surface smoothing technique which is documented in

"The SPHERIGON: A Simple Polygon Patch for Smoothing Quickly your Polygonal Meshes": P. Volino and N. Magnenat Thalmann, Computer Animation Proceedings (1998).

This implementation works in both a 2D manifold setting (for triangles and quadrilaterals embedded in 3-space) and in a full 3D enviroment (prisms, tetrahedra and hexahedra).

No additional information needs to be supplied directly to the module in order for it to work. However, 3D elements do rely on the mapping Mesh::spherigonSurfs to be populated by the relevant input modules.

Additionally, since the algorithm assumes normals are supplied which are perpendicular to the true surface defined at each vertex. If these are specified in Mesh::vertexNormals by the input module, better smoothing results can be obtained. Otherwise, normals are estimated by taking the average of all edge/face normals which connect to the vertex.

Definition at line 48 of file ProcessSpherigon.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessSpherigon::ProcessSpherigon ( NekMeshUtils::MeshSharedPtr m )

Default constructor.

Definition at line 99 of file ProcessSpherigon.cpp.

References Nektar::NekMeshUtils::Module::m_config.

                                                   : ProcessModule(m)
 {
     m_config["N"] =
         ConfigOption(false, "5", "Number of points to add to face edges.");
     m_config["surf"] =
         ConfigOption(false, "-1", "Tag identifying surface to process.");
     m_config["BothTriFacesOnPrism"] = ConfigOption(
         true, "-1", "Curve both triangular faces of prism on boundary.");
     m_config["usenormalfile"] = ConfigOption(
         false, "NoFile", "Use alternative file for Spherigon definition");
     m_config["scalefile"] = ConfigOption(
         false, "1.0", "Apply scaling factor to coordinates in file ");
     m_config["normalnoise"] =
         ConfigOption(false,
                      "NotSpecified",
                      "Add randowm noise to normals of amplitude AMP "
                      "in specified region. input string is "
                      "Amp,xmin,xmax,ymin,ymax,zmin,zmax");
 }

Nektar::Utilities::ProcessSpherigon::~ProcessSpherigon ( )

virtual

Destructor.

Definition at line 122 of file ProcessSpherigon.cpp.

123 {

124 }

Member Function Documentation

double Nektar::Utilities::ProcessSpherigon::Blend ( NekDouble r )

protected

Calculate the $ C^0 $ blending function for spherigon implementation.

See equation (5) of the paper.

Parameters

r	Barycentric coordinate.

Definition at line 165 of file ProcessSpherigon.cpp.

 {
     return r * r;
 }

static boost::shared_ptr<Module> Nektar::Utilities::ProcessSpherigon::create ( NekMeshUtils::MeshSharedPtr m )

inlinestatic

Creates an instance of this class.

Definition at line 52 of file ProcessSpherigon.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

     {
         return MemoryManager<ProcessSpherigon>::AllocateSharedPtr(m);
     }

double Nektar::Utilities::ProcessSpherigon::CrossProdMag	(	NekMeshUtils::Node &	a,
		NekMeshUtils::Node &	b
	)

protected

Calculate the magnitude of the cross product $\vec{a}\times\vec{b}$ .

Definition at line 149 of file ProcessSpherigon.cpp.

References Nektar::NekMeshUtils::Node::m_x, Nektar::NekMeshUtils::Node::m_y, and Nektar::NekMeshUtils::Node::m_z.

Referenced by Process().

 {
     double tmp1 = a.m_y * b.m_z - a.m_z * b.m_y;
     double tmp2 = a.m_z * b.m_x - a.m_x * b.m_z;
     double tmp3 = a.m_x * b.m_y - a.m_y * b.m_x;
     return sqrt(tmp1 * tmp1 + tmp2 * tmp2 + tmp3 * tmp3);
 }

void Nektar::Utilities::ProcessSpherigon::FindNormalFromPlyFile	(	NekMeshUtils::MeshSharedPtr &	plymesh,
		std::map< int, NekMeshUtils::NodeSharedPtr > &	surfverts
	)

protected

Definition at line 224 of file ProcessSpherigon.cpp.

References ASSERTL1, Nektar::iterator, Nektar::NekMeshUtils::Module::m_mesh, and Nektar::LibUtilities::PrintProgressbar().

Referenced by Process().

 {
     int      cnt = 0;
     int      j;
     int      prog=0,cntmin;
     NodeSet::iterator it;
     map<int, NodeSharedPtr>::iterator vIt;
 
     typedef bg::model::point<NekDouble, 3, bg::cs::cartesian> Point;
     typedef pair<Point, unsigned int> PointI;
 
     int n_neighbs = 5;
 
     map<int,int>  TreeidtoPlyid;
 
     //Fill vertex array into tree format
     vector<PointI> dataPts;
     for (j = 0, it = plymesh->m_vertexSet.begin();
          it != plymesh->m_vertexSet.end();
          ++it, ++j)
     {
         dataPts.push_back(make_pair(Point( (*it)->m_x,
                                            (*it)->m_y,
                                            (*it)->m_z), j));
         TreeidtoPlyid[j] = (*it)->m_id;
     }
 
     //Build tree
     bgi::rtree<PointI, bgi::rstar<16> > rtree;
     rtree.insert(dataPts.begin(), dataPts.end());
 
     //Find neipghbours
     for (cnt = 0, vIt = surfverts.begin(); vIt != surfverts.end();
          ++vIt, ++cnt)
     {
         if(m_mesh->m_verbose)
         {
             prog = LibUtilities::PrintProgressbar(cnt,surfverts.size(),
                                                   "Nearest ply verts",prog);
         }
 
 
         //I dont know why 5 nearest points are searched for when
         //only the nearest point is used for the data
         //was left like this in the ann->boost rewrite (MT 6/11/16)
         Point queryPt(vIt->second->m_x, vIt->second->m_y, vIt->second->m_z);
         n_neighbs  = 5;
         vector<PointI> result;
         rtree.query(bgi::nearest(queryPt, n_neighbs), std::back_inserter(result));
 
         ASSERTL1(bg::distance(result[0].first,queryPt) < bg::distance(result[1].first,queryPt),
             "Assumption that dist values are ordered from smallest to largest is not correct");
 
         cntmin = TreeidtoPlyid[result[0].second];
 
         ASSERTL1(cntmin < plymesh->m_vertexNormals.size(),
                  "cntmin is out of range");
 
         m_mesh->m_vertexNormals[vIt->first] =
             plymesh->m_vertexNormals[cntmin];
     }
 }

void Nektar::Utilities::ProcessSpherigon::GenerateNormals	(	std::vector< NekMeshUtils::ElementSharedPtr > &	el,
		NekMeshUtils::MeshSharedPtr &	mesh
	)

protected

Generate a set of approximate vertex normals to a surface represented by line segments in 2D and a hybrid triangular/quadrilateral mesh in 3D.

This routine approximates the true vertex normals to a surface by averaging the normals of all edges/faces which connect to the vertex. It is better to use the exact surface normals which can be set in Mesh::vertexNormals, but where they are not supplied this routine calculates an approximation for the spherigon implementation.

Parameters

el	Vector of elements denoting the surface mesh.

Definition at line 301 of file ProcessSpherigon.cpp.

References Nektar::NekMeshUtils::Node::abs2(), ASSERTL0, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eSegment, Nektar::LibUtilities::eTriangle, Nektar::iterator, Nektar::NekMeshUtils::Node::m_x, Nektar::NekMeshUtils::Node::m_y, Nektar::NekMeshUtils::Node::m_z, and UnitCrossProd().

Referenced by Process().

 {
     boost::unordered_map<int, Node>::iterator nIt;
 
     for (int i = 0; i < el.size(); ++i)
     {
         ElementSharedPtr e = el[i];
 
         // Ensure that element is a line, triangle or quad.
         ASSERTL0(e->GetConf().m_e == LibUtilities::eSegment ||
                      e->GetConf().m_e == LibUtilities::eTriangle ||
                      e->GetConf().m_e == LibUtilities::eQuadrilateral,
                  "Spherigon expansions must be lines, triangles or "
                  "quadrilaterals.");
 
         // Calculate normal for this element.
         int nV = e->GetVertexCount();
         vector<NodeSharedPtr> node(nV);
 
         for (int j = 0; j < nV; ++j)
         {
             node[j] = e->GetVertex(j);
         }
 
         Node n;
 
         if (mesh->m_spaceDim == 3)
         {
             // Create two tangent vectors and take unit cross product.
             Node v1 = *(node[1]) - *(node[0]);
             Node v2 = *(node[2]) - *(node[0]);
             UnitCrossProd(v1, v2, n);
         }
         else
         {
             // Calculate gradient vector and invert.
             Node dx = *(node[1]) - *(node[0]);
             NekDouble diff = dx.abs2();
             dx /= sqrt(diff);
             n.m_x = -dx.m_y;
             n.m_y = dx.m_x;
             n.m_z = 0;
         }
 
         // Insert face normal into vertex normal list or add to existing
         // value.
         for (int j = 0; j < nV; ++j)
         {
             nIt = mesh->m_vertexNormals.find(e->GetVertex(j)->m_id);
             if (nIt == mesh->m_vertexNormals.end())
             {
                 mesh->m_vertexNormals[e->GetVertex(j)->m_id] = n;
             }
             else
             {
                 nIt->second += n;
             }
         }
     }
 
     // Normalize resulting vectors.
     for (nIt = mesh->m_vertexNormals.begin();
          nIt != mesh->m_vertexNormals.end();
          ++nIt)
     {
         Node &n = mesh->m_vertexNormals[nIt->first];
         n /= sqrt(n.abs2());
     }
 
 }

void Nektar::Utilities::ProcessSpherigon::Process ( )

virtual

Write mesh to output file.

Perform the spherigon smoothing technique on the mesh.

Implements Nektar::NekMeshUtils::Module.

Definition at line 376 of file ProcessSpherigon.cpp.

 {
     ASSERTL0(m_mesh->m_spaceDim == 3 || m_mesh->m_spaceDim == 2,
              "Spherigon implementation only valid in 2D/3D.");
 
     boost::unordered_set<int>::iterator eIt;
     boost::unordered_set<int> visitedEdges;
 
     // First construct vector of elements to process.
     vector<ElementSharedPtr> el;
 
     if (m_mesh->m_verbose)
     {
         cout << "ProcessSpherigon: Smoothing mesh..." << endl;
     }
 
     if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3)
     {
         // Manifold case - copy expansion dimension.
         el = m_mesh->m_element[m_mesh->m_expDim];
     }
     else if (m_mesh->m_spaceDim == m_mesh->m_expDim)
     {
         // Full 2D or 3D case - iterate over stored edges/faces and
         // create segments/triangles representing those edges/faces.
         set<pair<int, int> >::iterator it;
         vector<int> t;
         t.push_back(0);
 
         // Construct list of spherigon edges/faces from a tag.
         string surfTag = m_config["surf"].as<string>();
         bool prismTag  = m_config["BothTriFacesOnPrism"].beenSet;
 
         if (surfTag != "")
         {
             vector<unsigned int> surfs;
             ParseUtils::GenerateSeqVector(surfTag.c_str(), surfs);
             sort(surfs.begin(), surfs.end());
 
             m_mesh->m_spherigonSurfs.clear();
             for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
             {
                 ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
                 int nSurf = m_mesh->m_expDim == 3 ? el->GetFaceCount()
                                                   : el->GetEdgeCount();
 
                 for (int j = 0; j < nSurf; ++j)
                 {
                     int bl = el->GetBoundaryLink(j);
                     if (bl == -1)
                     {
                         continue;
                     }
 
                     ElementSharedPtr bEl =
                         m_mesh->m_element[m_mesh->m_expDim - 1][bl];
                     vector<int> tags = bEl->GetTagList();
                     vector<int> inter;
 
                     sort(tags.begin(), tags.end());
                     set_intersection(surfs.begin(),
                                      surfs.end(),
                                      tags.begin(),
                                      tags.end(),
                                      back_inserter(inter));
 
                     if (inter.size() == 1)
                     {
                         m_mesh->m_spherigonSurfs.insert(make_pair(i, j));
 
                         // Curve other tri face on Prism. Note could be
                         // problem on pyramid when implemented.
                         if (nSurf == 5 && prismTag)
                         {
                             // add other end of prism on boundary for
                             // smoothing
                             int triFace = j == 1 ? 3 : 1;
                             m_mesh->m_spherigonSurfs.insert(
                                 make_pair(i, triFace));
                         }
                     }
                 }
             }
         }
 
         if (m_mesh->m_spherigonSurfs.size() == 0)
         {
             cerr << "WARNING: Spherigon surfaces have not been defined "
                  << "-- ignoring smoothing." << endl;
             return;
         }
 
         if (m_mesh->m_expDim == 3)
         {
             for (it = m_mesh->m_spherigonSurfs.begin();
                  it != m_mesh->m_spherigonSurfs.end();
                  ++it)
             {
                 FaceSharedPtr f =
                     m_mesh->m_element[m_mesh->m_expDim][it->first]->GetFace(
                         it->second);
                 vector<NodeSharedPtr> nodes = f->m_vertexList;
                 LibUtilities::ShapeType eType =
                     (LibUtilities::ShapeType)(nodes.size());
                 ElmtConfig conf(eType, 1, false, false);
 
                 // Create 2D element.
                 ElementSharedPtr elmt =
                     GetElementFactory().CreateInstance(eType, conf, nodes, t);
 
                 // Copy vertices/edges from face.
                 for (int i = 0; i < f->m_vertexList.size(); ++i)
                 {
                     elmt->SetVertex(i, f->m_vertexList[i]);
                 }
                 for (int i = 0; i < f->m_edgeList.size(); ++i)
                 {
                     elmt->SetEdge(i, f->m_edgeList[i]);
                 }
 
                 el.push_back(elmt);
             }
         }
         else
         {
             for (it = m_mesh->m_spherigonSurfs.begin();
                  it != m_mesh->m_spherigonSurfs.end();
                  ++it)
             {
                 EdgeSharedPtr edge =
                     m_mesh->m_element[m_mesh->m_expDim][it->first]->GetEdge(
                         it->second);
                 vector<NodeSharedPtr> nodes;
                 LibUtilities::ShapeType eType = LibUtilities::eSegment;
                 ElmtConfig conf(eType, 1, false, false);
 
                 nodes.push_back(edge->m_n1);
                 nodes.push_back(edge->m_n2);
 
                 // Create 2D element.
                 ElementSharedPtr elmt =
                     GetElementFactory().CreateInstance(eType, conf, nodes, t);
 
                 // Copy vertices/edges from original element.
                 elmt->SetVertex(0, nodes[0]);
                 elmt->SetVertex(1, nodes[1]);
                 elmt->SetEdge(
                     0,
                     m_mesh->m_element[m_mesh->m_expDim][it->first]->GetEdge(
                         it->second));
                 el.push_back(elmt);
             }
         }
     }
     else
     {
         ASSERTL0(false, "Spherigon expansions must be 2/3 dimensional");
     }
 
     // See if vertex normals have been generated. If they have not,
     // approximate them by summing normals of surrounding elements.
     bool normalsGenerated = false;
 
     // Read Normal file if one exists
     std::string normalfile = m_config["usenormalfile"].as<string>();
     if (normalfile.compare("NoFile") != 0)
     {
         NekDouble scale = m_config["scalefile"].as<NekDouble>();
 
         if (m_mesh->m_verbose)
         {
             cout << "Inputing normal file: " << normalfile
                  << " with scaling of " << scale << endl;
         }
 
         ifstream inplyTmp;
         io::filtering_istream inply;
         InputPlySharedPtr plyfile;
 
         inplyTmp.open(normalfile.c_str());
         ASSERTL0(inplyTmp,
                  string("Could not open input ply file: ") + normalfile);
 
         inply.push(inplyTmp);
 
         MeshSharedPtr m = boost::shared_ptr<Mesh>(new Mesh());
         plyfile = boost::shared_ptr<InputPly>(new InputPly(m));
         plyfile->ReadPly(inply, scale);
         plyfile->ProcessVertices();
 
         MeshSharedPtr plymesh = plyfile->GetMesh();
         if (m_mesh->m_verbose)
         {
             cout << "\t Generating ply normals" << endl;
         }
         GenerateNormals(plymesh->m_element[plymesh->m_expDim], plymesh);
 
         // finally find nearest vertex and set normal to mesh surface file
         // normal.  probably should have a hex tree search ?
         Node minx(0, 0.0, 0.0, 0.0), tmp, tmpsav;
         NodeSet::iterator it;
         map<int, NodeSharedPtr>::iterator vIt;
         map<int, NodeSharedPtr> surfverts;
 
         // make a map of normal vertices to visit based on elements el
         for (int i = 0; i < el.size(); ++i)
         {
             ElementSharedPtr e = el[i];
             int nV = e->GetVertexCount();
             for (int j = 0; j < nV; ++j)
             {
                 int id        = e->GetVertex(j)->m_id;
                 surfverts[id] = e->GetVertex(j);
             }
         }
 
 
         if (m_mesh->m_verbose)
         {
             cout << "\t Processing surface normals " << endl;
         }
 
         // loop over all element in ply mesh and determine
         // and set normal to nearest point in ply mesh
         FindNormalFromPlyFile(plymesh,surfverts);
 
         normalsGenerated = true;
     }
     else if (m_mesh->m_vertexNormals.size() == 0)
     {
         GenerateNormals(el, m_mesh);
         normalsGenerated = true;
     }
 
     // See if we should add noise to normals
     std::string normalnoise = m_config["normalnoise"].as<string>();
     if (normalnoise.compare("NotSpecified") != 0)
     {
         vector<NekDouble> values;
         ASSERTL0(
             ParseUtils::GenerateUnOrderedVector(normalnoise.c_str(), values),
             "Failed to interpret normal noise string");
 
         int nvalues   = values.size() / 2;
         NekDouble amp = values[0];
 
         if (m_mesh->m_verbose)
         {
             cout << "\t adding noise to normals of amplitude " << amp
                  << " in range: ";
             for (int i = 0; i < nvalues; ++i)
             {
                 cout << values[2 * i + 1] << "," << values[2 * i + 2] << " ";
             }
             cout << endl;
         }
 
         map<int, NodeSharedPtr>::iterator vIt;
         map<int, NodeSharedPtr> surfverts;
 
         // make a map of normal vertices to visit based on elements el
         for (int i = 0; i < el.size(); ++i)
         {
             ElementSharedPtr e = el[i];
             int nV = e->GetVertexCount();
             for (int j = 0; j < nV; ++j)
             {
                 int id        = e->GetVertex(j)->m_id;
                 surfverts[id] = e->GetVertex(j);
             }
         }
 
         for (vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt)
         {
             bool AddNoise = false;
 
             for (int i = 0; i < nvalues; ++i)
             {
                 // check to see if point is in range
                 switch (nvalues)
                 {
                     case 1:
                     {
                         if (((vIt->second)->m_x > values[2 * i + 1]) &&
                             ((vIt->second)->m_x < values[2 * i + 2]))
                         {
                             AddNoise = true;
                         }
                     }
                     break;
                     case 2:
                     {
                         if (((vIt->second)->m_x > values[2 * i + 1]) &&
                             ((vIt->second)->m_x < values[2 * i + 2]) &&
                             ((vIt->second)->m_y > values[2 * i + 3]) &&
                             ((vIt->second)->m_y < values[2 * i + 4]))
                         {
                             AddNoise = true;
                         }
                     }
                     break;
                     case 3:
                     {
                         if (((vIt->second)->m_x > values[2 * i + 1]) &&
                             ((vIt->second)->m_x < values[2 * i + 2]) &&
                             ((vIt->second)->m_y > values[2 * i + 3]) &&
                             ((vIt->second)->m_y < values[2 * i + 4]) &&
                             ((vIt->second)->m_z > values[2 * i + 5]) &&
                             ((vIt->second)->m_z < values[2 * i + 6]))
 
                         {
                             AddNoise = true;
                         }
                     }
                     break;
                 }
 
                 if (AddNoise)
                 {
                     // generate random unit vector;
                     Node rvec(0, rand(), rand(), rand());
                     rvec *= values[0] / sqrt(rvec.abs2());
 
                     Node normal = m_mesh->m_vertexNormals[vIt->first];
 
                     normal += rvec;
                     normal /= sqrt(normal.abs2());
 
                     m_mesh->m_vertexNormals[vIt->first] = normal;
                 }
             }
         }
     }
 
     // Allocate storage for interior points.
     int nq    = m_config["N"].as<int>();
     int nquad = m_mesh->m_spaceDim == 3 ? nq * nq : nq;
     Array<OneD, NekDouble> x(nq * nq);
     Array<OneD, NekDouble> y(nq * nq);
     Array<OneD, NekDouble> z(nq * nq);
 
     Array<OneD, NekDouble> xc(nq * nq);
     Array<OneD, NekDouble> yc(nq * nq);
     Array<OneD, NekDouble> zc(nq * nq);
 
     ASSERTL0(nq > 2, "Number of points must be greater than 2.");
 
     LibUtilities::BasisKey B0(
         LibUtilities::eOrtho_A,
         nq,
         LibUtilities::PointsKey(nq, LibUtilities::eGaussLobattoLegendre));
     LibUtilities::BasisKey B1(
         LibUtilities::eOrtho_B,
         nq,
         LibUtilities::PointsKey(nq, LibUtilities::eGaussRadauMAlpha1Beta0));
     StdRegions::StdNodalTriExpSharedPtr stdtri =
         MemoryManager<StdRegions::StdNodalTriExp>::AllocateSharedPtr(
             B0, B1, LibUtilities::eNodalTriElec);
 
     Array<OneD, NekDouble> xnodal(nq * (nq + 1) / 2), ynodal(nq * (nq + 1) / 2);
     stdtri->GetNodalPoints(xnodal, ynodal);
 
     int edgeMap[3][4][2] = {
         {{0, 1}, {-1, -1}, {-1, -1}, {-1, -1}},                         // seg
         {{0, 1}, {nq - 1, nq}, {nq * (nq - 1), -nq}, {-1, -1}},         // tri
         {{0, 1}, {nq - 1, nq}, {nq * nq - 1, -1}, {nq * (nq - 1), -nq}} // quad
     };
 
     int vertMap[3][4][2] = {
         {{0, 1}, {0, 0}, {0, 0}, {0, 0}}, // seg
         {{0, 1}, {1, 2}, {2, 3}, {0, 0}}, // tri
         {{0, 1}, {1, 2}, {2, 3}, {3, 0}}, // quad
     };
 
     for (int i = 0; i < el.size(); ++i)
     {
         // Construct a Nektar++ element to obtain coordinate points
         // inside the element. TODO: Add options for various
         // nodal/tensor point distributions + number of points to add.
         ElementSharedPtr e = el[i];
 
         LibUtilities::BasisKey B2(
             LibUtilities::eModified_A,
             nq,
             LibUtilities::PointsKey(nq, LibUtilities::eGaussLobattoLegendre));
 
         if (e->GetConf().m_e == LibUtilities::eSegment)
         {
             SpatialDomains::SegGeomSharedPtr geom =
                 boost::dynamic_pointer_cast<SpatialDomains::SegGeom>(
                     e->GetGeom(m_mesh->m_spaceDim));
             LocalRegions::SegExpSharedPtr seg =
                 MemoryManager<LocalRegions::SegExp>::AllocateSharedPtr(B2,
                                                                        geom);
             seg->GetCoords(x, y, z);
             nquad = nq;
         }
         else if (e->GetConf().m_e == LibUtilities::eTriangle)
         {
             SpatialDomains::TriGeomSharedPtr geom =
                 boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
                     e->GetGeom(3));
             LocalRegions::NodalTriExpSharedPtr tri =
                 MemoryManager<LocalRegions::NodalTriExp>::AllocateSharedPtr(
                     B0, B1, LibUtilities::eNodalTriElec, geom);
 
             Array<OneD, NekDouble> coord(2);
             tri->GetCoords(xc, yc, zc);
             nquad = nq * (nq + 1) / 2;
 
             for (int j = 0; j < nquad; ++j)
             {
                 coord[0] = xnodal[j];
                 coord[1] = ynodal[j];
                 x[j]     = stdtri->PhysEvaluate(coord, xc);
             }
 
             for (int j = 0; j < nquad; ++j)
             {
                 coord[0] = xnodal[j];
                 coord[1] = ynodal[j];
                 y[j]     = stdtri->PhysEvaluate(coord, yc);
             }
 
             for (int j = 0; j < nquad; ++j)
             {
                 coord[0] = xnodal[j];
                 coord[1] = ynodal[j];
                 z[j]     = stdtri->PhysEvaluate(coord, zc);
             }
         }
         else if (e->GetConf().m_e == LibUtilities::eQuadrilateral)
         {
             SpatialDomains::QuadGeomSharedPtr geom =
                 boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
                     e->GetGeom(3));
             LocalRegions::QuadExpSharedPtr quad =
                 MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(
                     B2, B2, geom);
             quad->GetCoords(x, y, z);
             nquad = nq * nq;
         }
         else
         {
             ASSERTL0(false, "Unknown expansion type.");
         }
 
         // Zero z-coordinate in 2D.
         if (m_mesh->m_spaceDim == 2)
         {
             Vmath::Zero(nquad, z, 1);
         }
 
         // Find vertex normals.
         int nV = e->GetVertexCount();
         vector<Node> v, vN;
         for (int j = 0; j < nV; ++j)
         {
             v.push_back(*(e->GetVertex(j)));
             ASSERTL1(m_mesh->m_vertexNormals.count(v[j].m_id) != 0,
                      "Normal has not been defined");
             vN.push_back(m_mesh->m_vertexNormals[v[j].m_id]);
         }
 
         vector<Node> tmp(nV);
         vector<NekDouble> r(nV);
         vector<Node> K(nV);
         vector<Node> Q(nV);
         vector<Node> Qp(nV);
         vector<NekDouble> blend(nV);
         vector<Node> out(nquad);
 
         // Calculate segment length for 2D spherigon routine.
         NekDouble segLength = sqrt((v[0] - v[1]).abs2());
 
         // Perform Spherigon method to smooth manifold.
         for (int j = 0; j < nquad; ++j)
         {
             Node P(0, x[j], y[j], z[j]);
             Node N(0, 0, 0, 0);
 
             // Calculate generalised barycentric coordinates r[] and the
             // Phong normal N = vN . r for this point of the element.
             if (m_mesh->m_spaceDim == 2)
             {
                 // In 2D the coordinates are given by a ratio of the
                 // segment length to the distance from one of the
                 // endpoints.
                 r[0] = sqrt((P - v[0]).abs2()) / segLength;
                 r[0] = max(min(1.0, r[0]), 0.0);
                 r[1] = 1.0 - r[0];
 
                 // Calculate Phong normal.
                 N = vN[0] * r[0] + vN[1] * r[1];
             }
             else if (m_mesh->m_spaceDim == 3)
             {
                 for (int k = 0; k < nV; ++k)
                 {
                     tmp[k] = P - v[k];
                 }
 
                 // Calculate generalized barycentric coordinate system
                 // (see equation 6 of paper).
                 NekDouble weight = 0.0;
                 for (int k = 0; k < nV; ++k)
                 {
                     r[k] = 1.0;
                     for (int l = 0; l < nV - 2; ++l)
                     {
                         r[k] *= CrossProdMag(tmp[(k + l + 1) % nV],
                                              tmp[(k + l + 2) % nV]);
                     }
                     weight += r[k];
                 }
 
                 // Calculate Phong normal (equation 1).
                 for (int k = 0; k < nV; ++k)
                 {
                     r[k] /= weight;
                     N += vN[k] * r[k];
                 }
             }
 
             // Normalise Phong normal.
             N /= sqrt(N.abs2());
 
             for (int k = 0; k < nV; ++k)
             {
                 // Perform steps denoted in equations 2, 3, 8 for C1
                 // smoothing.
                 NekDouble tmp1;
                 K[k]  = P + N * ((v[k] - P).dot(N));
                 tmp1  = (v[k] - K[k]).dot(vN[k]) / (1.0 + N.dot(vN[k]));
                 Q[k]  = K[k] + N * tmp1;
                 Qp[k] = v[k] - N * ((v[k] - P).dot(N));
             }
 
             // Apply C1 blending function to the surface. TODO: Add
             // option to do (more efficient) C0 blending function.
             SuperBlend(r, Qp, P, blend);
             P.m_x = P.m_y = P.m_z = 0.0;
 
             // Apply blending (equation 4).
             for (int k = 0; k < nV; ++k)
             {
                 P += Q[k] * blend[k];
             }
 
             if ((boost::math::isnan)(P.m_x) || (boost::math::isnan)(P.m_y) ||
                 (boost::math::isnan)(P.m_z))
             {
                 ASSERTL0(false,
                          "spherigon point is a nan. Check to see if "
                          "ply file is correct if using input normal file");
             }
             else
             {
                 out[j] = P;
             }
         }
 
         // Push nodes into lines - TODO: face interior nodes.
         // offset = 0 (seg), 1 (tri) or 2 (quad)
         int offset = (int)e->GetConf().m_e - 2;
 
         for (int edge = 0; edge < e->GetEdgeCount(); ++edge)
         {
             eIt = visitedEdges.find(e->GetEdge(edge)->m_id);
             if (eIt == visitedEdges.end())
             {
                 bool reverseEdge =
                     !(v[vertMap[offset][edge][0]] == *(e->GetEdge(edge)->m_n1));
 
                 // Clear existing curvature.
                 e->GetEdge(edge)->m_edgeNodes.clear();
 
                 if (e->GetConf().m_e != LibUtilities::eTriangle)
                 {
                     for (int j = 1; j < nq - 1; ++j)
                     {
                         int v = edgeMap[offset][edge][0] +
                                 j * edgeMap[offset][edge][1];
                         e->GetEdge(edge)->m_edgeNodes.push_back(
                             NodeSharedPtr(new Node(out[v])));
                     }
                 }
                 else
                 {
                     for (int j = 0; j < nq - 2; ++j)
                     {
                         int v = 3 + edge * (nq - 2) + j;
                         e->GetEdge(edge)->m_edgeNodes.push_back(
                             NodeSharedPtr(new Node(out[v])));
                     }
                 }
 
                 if (reverseEdge)
                 {
                     reverse(e->GetEdge(edge)->m_edgeNodes.begin(),
                             e->GetEdge(edge)->m_edgeNodes.end());
                 }
 
                 e->GetEdge(edge)->m_curveType =
                     LibUtilities::eGaussLobattoLegendre;
 
                 visitedEdges.insert(e->GetEdge(edge)->m_id);
             }
         }
 
         // Add face nodes in manifold and full 3D case, but not for 2D.
         if (m_mesh->m_spaceDim == 3)
         {
             vector<NodeSharedPtr> volNodes;
 
             if (e->GetConf().m_e == LibUtilities::eQuadrilateral)
             {
                 volNodes.resize((nq - 2) * (nq - 2));
                 for (int j = 1; j < nq - 1; ++j)
                 {
                     for (int k = 1; k < nq - 1; ++k)
                     {
                         int v = j * nq + k;
                         volNodes[(j - 1) * (nq - 2) + (k - 1)] =
                             NodeSharedPtr(new Node(out[v]));
                     }
                 }
             }
             else
             {
                 for (int j = 3 + 3 * (nq - 2); j < nquad; ++j)
                 {
                     volNodes.push_back(NodeSharedPtr(new Node(out[j])));
                 }
             }
 
             e->SetVolumeNodes(volNodes);
         }
     }
 
     // Copy face nodes back into 3D element faces
     if (m_mesh->m_expDim == 3)
     {
         set<pair<int, int> >::iterator it;
         int elmt = 0;
         for (it = m_mesh->m_spherigonSurfs.begin();
              it != m_mesh->m_spherigonSurfs.end();
              ++it, ++elmt)
         {
             FaceSharedPtr f =
                 m_mesh->m_element[m_mesh->m_expDim][it->first]->GetFace(
                     it->second);
 
             f->m_faceNodes = el[elmt]->GetVolumeNodes();
             f->m_curveType = f->m_vertexList.size() == 3
                                  ? LibUtilities::eNodalTriElec
                                  : LibUtilities::eGaussLobattoLegendre;
         }
     }
 
     if (normalsGenerated)
     {
         m_mesh->m_vertexNormals.clear();
     }
 }

void Nektar::Utilities::ProcessSpherigon::SuperBlend	(	std::vector< NekDouble > &	r,
		std::vector< NekMeshUtils::Node > &	Q,
		NekMeshUtils::Node &	P,
		std::vector< NekDouble > &	blend
	)

protected

Calculate the $ C^1 $ blending function for spherigon implementation.

See equation (10) of the paper.

Parameters

r	Generalised barycentric coordinates of the point P.
Q	Vector of vertices denoting this triangle/quad.
P	Point in the triangle to apply blending to.
blend	The resulting blending components for each vertex.

Definition at line 181 of file ProcessSpherigon.cpp.

References class_topology::P, and TOL_BLEND.

Referenced by Process().

 {
     vector<double> tmp(r.size());
     double totBlend = 0.0;
     int i;
     int nV = r.size();
 
     for (i = 0; i < nV; ++i)
     {
         blend[i] = 0.0;
         tmp[i]   = (Q[i] - P).abs2();
     }
 
     for (i = 0; i < nV; ++i)
     {
         int ip = (i + 1) % nV, im = (i - 1 + nV) % nV;
 
         if (r[i] > TOL_BLEND && r[i] < (1 - TOL_BLEND))
         {
             blend[i] =
                 r[i] * r[i] * (r[im] * r[im] * tmp[im] / (tmp[im] + tmp[i]) +
                                r[ip] * r[ip] * tmp[ip] / (tmp[ip] + tmp[i]));
             totBlend += blend[i];
         }
     }
 
     for (i = 0; i < nV; ++i)
     {
         blend[i] /= totBlend;
         if (r[i] >= (1 - TOL_BLEND))
         {
             blend[i] = 1.0;
         }
         if (r[i] <= TOL_BLEND)
         {
             blend[i] = 0.0;
         }
     }
 }

void Nektar::Utilities::ProcessSpherigon::UnitCrossProd	(	NekMeshUtils::Node &	a,
		NekMeshUtils::Node &	b,
		NekMeshUtils::Node &	c
	)

protected

Definition at line 130 of file ProcessSpherigon.cpp.

References Nektar::NekMeshUtils::Node::m_x, Nektar::NekMeshUtils::Node::m_y, and Nektar::NekMeshUtils::Node::m_z.

Referenced by GenerateNormals().

 {
     double inv_mag;
 
     c.m_x = a.m_y * b.m_z - a.m_z * b.m_y;
     c.m_y = a.m_z * b.m_x - a.m_x * b.m_z;
     c.m_z = a.m_x * b.m_y - a.m_y * b.m_x;
 
     inv_mag = 1.0 / sqrt(c.m_x * c.m_x + c.m_y * c.m_y + c.m_z * c.m_z);
 
     c.m_x = c.m_x * inv_mag;
     c.m_y = c.m_y * inv_mag;
     c.m_z = c.m_z * inv_mag;
 }

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessSpherigon::className

static

Initial value:

=
    GetModuleFactory().RegisterCreatorFunction(
        ModuleKey(eProcessModule, "spherigon"), ProcessSpherigon::create)

Definition at line 56 of file ProcessSpherigon.h.

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation