#include <ProcessSpherigon.h>

Inheritance diagram for Nektar::Utilities::ProcessSpherigon:

Collaboration diagram for Nektar::Utilities::ProcessSpherigon:

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

virtual	~ProcessSpherigon ()
	Destructor. More...

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

Public Member Functions inherited from Nektar::Utilities::ProcessModule
	ProcessModule ()

	ProcessModule (FieldSharedPtr p_f)

	ProcessModule (MeshSharedPtr p_m)

Public Member Functions inherited from Nektar::Utilities::Module
	Module (FieldSharedPtr p_f)

virtual void	Process (po::variables_map &vm)=0

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

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

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

bool	GetRequireEquiSpaced (void)

void	SetRequireEquiSpaced (bool pVal)

void	EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)

	Module (MeshSharedPtr p_m)

void	RegisterConfig (string key, string value)

void	PrintConfig ()

void	SetDefaults ()

MeshSharedPtr	GetMesh ()

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

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

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

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

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

virtual void	ClearElementLinks ()

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

Static Public Attributes
static ModuleKey	className

Protected Member Functions
void	GenerateNormals (vector< ElementSharedPtr > &el, 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...

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

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

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

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

Protected Member Functions inherited from Nektar::Utilities::Module
	Module ()

void	ReorderPrisms (PerMap &perFaces)
	Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...

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

Additional Inherited Members
Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr	m_f
	Field object. More...

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

bool	m_requireEquiSpaced

MeshSharedPtr	m_mesh
	Mesh object. 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 ( MeshSharedPtr m )

Default constructor.

Definition at line 91 of file ProcessSpherigon.cpp.

References Nektar::Utilities::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 114 of file ProcessSpherigon.cpp.

115 {

116 }

Member Function Documentation

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

protected

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

See equation (5) of the paper.

Parameters

r	Barycentric coordinate.

Definition at line 157 of file ProcessSpherigon.cpp.

 {
     return r * r;
 }

static boost::shared_ptr<Module> Nektar::Utilities::ProcessSpherigon::create ( 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	(	Node &	a,
		Node &	b
	)

protected

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

Definition at line 141 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::GenerateNormals	(	std::vector< ElementSharedPtr > &	el,
		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 229 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]);
             dx /= sqrt(dx.abs2());
             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::Utilities::Module.

Definition at line 302 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 inply;
         InputPlySharedPtr plyfile;
 
         inply.open(normalfile.c_str());
         ASSERTL0(inply, string("Could not open input ply file: ") + normalfile);
 
         int j;
         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();
         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 ?
         Array<OneD, NekDouble> len2(plymesh->m_vertexSet.size());
         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);
             }
         }
 
         // loop over all element in ply mesh and determine
         // xmin,xmax,ymin,ymax as search criterion
 
         NekDouble mindiff, diff;
         int cntmin;
 
         if (m_mesh->m_verbose)
         {
             cout << "\t Processing surface normals " << endl;
         }
         int cnt = 0;
         map<int, int> locnorm;
         for (vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt, ++cnt)
         {
             mindiff = 1e12;
 
             for (j = 0, it = plymesh->m_vertexSet.begin();
                  it != plymesh->m_vertexSet.end();
                  ++it, ++j)
             {
                 tmp  = *(vIt->second) - *(*it);
                 diff = tmp.abs2();
 
                 if (diff < mindiff)
                 {
                     mindiff = diff;
                     cntmin  = (*it)->m_id;
                     tmpsav  = tmp;
                 }
             }
             locnorm[cntmin] = vIt->first;
 
             ASSERTL1(cntmin < plymesh->m_vertexNormals.size(),
                      "cntmin is out of range");
             m_mesh->m_vertexNormals[vIt->first] =
                 plymesh->m_vertexNormals[cntmin];
         }
         if (m_mesh->m_verbose)
         {
             cout << "\t end of processing surface normals " << endl;
         }
         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<double> r(nV);
         vector<Node> K(nV);
         vector<Node> Q(nV);
         vector<Node> Qp(nV);
         vector<double> blend(nV);
         vector<Node> out(nquad);
 
         // Calculate segment length for 2D spherigon routine.
         double 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).
                 double 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.
                 double 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	(	vector< double > &	r,
		vector< Node > &	Q,
		Node &	P,
		vector< double > &	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 173 of file ProcessSpherigon.cpp.

References 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	(	Node &	a,
		Node &	b,
		Node &	c
	)

protected

Definition at line 122 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