Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | List of all members
Nektar::Utilities::ProcessSpherigon Class Reference

#include <ProcessSpherigon.h>

Inheritance diagram for Nektar::Utilities::ProcessSpherigon:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::ProcessSpherigon:
Collaboration graph
[legend]

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< Modulecreate (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 $ C^0 $ 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 $ C^1 $ 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.

99  : ProcessModule(m)
100 {
101  m_config["N"] =
102  ConfigOption(false, "5", "Number of points to add to face edges.");
103  m_config["surf"] =
104  ConfigOption(false, "-1", "Tag identifying surface to process.");
105  m_config["BothTriFacesOnPrism"] = ConfigOption(
106  true, "-1", "Curve both triangular faces of prism on boundary.");
107  m_config["usenormalfile"] = ConfigOption(
108  false, "NoFile", "Use alternative file for Spherigon definition");
109  m_config["scalefile"] = ConfigOption(
110  false, "1.0", "Apply scaling factor to coordinates in file ");
111  m_config["normalnoise"] =
112  ConfigOption(false,
113  "NotSpecified",
114  "Add randowm noise to normals of amplitude AMP "
115  "in specified region. input string is "
116  "Amp,xmin,xmax,ymin,ymax,zmin,zmax");
117 }
Represents a command-line configuration option.
NEKMESHUTILS_EXPORT ProcessModule(MeshSharedPtr p_m)
std::map< std::string, ConfigOption > m_config
List of configuration values.
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
rBarycentric coordinate.

Definition at line 165 of file ProcessSpherigon.cpp.

166 {
167  return r * r;
168 }
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().

53  {
55  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
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().

150 {
151  double tmp1 = a.m_y * b.m_z - a.m_z * b.m_y;
152  double tmp2 = a.m_z * b.m_x - a.m_x * b.m_z;
153  double tmp3 = a.m_x * b.m_y - a.m_y * b.m_x;
154  return sqrt(tmp1 * tmp1 + tmp2 * tmp2 + tmp3 * tmp3);
155 }
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().

226 {
227  int cnt = 0;
228  int j;
229  int prog=0,cntmin;
232 
233  typedef bg::model::point<NekDouble, 3, bg::cs::cartesian> Point;
234  typedef pair<Point, unsigned int> PointI;
235 
236  int n_neighbs = 5;
237 
238  map<int,int> TreeidtoPlyid;
239 
240  //Fill vertex array into tree format
241  vector<PointI> dataPts;
242  for (j = 0, it = plymesh->m_vertexSet.begin();
243  it != plymesh->m_vertexSet.end();
244  ++it, ++j)
245  {
246  dataPts.push_back(make_pair(Point( (*it)->m_x,
247  (*it)->m_y,
248  (*it)->m_z), j));
249  TreeidtoPlyid[j] = (*it)->m_id;
250  }
251 
252  //Build tree
253  bgi::rtree<PointI, bgi::rstar<16> > rtree;
254  rtree.insert(dataPts.begin(), dataPts.end());
255 
256  //Find neipghbours
257  for (cnt = 0, vIt = surfverts.begin(); vIt != surfverts.end();
258  ++vIt, ++cnt)
259  {
260  if(m_mesh->m_verbose)
261  {
262  prog = LibUtilities::PrintProgressbar(cnt,surfverts.size(),
263  "Nearest ply verts",prog);
264  }
265 
266 
267  //I dont know why 5 nearest points are searched for when
268  //only the nearest point is used for the data
269  //was left like this in the ann->boost rewrite (MT 6/11/16)
270  Point queryPt(vIt->second->m_x, vIt->second->m_y, vIt->second->m_z);
271  n_neighbs = 5;
272  vector<PointI> result;
273  rtree.query(bgi::nearest(queryPt, n_neighbs), std::back_inserter(result));
274 
275  ASSERTL1(bg::distance(result[0].first,queryPt) < bg::distance(result[1].first,queryPt),
276  "Assumption that dist values are ordered from smallest to largest is not correct");
277 
278  cntmin = TreeidtoPlyid[result[0].second];
279 
280  ASSERTL1(cntmin < plymesh->m_vertexNormals.size(),
281  "cntmin is out of range");
282 
283  m_mesh->m_vertexNormals[vIt->first] =
284  plymesh->m_vertexNormals[cntmin];
285  }
286 }
int PrintProgressbar(const int position, const int goal, const string message, int lastprogress=-1)
Prints a progressbar.
Definition: Progressbar.hpp:69
A 0-dimensional vertex.
Definition: Point.h:49
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
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
elVector 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().

303 {
305 
306  for (int i = 0; i < el.size(); ++i)
307  {
308  ElementSharedPtr e = el[i];
309 
310  // Ensure that element is a line, triangle or quad.
311  ASSERTL0(e->GetConf().m_e == LibUtilities::eSegment ||
312  e->GetConf().m_e == LibUtilities::eTriangle ||
313  e->GetConf().m_e == LibUtilities::eQuadrilateral,
314  "Spherigon expansions must be lines, triangles or "
315  "quadrilaterals.");
316 
317  // Calculate normal for this element.
318  int nV = e->GetVertexCount();
319  vector<NodeSharedPtr> node(nV);
320 
321  for (int j = 0; j < nV; ++j)
322  {
323  node[j] = e->GetVertex(j);
324  }
325 
326  Node n;
327 
328  if (mesh->m_spaceDim == 3)
329  {
330  // Create two tangent vectors and take unit cross product.
331  Node v1 = *(node[1]) - *(node[0]);
332  Node v2 = *(node[2]) - *(node[0]);
333  UnitCrossProd(v1, v2, n);
334  }
335  else
336  {
337  // Calculate gradient vector and invert.
338  Node dx = *(node[1]) - *(node[0]);
339  NekDouble diff = dx.abs2();
340  dx /= sqrt(diff);
341  n.m_x = -dx.m_y;
342  n.m_y = dx.m_x;
343  n.m_z = 0;
344  }
345 
346  // Insert face normal into vertex normal list or add to existing
347  // value.
348  for (int j = 0; j < nV; ++j)
349  {
350  nIt = mesh->m_vertexNormals.find(e->GetVertex(j)->m_id);
351  if (nIt == mesh->m_vertexNormals.end())
352  {
353  mesh->m_vertexNormals[e->GetVertex(j)->m_id] = n;
354  }
355  else
356  {
357  nIt->second += n;
358  }
359  }
360  }
361 
362  // Normalize resulting vectors.
363  for (nIt = mesh->m_vertexNormals.begin();
364  nIt != mesh->m_vertexNormals.end();
365  ++nIt)
366  {
367  Node &n = mesh->m_vertexNormals[nIt->first];
368  n /= sqrt(n.abs2());
369  }
370 
371 }
void UnitCrossProd(NekMeshUtils::Node &a, NekMeshUtils::Node &b, NekMeshUtils::Node &c)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
NekDouble m_y
Y-coordinate.
Definition: Node.h:401
Represents a point in the domain.
Definition: Node.h:60
NEKMESHUTILS_EXPORT NekDouble abs2() const
Definition: Node.h:154
double NekDouble
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
NekDouble m_x
X-coordinate.
Definition: Node.h:399
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
NekDouble m_z
Z-coordinate.
Definition: Node.h:403
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.

References Nektar::NekMeshUtils::Node::abs2(), Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ASSERTL1, Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), CrossProdMag(), Nektar::NekMeshUtils::Node::dot(), Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauMAlpha1Beta0, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eNodalTriElec, Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eSegment, Nektar::LibUtilities::eTriangle, FindNormalFromPlyFile(), GenerateNormals(), Nektar::ParseUtils::GenerateSeqVector(), Nektar::ParseUtils::GenerateUnOrderedVector(), Nektar::NekMeshUtils::GetElementFactory(), Nektar::iterator, Nektar::NekMeshUtils::Module::m_config, Nektar::NekMeshUtils::Module::m_mesh, Nektar::NekMeshUtils::Node::m_x, Nektar::NekMeshUtils::Node::m_y, Nektar::NekMeshUtils::Node::m_z, CG_Iterations::Mesh, class_topology::Node, class_topology::P, SuperBlend(), and Vmath::Zero().

377 {
378  ASSERTL0(m_mesh->m_spaceDim == 3 || m_mesh->m_spaceDim == 2,
379  "Spherigon implementation only valid in 2D/3D.");
380 
382  boost::unordered_set<int> visitedEdges;
383 
384  // First construct vector of elements to process.
385  vector<ElementSharedPtr> el;
386 
387  if (m_mesh->m_verbose)
388  {
389  cout << "ProcessSpherigon: Smoothing mesh..." << endl;
390  }
391 
392  if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3)
393  {
394  // Manifold case - copy expansion dimension.
395  el = m_mesh->m_element[m_mesh->m_expDim];
396  }
397  else if (m_mesh->m_spaceDim == m_mesh->m_expDim)
398  {
399  // Full 2D or 3D case - iterate over stored edges/faces and
400  // create segments/triangles representing those edges/faces.
401  set<pair<int, int> >::iterator it;
402  vector<int> t;
403  t.push_back(0);
404 
405  // Construct list of spherigon edges/faces from a tag.
406  string surfTag = m_config["surf"].as<string>();
407  bool prismTag = m_config["BothTriFacesOnPrism"].beenSet;
408 
409  if (surfTag != "")
410  {
411  vector<unsigned int> surfs;
412  ParseUtils::GenerateSeqVector(surfTag.c_str(), surfs);
413  sort(surfs.begin(), surfs.end());
414 
415  m_mesh->m_spherigonSurfs.clear();
416  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
417  {
418  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
419  int nSurf = m_mesh->m_expDim == 3 ? el->GetFaceCount()
420  : el->GetEdgeCount();
421 
422  for (int j = 0; j < nSurf; ++j)
423  {
424  int bl = el->GetBoundaryLink(j);
425  if (bl == -1)
426  {
427  continue;
428  }
429 
430  ElementSharedPtr bEl =
431  m_mesh->m_element[m_mesh->m_expDim - 1][bl];
432  vector<int> tags = bEl->GetTagList();
433  vector<int> inter;
434 
435  sort(tags.begin(), tags.end());
436  set_intersection(surfs.begin(),
437  surfs.end(),
438  tags.begin(),
439  tags.end(),
440  back_inserter(inter));
441 
442  if (inter.size() == 1)
443  {
444  m_mesh->m_spherigonSurfs.insert(make_pair(i, j));
445 
446  // Curve other tri face on Prism. Note could be
447  // problem on pyramid when implemented.
448  if (nSurf == 5 && prismTag)
449  {
450  // add other end of prism on boundary for
451  // smoothing
452  int triFace = j == 1 ? 3 : 1;
453  m_mesh->m_spherigonSurfs.insert(
454  make_pair(i, triFace));
455  }
456  }
457  }
458  }
459  }
460 
461  if (m_mesh->m_spherigonSurfs.size() == 0)
462  {
463  cerr << "WARNING: Spherigon surfaces have not been defined "
464  << "-- ignoring smoothing." << endl;
465  return;
466  }
467 
468  if (m_mesh->m_expDim == 3)
469  {
470  for (it = m_mesh->m_spherigonSurfs.begin();
471  it != m_mesh->m_spherigonSurfs.end();
472  ++it)
473  {
474  FaceSharedPtr f =
475  m_mesh->m_element[m_mesh->m_expDim][it->first]->GetFace(
476  it->second);
477  vector<NodeSharedPtr> nodes = f->m_vertexList;
479  (LibUtilities::ShapeType)(nodes.size());
480  ElmtConfig conf(eType, 1, false, false);
481 
482  // Create 2D element.
483  ElementSharedPtr elmt =
484  GetElementFactory().CreateInstance(eType, conf, nodes, t);
485 
486  // Copy vertices/edges from face.
487  for (int i = 0; i < f->m_vertexList.size(); ++i)
488  {
489  elmt->SetVertex(i, f->m_vertexList[i]);
490  }
491  for (int i = 0; i < f->m_edgeList.size(); ++i)
492  {
493  elmt->SetEdge(i, f->m_edgeList[i]);
494  }
495 
496  el.push_back(elmt);
497  }
498  }
499  else
500  {
501  for (it = m_mesh->m_spherigonSurfs.begin();
502  it != m_mesh->m_spherigonSurfs.end();
503  ++it)
504  {
505  EdgeSharedPtr edge =
506  m_mesh->m_element[m_mesh->m_expDim][it->first]->GetEdge(
507  it->second);
508  vector<NodeSharedPtr> nodes;
510  ElmtConfig conf(eType, 1, false, false);
511 
512  nodes.push_back(edge->m_n1);
513  nodes.push_back(edge->m_n2);
514 
515  // Create 2D element.
516  ElementSharedPtr elmt =
517  GetElementFactory().CreateInstance(eType, conf, nodes, t);
518 
519  // Copy vertices/edges from original element.
520  elmt->SetVertex(0, nodes[0]);
521  elmt->SetVertex(1, nodes[1]);
522  elmt->SetEdge(
523  0,
524  m_mesh->m_element[m_mesh->m_expDim][it->first]->GetEdge(
525  it->second));
526  el.push_back(elmt);
527  }
528  }
529  }
530  else
531  {
532  ASSERTL0(false, "Spherigon expansions must be 2/3 dimensional");
533  }
534 
535  // See if vertex normals have been generated. If they have not,
536  // approximate them by summing normals of surrounding elements.
537  bool normalsGenerated = false;
538 
539  // Read Normal file if one exists
540  std::string normalfile = m_config["usenormalfile"].as<string>();
541  if (normalfile.compare("NoFile") != 0)
542  {
543  NekDouble scale = m_config["scalefile"].as<NekDouble>();
544 
545  if (m_mesh->m_verbose)
546  {
547  cout << "Inputing normal file: " << normalfile
548  << " with scaling of " << scale << endl;
549  }
550 
551  ifstream inplyTmp;
552  io::filtering_istream inply;
553  InputPlySharedPtr plyfile;
554 
555  inplyTmp.open(normalfile.c_str());
556  ASSERTL0(inplyTmp,
557  string("Could not open input ply file: ") + normalfile);
558 
559  inply.push(inplyTmp);
560 
561  MeshSharedPtr m = boost::shared_ptr<Mesh>(new Mesh());
562  plyfile = boost::shared_ptr<InputPly>(new InputPly(m));
563  plyfile->ReadPly(inply, scale);
564  plyfile->ProcessVertices();
565 
566  MeshSharedPtr plymesh = plyfile->GetMesh();
567  if (m_mesh->m_verbose)
568  {
569  cout << "\t Generating ply normals" << endl;
570  }
571  GenerateNormals(plymesh->m_element[plymesh->m_expDim], plymesh);
572 
573  // finally find nearest vertex and set normal to mesh surface file
574  // normal. probably should have a hex tree search ?
575  Node minx(0, 0.0, 0.0, 0.0), tmp, tmpsav;
578  map<int, NodeSharedPtr> surfverts;
579 
580  // make a map of normal vertices to visit based on elements el
581  for (int i = 0; i < el.size(); ++i)
582  {
583  ElementSharedPtr e = el[i];
584  int nV = e->GetVertexCount();
585  for (int j = 0; j < nV; ++j)
586  {
587  int id = e->GetVertex(j)->m_id;
588  surfverts[id] = e->GetVertex(j);
589  }
590  }
591 
592 
593  if (m_mesh->m_verbose)
594  {
595  cout << "\t Processing surface normals " << endl;
596  }
597 
598  // loop over all element in ply mesh and determine
599  // and set normal to nearest point in ply mesh
600  FindNormalFromPlyFile(plymesh,surfverts);
601 
602  normalsGenerated = true;
603  }
604  else if (m_mesh->m_vertexNormals.size() == 0)
605  {
606  GenerateNormals(el, m_mesh);
607  normalsGenerated = true;
608  }
609 
610  // See if we should add noise to normals
611  std::string normalnoise = m_config["normalnoise"].as<string>();
612  if (normalnoise.compare("NotSpecified") != 0)
613  {
614  vector<NekDouble> values;
615  ASSERTL0(
616  ParseUtils::GenerateUnOrderedVector(normalnoise.c_str(), values),
617  "Failed to interpret normal noise string");
618 
619  int nvalues = values.size() / 2;
620  NekDouble amp = values[0];
621 
622  if (m_mesh->m_verbose)
623  {
624  cout << "\t adding noise to normals of amplitude " << amp
625  << " in range: ";
626  for (int i = 0; i < nvalues; ++i)
627  {
628  cout << values[2 * i + 1] << "," << values[2 * i + 2] << " ";
629  }
630  cout << endl;
631  }
632 
634  map<int, NodeSharedPtr> surfverts;
635 
636  // make a map of normal vertices to visit based on elements el
637  for (int i = 0; i < el.size(); ++i)
638  {
639  ElementSharedPtr e = el[i];
640  int nV = e->GetVertexCount();
641  for (int j = 0; j < nV; ++j)
642  {
643  int id = e->GetVertex(j)->m_id;
644  surfverts[id] = e->GetVertex(j);
645  }
646  }
647 
648  for (vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt)
649  {
650  bool AddNoise = false;
651 
652  for (int i = 0; i < nvalues; ++i)
653  {
654  // check to see if point is in range
655  switch (nvalues)
656  {
657  case 1:
658  {
659  if (((vIt->second)->m_x > values[2 * i + 1]) &&
660  ((vIt->second)->m_x < values[2 * i + 2]))
661  {
662  AddNoise = true;
663  }
664  }
665  break;
666  case 2:
667  {
668  if (((vIt->second)->m_x > values[2 * i + 1]) &&
669  ((vIt->second)->m_x < values[2 * i + 2]) &&
670  ((vIt->second)->m_y > values[2 * i + 3]) &&
671  ((vIt->second)->m_y < values[2 * i + 4]))
672  {
673  AddNoise = true;
674  }
675  }
676  break;
677  case 3:
678  {
679  if (((vIt->second)->m_x > values[2 * i + 1]) &&
680  ((vIt->second)->m_x < values[2 * i + 2]) &&
681  ((vIt->second)->m_y > values[2 * i + 3]) &&
682  ((vIt->second)->m_y < values[2 * i + 4]) &&
683  ((vIt->second)->m_z > values[2 * i + 5]) &&
684  ((vIt->second)->m_z < values[2 * i + 6]))
685 
686  {
687  AddNoise = true;
688  }
689  }
690  break;
691  }
692 
693  if (AddNoise)
694  {
695  // generate random unit vector;
696  Node rvec(0, rand(), rand(), rand());
697  rvec *= values[0] / sqrt(rvec.abs2());
698 
699  Node normal = m_mesh->m_vertexNormals[vIt->first];
700 
701  normal += rvec;
702  normal /= sqrt(normal.abs2());
703 
704  m_mesh->m_vertexNormals[vIt->first] = normal;
705  }
706  }
707  }
708  }
709 
710  // Allocate storage for interior points.
711  int nq = m_config["N"].as<int>();
712  int nquad = m_mesh->m_spaceDim == 3 ? nq * nq : nq;
713  Array<OneD, NekDouble> x(nq * nq);
714  Array<OneD, NekDouble> y(nq * nq);
715  Array<OneD, NekDouble> z(nq * nq);
716 
717  Array<OneD, NekDouble> xc(nq * nq);
718  Array<OneD, NekDouble> yc(nq * nq);
719  Array<OneD, NekDouble> zc(nq * nq);
720 
721  ASSERTL0(nq > 2, "Number of points must be greater than 2.");
722 
723  LibUtilities::BasisKey B0(
725  nq,
726  LibUtilities::PointsKey(nq, LibUtilities::eGaussLobattoLegendre));
727  LibUtilities::BasisKey B1(
729  nq,
730  LibUtilities::PointsKey(nq, LibUtilities::eGaussRadauMAlpha1Beta0));
734 
735  Array<OneD, NekDouble> xnodal(nq * (nq + 1) / 2), ynodal(nq * (nq + 1) / 2);
736  stdtri->GetNodalPoints(xnodal, ynodal);
737 
738  int edgeMap[3][4][2] = {
739  {{0, 1}, {-1, -1}, {-1, -1}, {-1, -1}}, // seg
740  {{0, 1}, {nq - 1, nq}, {nq * (nq - 1), -nq}, {-1, -1}}, // tri
741  {{0, 1}, {nq - 1, nq}, {nq * nq - 1, -1}, {nq * (nq - 1), -nq}} // quad
742  };
743 
744  int vertMap[3][4][2] = {
745  {{0, 1}, {0, 0}, {0, 0}, {0, 0}}, // seg
746  {{0, 1}, {1, 2}, {2, 3}, {0, 0}}, // tri
747  {{0, 1}, {1, 2}, {2, 3}, {3, 0}}, // quad
748  };
749 
750  for (int i = 0; i < el.size(); ++i)
751  {
752  // Construct a Nektar++ element to obtain coordinate points
753  // inside the element. TODO: Add options for various
754  // nodal/tensor point distributions + number of points to add.
755  ElementSharedPtr e = el[i];
756 
757  LibUtilities::BasisKey B2(
759  nq,
760  LibUtilities::PointsKey(nq, LibUtilities::eGaussLobattoLegendre));
761 
762  if (e->GetConf().m_e == LibUtilities::eSegment)
763  {
765  boost::dynamic_pointer_cast<SpatialDomains::SegGeom>(
766  e->GetGeom(m_mesh->m_spaceDim));
769  geom);
770  seg->GetCoords(x, y, z);
771  nquad = nq;
772  }
773  else if (e->GetConf().m_e == LibUtilities::eTriangle)
774  {
776  boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
777  e->GetGeom(3));
780  B0, B1, LibUtilities::eNodalTriElec, geom);
781 
782  Array<OneD, NekDouble> coord(2);
783  tri->GetCoords(xc, yc, zc);
784  nquad = nq * (nq + 1) / 2;
785 
786  for (int j = 0; j < nquad; ++j)
787  {
788  coord[0] = xnodal[j];
789  coord[1] = ynodal[j];
790  x[j] = stdtri->PhysEvaluate(coord, xc);
791  }
792 
793  for (int j = 0; j < nquad; ++j)
794  {
795  coord[0] = xnodal[j];
796  coord[1] = ynodal[j];
797  y[j] = stdtri->PhysEvaluate(coord, yc);
798  }
799 
800  for (int j = 0; j < nquad; ++j)
801  {
802  coord[0] = xnodal[j];
803  coord[1] = ynodal[j];
804  z[j] = stdtri->PhysEvaluate(coord, zc);
805  }
806  }
807  else if (e->GetConf().m_e == LibUtilities::eQuadrilateral)
808  {
810  boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
811  e->GetGeom(3));
814  B2, B2, geom);
815  quad->GetCoords(x, y, z);
816  nquad = nq * nq;
817  }
818  else
819  {
820  ASSERTL0(false, "Unknown expansion type.");
821  }
822 
823  // Zero z-coordinate in 2D.
824  if (m_mesh->m_spaceDim == 2)
825  {
826  Vmath::Zero(nquad, z, 1);
827  }
828 
829  // Find vertex normals.
830  int nV = e->GetVertexCount();
831  vector<Node> v, vN;
832  for (int j = 0; j < nV; ++j)
833  {
834  v.push_back(*(e->GetVertex(j)));
835  ASSERTL1(m_mesh->m_vertexNormals.count(v[j].m_id) != 0,
836  "Normal has not been defined");
837  vN.push_back(m_mesh->m_vertexNormals[v[j].m_id]);
838  }
839 
840  vector<Node> tmp(nV);
841  vector<NekDouble> r(nV);
842  vector<Node> K(nV);
843  vector<Node> Q(nV);
844  vector<Node> Qp(nV);
845  vector<NekDouble> blend(nV);
846  vector<Node> out(nquad);
847 
848  // Calculate segment length for 2D spherigon routine.
849  NekDouble segLength = sqrt((v[0] - v[1]).abs2());
850 
851  // Perform Spherigon method to smooth manifold.
852  for (int j = 0; j < nquad; ++j)
853  {
854  Node P(0, x[j], y[j], z[j]);
855  Node N(0, 0, 0, 0);
856 
857  // Calculate generalised barycentric coordinates r[] and the
858  // Phong normal N = vN . r for this point of the element.
859  if (m_mesh->m_spaceDim == 2)
860  {
861  // In 2D the coordinates are given by a ratio of the
862  // segment length to the distance from one of the
863  // endpoints.
864  r[0] = sqrt((P - v[0]).abs2()) / segLength;
865  r[0] = max(min(1.0, r[0]), 0.0);
866  r[1] = 1.0 - r[0];
867 
868  // Calculate Phong normal.
869  N = vN[0] * r[0] + vN[1] * r[1];
870  }
871  else if (m_mesh->m_spaceDim == 3)
872  {
873  for (int k = 0; k < nV; ++k)
874  {
875  tmp[k] = P - v[k];
876  }
877 
878  // Calculate generalized barycentric coordinate system
879  // (see equation 6 of paper).
880  NekDouble weight = 0.0;
881  for (int k = 0; k < nV; ++k)
882  {
883  r[k] = 1.0;
884  for (int l = 0; l < nV - 2; ++l)
885  {
886  r[k] *= CrossProdMag(tmp[(k + l + 1) % nV],
887  tmp[(k + l + 2) % nV]);
888  }
889  weight += r[k];
890  }
891 
892  // Calculate Phong normal (equation 1).
893  for (int k = 0; k < nV; ++k)
894  {
895  r[k] /= weight;
896  N += vN[k] * r[k];
897  }
898  }
899 
900  // Normalise Phong normal.
901  N /= sqrt(N.abs2());
902 
903  for (int k = 0; k < nV; ++k)
904  {
905  // Perform steps denoted in equations 2, 3, 8 for C1
906  // smoothing.
907  NekDouble tmp1;
908  K[k] = P + N * ((v[k] - P).dot(N));
909  tmp1 = (v[k] - K[k]).dot(vN[k]) / (1.0 + N.dot(vN[k]));
910  Q[k] = K[k] + N * tmp1;
911  Qp[k] = v[k] - N * ((v[k] - P).dot(N));
912  }
913 
914  // Apply C1 blending function to the surface. TODO: Add
915  // option to do (more efficient) C0 blending function.
916  SuperBlend(r, Qp, P, blend);
917  P.m_x = P.m_y = P.m_z = 0.0;
918 
919  // Apply blending (equation 4).
920  for (int k = 0; k < nV; ++k)
921  {
922  P += Q[k] * blend[k];
923  }
924 
925  if ((boost::math::isnan)(P.m_x) || (boost::math::isnan)(P.m_y) ||
926  (boost::math::isnan)(P.m_z))
927  {
928  ASSERTL0(false,
929  "spherigon point is a nan. Check to see if "
930  "ply file is correct if using input normal file");
931  }
932  else
933  {
934  out[j] = P;
935  }
936  }
937 
938  // Push nodes into lines - TODO: face interior nodes.
939  // offset = 0 (seg), 1 (tri) or 2 (quad)
940  int offset = (int)e->GetConf().m_e - 2;
941 
942  for (int edge = 0; edge < e->GetEdgeCount(); ++edge)
943  {
944  eIt = visitedEdges.find(e->GetEdge(edge)->m_id);
945  if (eIt == visitedEdges.end())
946  {
947  bool reverseEdge =
948  !(v[vertMap[offset][edge][0]] == *(e->GetEdge(edge)->m_n1));
949 
950  // Clear existing curvature.
951  e->GetEdge(edge)->m_edgeNodes.clear();
952 
953  if (e->GetConf().m_e != LibUtilities::eTriangle)
954  {
955  for (int j = 1; j < nq - 1; ++j)
956  {
957  int v = edgeMap[offset][edge][0] +
958  j * edgeMap[offset][edge][1];
959  e->GetEdge(edge)->m_edgeNodes.push_back(
960  NodeSharedPtr(new Node(out[v])));
961  }
962  }
963  else
964  {
965  for (int j = 0; j < nq - 2; ++j)
966  {
967  int v = 3 + edge * (nq - 2) + j;
968  e->GetEdge(edge)->m_edgeNodes.push_back(
969  NodeSharedPtr(new Node(out[v])));
970  }
971  }
972 
973  if (reverseEdge)
974  {
975  reverse(e->GetEdge(edge)->m_edgeNodes.begin(),
976  e->GetEdge(edge)->m_edgeNodes.end());
977  }
978 
979  e->GetEdge(edge)->m_curveType =
981 
982  visitedEdges.insert(e->GetEdge(edge)->m_id);
983  }
984  }
985 
986  // Add face nodes in manifold and full 3D case, but not for 2D.
987  if (m_mesh->m_spaceDim == 3)
988  {
989  vector<NodeSharedPtr> volNodes;
990 
991  if (e->GetConf().m_e == LibUtilities::eQuadrilateral)
992  {
993  volNodes.resize((nq - 2) * (nq - 2));
994  for (int j = 1; j < nq - 1; ++j)
995  {
996  for (int k = 1; k < nq - 1; ++k)
997  {
998  int v = j * nq + k;
999  volNodes[(j - 1) * (nq - 2) + (k - 1)] =
1000  NodeSharedPtr(new Node(out[v]));
1001  }
1002  }
1003  }
1004  else
1005  {
1006  for (int j = 3 + 3 * (nq - 2); j < nquad; ++j)
1007  {
1008  volNodes.push_back(NodeSharedPtr(new Node(out[j])));
1009  }
1010  }
1011 
1012  e->SetVolumeNodes(volNodes);
1013  }
1014  }
1015 
1016  // Copy face nodes back into 3D element faces
1017  if (m_mesh->m_expDim == 3)
1018  {
1019  set<pair<int, int> >::iterator it;
1020  int elmt = 0;
1021  for (it = m_mesh->m_spherigonSurfs.begin();
1022  it != m_mesh->m_spherigonSurfs.end();
1023  ++it, ++elmt)
1024  {
1025  FaceSharedPtr f =
1026  m_mesh->m_element[m_mesh->m_expDim][it->first]->GetFace(
1027  it->second);
1028 
1029  f->m_faceNodes = el[elmt]->GetVolumeNodes();
1030  f->m_curveType = f->m_vertexList.size() == 3
1033  }
1034  }
1035 
1036  if (normalsGenerated)
1037  {
1038  m_mesh->m_vertexNormals.clear();
1039  }
1040 }
NekDouble CrossProdMag(NekMeshUtils::Node &a, NekMeshUtils::Node &b)
Calculate the magnitude of the cross product .
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Basic information about an element.
Definition: ElementConfig.h:50
void FindNormalFromPlyFile(NekMeshUtils::MeshSharedPtr &plymesh, std::map< int, NekMeshUtils::NodeSharedPtr > &surfverts)
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
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.
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 hy...
Principle Modified Functions .
Definition: BasisType.h:49
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
boost::shared_ptr< StdNodalTriExp > StdNodalTriExpSharedPtr
Represents a point in the domain.
Definition: Node.h:60
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
Principle Orthogonal Functions .
Definition: BasisType.h:47
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:270
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
Principle Orthogonal Functions .
Definition: BasisType.h:46
boost::shared_ptr< QuadExp > QuadExpSharedPtr
Definition: QuadExp.h:297
double NekDouble
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< InputPly > InputPlySharedPtr
Definition: InputPly.h:68
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:136
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
boost::shared_ptr< Face > FaceSharedPtr
Definition: Face.h:153
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
static bool GenerateUnOrderedVector(const char *const str, std::vector< NekDouble > &vec)
Definition: ParseUtils.hpp:128
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:70
boost::shared_ptr< NodalTriExp > NodalTriExpSharedPtr
Definition: NodalTriExp.h:381
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
rGeneralised barycentric coordinates of the point P.
QVector of vertices denoting this triangle/quad.
PPoint in the triangle to apply blending to.
blendThe resulting blending components for each vertex.

Definition at line 181 of file ProcessSpherigon.cpp.

References class_topology::P, and TOL_BLEND.

Referenced by Process().

185 {
186  vector<double> tmp(r.size());
187  double totBlend = 0.0;
188  int i;
189  int nV = r.size();
190 
191  for (i = 0; i < nV; ++i)
192  {
193  blend[i] = 0.0;
194  tmp[i] = (Q[i] - P).abs2();
195  }
196 
197  for (i = 0; i < nV; ++i)
198  {
199  int ip = (i + 1) % nV, im = (i - 1 + nV) % nV;
200 
201  if (r[i] > TOL_BLEND && r[i] < (1 - TOL_BLEND))
202  {
203  blend[i] =
204  r[i] * r[i] * (r[im] * r[im] * tmp[im] / (tmp[im] + tmp[i]) +
205  r[ip] * r[ip] * tmp[ip] / (tmp[ip] + tmp[i]));
206  totBlend += blend[i];
207  }
208  }
209 
210  for (i = 0; i < nV; ++i)
211  {
212  blend[i] /= totBlend;
213  if (r[i] >= (1 - TOL_BLEND))
214  {
215  blend[i] = 1.0;
216  }
217  if (r[i] <= TOL_BLEND)
218  {
219  blend[i] = 0.0;
220  }
221  }
222 }
#define TOL_BLEND
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().

131 {
132  double inv_mag;
133 
134  c.m_x = a.m_y * b.m_z - a.m_z * b.m_y;
135  c.m_y = a.m_z * b.m_x - a.m_x * b.m_z;
136  c.m_z = a.m_x * b.m_y - a.m_y * b.m_x;
137 
138  inv_mag = 1.0 / sqrt(c.m_x * c.m_x + c.m_y * c.m_y + c.m_z * c.m_z);
139 
140  c.m_x = c.m_x * inv_mag;
141  c.m_y = c.m_y * inv_mag;
142  c.m_z = c.m_z * inv_mag;
143 }

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessSpherigon::className
static
Initial value:

Definition at line 56 of file ProcessSpherigon.h.