Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 (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 (std::string key, std::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< Modulecreate (MeshSharedPtr m)
 Creates an instance of this class. More...
 

Static Public Attributes

static ModuleKey className
 

Protected Member Functions

void GenerateNormals (std::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 $ C^0 $ blending function for spherigon implementation. More...
 
void SuperBlend (std::vector< double > &r, std::vector< Node > &Q, Node &P, std::vector< double > &blend)
 Calculate the $ C^1 $ 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, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
 

Additional Inherited Members

- Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr m_f
 Field object. More...
 
map< string, ConfigOptionm_config
 List of configuration values. More...
 
bool m_requireEquiSpaced
 
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 ( MeshSharedPtr  m)

Default constructor.

Definition at line 91 of file ProcessSpherigon.cpp.

References Nektar::Utilities::Module::m_config.

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

Definition at line 157 of file ProcessSpherigon.cpp.

158 {
159  return r * r;
160 }
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().

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

142 {
143  double tmp1 = a.m_y * b.m_z - a.m_z * b.m_y;
144  double tmp2 = a.m_z * b.m_x - a.m_x * b.m_z;
145  double tmp3 = a.m_x * b.m_y - a.m_y * b.m_x;
146  return sqrt(tmp1 * tmp1 + tmp2 * tmp2 + tmp3 * tmp3);
147 }
NekDouble m_y
Y-coordinate.
Definition: Node.h:314
NekDouble m_x
X-coordinate.
Definition: Node.h:312
NekDouble m_z
Z-coordinate.
Definition: Node.h:316
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
elVector 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().

231 {
233 
234  for (int i = 0; i < el.size(); ++i)
235  {
236  ElementSharedPtr e = el[i];
237 
238  // Ensure that element is a line, triangle or quad.
239  ASSERTL0(e->GetConf().m_e == LibUtilities::eSegment ||
240  e->GetConf().m_e == LibUtilities::eTriangle ||
241  e->GetConf().m_e == LibUtilities::eQuadrilateral,
242  "Spherigon expansions must be lines, triangles or "
243  "quadrilaterals.");
244 
245  // Calculate normal for this element.
246  int nV = e->GetVertexCount();
247  vector<NodeSharedPtr> node(nV);
248 
249  for (int j = 0; j < nV; ++j)
250  {
251  node[j] = e->GetVertex(j);
252  }
253 
254  Node n;
255 
256  if (mesh->m_spaceDim == 3)
257  {
258  // Create two tangent vectors and take unit cross product.
259  Node v1 = *(node[1]) - *(node[0]);
260  Node v2 = *(node[2]) - *(node[0]);
261  UnitCrossProd(v1, v2, n);
262  }
263  else
264  {
265  // Calculate gradient vector and invert.
266  Node dx = *(node[1]) - *(node[0]);
267  dx /= sqrt(dx.abs2());
268  n.m_x = -dx.m_y;
269  n.m_y = dx.m_x;
270  n.m_z = 0;
271  }
272 
273  // Insert face normal into vertex normal list or add to existing
274  // value.
275  for (int j = 0; j < nV; ++j)
276  {
277  nIt = mesh->m_vertexNormals.find(e->GetVertex(j)->m_id);
278  if (nIt == mesh->m_vertexNormals.end())
279  {
280  mesh->m_vertexNormals[e->GetVertex(j)->m_id] = n;
281  }
282  else
283  {
284  nIt->second += n;
285  }
286  }
287  }
288 
289  // Normalize resulting vectors.
290  for (nIt = mesh->m_vertexNormals.begin();
291  nIt != mesh->m_vertexNormals.end();
292  ++nIt)
293  {
294  Node &n = mesh->m_vertexNormals[nIt->first];
295  n /= sqrt(n.abs2());
296  }
297 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
NekDouble m_y
Y-coordinate.
Definition: Node.h:314
Represents a point in the domain.
Definition: Node.h:60
NEKMESHUTILS_EXPORT NekDouble abs2() const
Definition: Node.h:148
void UnitCrossProd(Node &a, Node &b, Node &c)
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
NekDouble m_x
X-coordinate.
Definition: Node.h:312
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
NekDouble m_z
Z-coordinate.
Definition: Node.h:316
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.

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, GenerateNormals(), Nektar::ParseUtils::GenerateSeqVector(), Nektar::ParseUtils::GenerateUnOrderedVector(), Nektar::NekMeshUtils::GetElementFactory(), Nektar::iterator, Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_mesh, Nektar::NekMeshUtils::Node::m_x, Nektar::NekMeshUtils::Node::m_y, Nektar::NekMeshUtils::Node::m_z, SuperBlend(), and Vmath::Zero().

303 {
304  ASSERTL0(m_mesh->m_spaceDim == 3 || m_mesh->m_spaceDim == 2,
305  "Spherigon implementation only valid in 2D/3D.");
306 
308  boost::unordered_set<int> visitedEdges;
309 
310  // First construct vector of elements to process.
311  vector<ElementSharedPtr> el;
312 
313  if (m_mesh->m_verbose)
314  {
315  cout << "ProcessSpherigon: Smoothing mesh..." << endl;
316  }
317 
318  if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3)
319  {
320  // Manifold case - copy expansion dimension.
321  el = m_mesh->m_element[m_mesh->m_expDim];
322  }
323  else if (m_mesh->m_spaceDim == m_mesh->m_expDim)
324  {
325  // Full 2D or 3D case - iterate over stored edges/faces and
326  // create segments/triangles representing those edges/faces.
327  set<pair<int, int> >::iterator it;
328  vector<int> t;
329  t.push_back(0);
330 
331  // Construct list of spherigon edges/faces from a tag.
332  string surfTag = m_config["surf"].as<string>();
333  bool prismTag = m_config["BothTriFacesOnPrism"].beenSet;
334 
335  if (surfTag != "")
336  {
337  vector<unsigned int> surfs;
338  ParseUtils::GenerateSeqVector(surfTag.c_str(), surfs);
339  sort(surfs.begin(), surfs.end());
340 
341  m_mesh->m_spherigonSurfs.clear();
342  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
343  {
344  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
345  int nSurf = m_mesh->m_expDim == 3 ? el->GetFaceCount()
346  : el->GetEdgeCount();
347 
348  for (int j = 0; j < nSurf; ++j)
349  {
350  int bl = el->GetBoundaryLink(j);
351  if (bl == -1)
352  {
353  continue;
354  }
355 
356  ElementSharedPtr bEl =
357  m_mesh->m_element[m_mesh->m_expDim - 1][bl];
358  vector<int> tags = bEl->GetTagList();
359  vector<int> inter;
360 
361  sort(tags.begin(), tags.end());
362  set_intersection(surfs.begin(),
363  surfs.end(),
364  tags.begin(),
365  tags.end(),
366  back_inserter(inter));
367 
368  if (inter.size() == 1)
369  {
370  m_mesh->m_spherigonSurfs.insert(make_pair(i, j));
371 
372  // Curve other tri face on Prism. Note could be
373  // problem on pyramid when implemented.
374  if (nSurf == 5 && prismTag)
375  {
376  // add other end of prism on boundary for
377  // smoothing
378  int triFace = j == 1 ? 3 : 1;
379  m_mesh->m_spherigonSurfs.insert(
380  make_pair(i, triFace));
381  }
382  }
383  }
384  }
385  }
386 
387  if (m_mesh->m_spherigonSurfs.size() == 0)
388  {
389  cerr << "WARNING: Spherigon surfaces have not been defined "
390  << "-- ignoring smoothing." << endl;
391  return;
392  }
393 
394  if (m_mesh->m_expDim == 3)
395  {
396  for (it = m_mesh->m_spherigonSurfs.begin();
397  it != m_mesh->m_spherigonSurfs.end();
398  ++it)
399  {
400  FaceSharedPtr f =
401  m_mesh->m_element[m_mesh->m_expDim][it->first]->GetFace(
402  it->second);
403  vector<NodeSharedPtr> nodes = f->m_vertexList;
405  (LibUtilities::ShapeType)(nodes.size());
406  ElmtConfig conf(eType, 1, false, false);
407 
408  // Create 2D element.
409  ElementSharedPtr elmt =
410  GetElementFactory().CreateInstance(eType, conf, nodes, t);
411 
412  // Copy vertices/edges from face.
413  for (int i = 0; i < f->m_vertexList.size(); ++i)
414  {
415  elmt->SetVertex(i, f->m_vertexList[i]);
416  }
417  for (int i = 0; i < f->m_edgeList.size(); ++i)
418  {
419  elmt->SetEdge(i, f->m_edgeList[i]);
420  }
421 
422  el.push_back(elmt);
423  }
424  }
425  else
426  {
427  for (it = m_mesh->m_spherigonSurfs.begin();
428  it != m_mesh->m_spherigonSurfs.end();
429  ++it)
430  {
431  EdgeSharedPtr edge =
432  m_mesh->m_element[m_mesh->m_expDim][it->first]->GetEdge(
433  it->second);
434  vector<NodeSharedPtr> nodes;
436  ElmtConfig conf(eType, 1, false, false);
437 
438  nodes.push_back(edge->m_n1);
439  nodes.push_back(edge->m_n2);
440 
441  // Create 2D element.
442  ElementSharedPtr elmt =
443  GetElementFactory().CreateInstance(eType, conf, nodes, t);
444 
445  // Copy vertices/edges from original element.
446  elmt->SetVertex(0, nodes[0]);
447  elmt->SetVertex(1, nodes[1]);
448  elmt->SetEdge(
449  0,
450  m_mesh->m_element[m_mesh->m_expDim][it->first]->GetEdge(
451  it->second));
452  el.push_back(elmt);
453  }
454  }
455  }
456  else
457  {
458  ASSERTL0(false, "Spherigon expansions must be 2/3 dimensional");
459  }
460 
461  // See if vertex normals have been generated. If they have not,
462  // approximate them by summing normals of surrounding elements.
463  bool normalsGenerated = false;
464 
465  // Read Normal file if one exists
466  std::string normalfile = m_config["usenormalfile"].as<string>();
467  if (normalfile.compare("NoFile") != 0)
468  {
469  NekDouble scale = m_config["scalefile"].as<NekDouble>();
470 
471  if (m_mesh->m_verbose)
472  {
473  cout << "Inputing normal file: " << normalfile
474  << " with scaling of " << scale << endl;
475  }
476 
477  ifstream inply;
478  InputPlySharedPtr plyfile;
479 
480  inply.open(normalfile.c_str());
481  ASSERTL0(inply, string("Could not open input ply file: ") + normalfile);
482 
483  int j;
484  MeshSharedPtr m = boost::shared_ptr<Mesh>(new Mesh());
485  plyfile = boost::shared_ptr<InputPly>(new InputPly(m));
486  plyfile->ReadPly(inply, scale);
487  plyfile->ProcessVertices();
488 
489  MeshSharedPtr plymesh = plyfile->GetMesh();
490  GenerateNormals(plymesh->m_element[plymesh->m_expDim], plymesh);
491 
492  // finally find nearest vertex and set normal to mesh surface file
493  // normal. probably should have a hex tree search ?
494  Array<OneD, NekDouble> len2(plymesh->m_vertexSet.size());
495  Node minx(0, 0.0, 0.0, 0.0), tmp, tmpsav;
498  map<int, NodeSharedPtr> surfverts;
499 
500  // make a map of normal vertices to visit based on elements el
501  for (int i = 0; i < el.size(); ++i)
502  {
503  ElementSharedPtr e = el[i];
504  int nV = e->GetVertexCount();
505  for (int j = 0; j < nV; ++j)
506  {
507  int id = e->GetVertex(j)->m_id;
508  surfverts[id] = e->GetVertex(j);
509  }
510  }
511 
512  // loop over all element in ply mesh and determine
513  // xmin,xmax,ymin,ymax as search criterion
514 
515  NekDouble mindiff, diff;
516  int cntmin;
517 
518  if (m_mesh->m_verbose)
519  {
520  cout << "\t Processing surface normals " << endl;
521  }
522  int cnt = 0;
523  map<int, int> locnorm;
524  for (vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt, ++cnt)
525  {
526  mindiff = 1e12;
527 
528  for (j = 0, it = plymesh->m_vertexSet.begin();
529  it != plymesh->m_vertexSet.end();
530  ++it, ++j)
531  {
532  tmp = *(vIt->second) - *(*it);
533  diff = tmp.abs2();
534 
535  if (diff < mindiff)
536  {
537  mindiff = diff;
538  cntmin = (*it)->m_id;
539  tmpsav = tmp;
540  }
541  }
542  locnorm[cntmin] = vIt->first;
543 
544  ASSERTL1(cntmin < plymesh->m_vertexNormals.size(),
545  "cntmin is out of range");
546  m_mesh->m_vertexNormals[vIt->first] =
547  plymesh->m_vertexNormals[cntmin];
548  }
549  if (m_mesh->m_verbose)
550  {
551  cout << "\t end of processing surface normals " << endl;
552  }
553  normalsGenerated = true;
554  }
555  else if (m_mesh->m_vertexNormals.size() == 0)
556  {
557  GenerateNormals(el, m_mesh);
558  normalsGenerated = true;
559  }
560 
561  // See if we should add noise to normals
562  std::string normalnoise = m_config["normalnoise"].as<string>();
563  if (normalnoise.compare("NotSpecified") != 0)
564  {
565  vector<NekDouble> values;
566  ASSERTL0(
567  ParseUtils::GenerateUnOrderedVector(normalnoise.c_str(), values),
568  "Failed to interpret normal noise string");
569 
570  int nvalues = values.size() / 2;
571  NekDouble amp = values[0];
572 
573  if (m_mesh->m_verbose)
574  {
575  cout << "\t adding noise to normals of amplitude " << amp
576  << " in range: ";
577  for (int i = 0; i < nvalues; ++i)
578  {
579  cout << values[2 * i + 1] << "," << values[2 * i + 2] << " ";
580  }
581  cout << endl;
582  }
583 
585  map<int, NodeSharedPtr> surfverts;
586 
587  // make a map of normal vertices to visit based on elements el
588  for (int i = 0; i < el.size(); ++i)
589  {
590  ElementSharedPtr e = el[i];
591  int nV = e->GetVertexCount();
592  for (int j = 0; j < nV; ++j)
593  {
594  int id = e->GetVertex(j)->m_id;
595  surfverts[id] = e->GetVertex(j);
596  }
597  }
598 
599  for (vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt)
600  {
601  bool AddNoise = false;
602 
603  for (int i = 0; i < nvalues; ++i)
604  {
605  // check to see if point is in range
606  switch (nvalues)
607  {
608  case 1:
609  {
610  if (((vIt->second)->m_x > values[2 * i + 1]) &&
611  ((vIt->second)->m_x < values[2 * i + 2]))
612  {
613  AddNoise = true;
614  }
615  }
616  break;
617  case 2:
618  {
619  if (((vIt->second)->m_x > values[2 * i + 1]) &&
620  ((vIt->second)->m_x < values[2 * i + 2]) &&
621  ((vIt->second)->m_y > values[2 * i + 3]) &&
622  ((vIt->second)->m_y < values[2 * i + 4]))
623  {
624  AddNoise = true;
625  }
626  }
627  break;
628  case 3:
629  {
630  if (((vIt->second)->m_x > values[2 * i + 1]) &&
631  ((vIt->second)->m_x < values[2 * i + 2]) &&
632  ((vIt->second)->m_y > values[2 * i + 3]) &&
633  ((vIt->second)->m_y < values[2 * i + 4]) &&
634  ((vIt->second)->m_z > values[2 * i + 5]) &&
635  ((vIt->second)->m_z < values[2 * i + 6]))
636 
637  {
638  AddNoise = true;
639  }
640  }
641  break;
642  }
643 
644  if (AddNoise)
645  {
646  // generate random unit vector;
647  Node rvec(0, rand(), rand(), rand());
648  rvec *= values[0] / sqrt(rvec.abs2());
649 
650  Node normal = m_mesh->m_vertexNormals[vIt->first];
651 
652  normal += rvec;
653  normal /= sqrt(normal.abs2());
654 
655  m_mesh->m_vertexNormals[vIt->first] = normal;
656  }
657  }
658  }
659  }
660 
661  // Allocate storage for interior points.
662  int nq = m_config["N"].as<int>();
663  int nquad = m_mesh->m_spaceDim == 3 ? nq * nq : nq;
664  Array<OneD, NekDouble> x(nq * nq);
665  Array<OneD, NekDouble> y(nq * nq);
666  Array<OneD, NekDouble> z(nq * nq);
667 
668  Array<OneD, NekDouble> xc(nq * nq);
669  Array<OneD, NekDouble> yc(nq * nq);
670  Array<OneD, NekDouble> zc(nq * nq);
671 
672  ASSERTL0(nq > 2, "Number of points must be greater than 2.");
673 
674  LibUtilities::BasisKey B0(
676  nq,
677  LibUtilities::PointsKey(nq, LibUtilities::eGaussLobattoLegendre));
678  LibUtilities::BasisKey B1(
680  nq,
681  LibUtilities::PointsKey(nq, LibUtilities::eGaussRadauMAlpha1Beta0));
685 
686  Array<OneD, NekDouble> xnodal(nq * (nq + 1) / 2), ynodal(nq * (nq + 1) / 2);
687  stdtri->GetNodalPoints(xnodal, ynodal);
688 
689  int edgeMap[3][4][2] = {
690  {{0, 1}, {-1, -1}, {-1, -1}, {-1, -1}}, // seg
691  {{0, 1}, {nq - 1, nq}, {nq * (nq - 1), -nq}, {-1, -1}}, // tri
692  {{0, 1}, {nq - 1, nq}, {nq * nq - 1, -1}, {nq * (nq - 1), -nq}} // quad
693  };
694 
695  int vertMap[3][4][2] = {
696  {{0, 1}, {0, 0}, {0, 0}, {0, 0}}, // seg
697  {{0, 1}, {1, 2}, {2, 3}, {0, 0}}, // tri
698  {{0, 1}, {1, 2}, {2, 3}, {3, 0}}, // quad
699  };
700 
701  for (int i = 0; i < el.size(); ++i)
702  {
703  // Construct a Nektar++ element to obtain coordinate points
704  // inside the element. TODO: Add options for various
705  // nodal/tensor point distributions + number of points to add.
706  ElementSharedPtr e = el[i];
707 
708  LibUtilities::BasisKey B2(
710  nq,
711  LibUtilities::PointsKey(nq, LibUtilities::eGaussLobattoLegendre));
712 
713  if (e->GetConf().m_e == LibUtilities::eSegment)
714  {
716  boost::dynamic_pointer_cast<SpatialDomains::SegGeom>(
717  e->GetGeom(m_mesh->m_spaceDim));
720  geom);
721  seg->GetCoords(x, y, z);
722  nquad = nq;
723  }
724  else if (e->GetConf().m_e == LibUtilities::eTriangle)
725  {
727  boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
728  e->GetGeom(3));
731  B0, B1, LibUtilities::eNodalTriElec, geom);
732 
733  Array<OneD, NekDouble> coord(2);
734  tri->GetCoords(xc, yc, zc);
735  nquad = nq * (nq + 1) / 2;
736 
737  for (int j = 0; j < nquad; ++j)
738  {
739  coord[0] = xnodal[j];
740  coord[1] = ynodal[j];
741  x[j] = stdtri->PhysEvaluate(coord, xc);
742  }
743 
744  for (int j = 0; j < nquad; ++j)
745  {
746  coord[0] = xnodal[j];
747  coord[1] = ynodal[j];
748  y[j] = stdtri->PhysEvaluate(coord, yc);
749  }
750 
751  for (int j = 0; j < nquad; ++j)
752  {
753  coord[0] = xnodal[j];
754  coord[1] = ynodal[j];
755  z[j] = stdtri->PhysEvaluate(coord, zc);
756  }
757  }
758  else if (e->GetConf().m_e == LibUtilities::eQuadrilateral)
759  {
761  boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
762  e->GetGeom(3));
765  B2, B2, geom);
766  quad->GetCoords(x, y, z);
767  nquad = nq * nq;
768  }
769  else
770  {
771  ASSERTL0(false, "Unknown expansion type.");
772  }
773 
774  // Zero z-coordinate in 2D.
775  if (m_mesh->m_spaceDim == 2)
776  {
777  Vmath::Zero(nquad, z, 1);
778  }
779 
780  // Find vertex normals.
781  int nV = e->GetVertexCount();
782  vector<Node> v, vN;
783  for (int j = 0; j < nV; ++j)
784  {
785  v.push_back(*(e->GetVertex(j)));
786  ASSERTL1(m_mesh->m_vertexNormals.count(v[j].m_id) != 0,
787  "Normal has not been defined");
788  vN.push_back(m_mesh->m_vertexNormals[v[j].m_id]);
789  }
790 
791  vector<Node> tmp(nV);
792  vector<double> r(nV);
793  vector<Node> K(nV);
794  vector<Node> Q(nV);
795  vector<Node> Qp(nV);
796  vector<double> blend(nV);
797  vector<Node> out(nquad);
798 
799  // Calculate segment length for 2D spherigon routine.
800  double segLength = sqrt((v[0] - v[1]).abs2());
801 
802  // Perform Spherigon method to smooth manifold.
803  for (int j = 0; j < nquad; ++j)
804  {
805  Node P(0, x[j], y[j], z[j]);
806  Node N(0, 0, 0, 0);
807 
808  // Calculate generalised barycentric coordinates r[] and the
809  // Phong normal N = vN . r for this point of the element.
810  if (m_mesh->m_spaceDim == 2)
811  {
812  // In 2D the coordinates are given by a ratio of the
813  // segment length to the distance from one of the
814  // endpoints.
815  r[0] = sqrt((P - v[0]).abs2()) / segLength;
816  r[0] = max(min(1.0, r[0]), 0.0);
817  r[1] = 1.0 - r[0];
818 
819  // Calculate Phong normal.
820  N = vN[0] * r[0] + vN[1] * r[1];
821  }
822  else if (m_mesh->m_spaceDim == 3)
823  {
824  for (int k = 0; k < nV; ++k)
825  {
826  tmp[k] = P - v[k];
827  }
828 
829  // Calculate generalized barycentric coordinate system
830  // (see equation 6 of paper).
831  double weight = 0.0;
832  for (int k = 0; k < nV; ++k)
833  {
834  r[k] = 1.0;
835  for (int l = 0; l < nV - 2; ++l)
836  {
837  r[k] *= CrossProdMag(tmp[(k + l + 1) % nV],
838  tmp[(k + l + 2) % nV]);
839  }
840  weight += r[k];
841  }
842 
843  // Calculate Phong normal (equation 1).
844  for (int k = 0; k < nV; ++k)
845  {
846  r[k] /= weight;
847  N += vN[k] * r[k];
848  }
849  }
850 
851  // Normalise Phong normal.
852  N /= sqrt(N.abs2());
853 
854  for (int k = 0; k < nV; ++k)
855  {
856  // Perform steps denoted in equations 2, 3, 8 for C1
857  // smoothing.
858  double tmp1;
859  K[k] = P + N * ((v[k] - P).dot(N));
860  tmp1 = (v[k] - K[k]).dot(vN[k]) / (1.0 + N.dot(vN[k]));
861  Q[k] = K[k] + N * tmp1;
862  Qp[k] = v[k] - N * ((v[k] - P).dot(N));
863  }
864 
865  // Apply C1 blending function to the surface. TODO: Add
866  // option to do (more efficient) C0 blending function.
867  SuperBlend(r, Qp, P, blend);
868  P.m_x = P.m_y = P.m_z = 0.0;
869 
870  // Apply blending (equation 4).
871  for (int k = 0; k < nV; ++k)
872  {
873  P += Q[k] * blend[k];
874  }
875 
876  if ((boost::math::isnan)(P.m_x) || (boost::math::isnan)(P.m_y) ||
877  (boost::math::isnan)(P.m_z))
878  {
879  ASSERTL0(false,
880  "spherigon point is a nan. Check to see if "
881  "ply file is correct if using input normal file");
882  }
883  else
884  {
885  out[j] = P;
886  }
887  }
888 
889  // Push nodes into lines - TODO: face interior nodes.
890  // offset = 0 (seg), 1 (tri) or 2 (quad)
891  int offset = (int)e->GetConf().m_e - 2;
892 
893  for (int edge = 0; edge < e->GetEdgeCount(); ++edge)
894  {
895  eIt = visitedEdges.find(e->GetEdge(edge)->m_id);
896  if (eIt == visitedEdges.end())
897  {
898  bool reverseEdge =
899  !(v[vertMap[offset][edge][0]] == *(e->GetEdge(edge)->m_n1));
900 
901  // Clear existing curvature.
902  e->GetEdge(edge)->m_edgeNodes.clear();
903 
904  if (e->GetConf().m_e != LibUtilities::eTriangle)
905  {
906  for (int j = 1; j < nq - 1; ++j)
907  {
908  int v = edgeMap[offset][edge][0] +
909  j * edgeMap[offset][edge][1];
910  e->GetEdge(edge)->m_edgeNodes.push_back(
911  NodeSharedPtr(new Node(out[v])));
912  }
913  }
914  else
915  {
916  for (int j = 0; j < nq - 2; ++j)
917  {
918  int v = 3 + edge * (nq - 2) + j;
919  e->GetEdge(edge)->m_edgeNodes.push_back(
920  NodeSharedPtr(new Node(out[v])));
921  }
922  }
923 
924  if (reverseEdge)
925  {
926  reverse(e->GetEdge(edge)->m_edgeNodes.begin(),
927  e->GetEdge(edge)->m_edgeNodes.end());
928  }
929 
930  e->GetEdge(edge)->m_curveType =
932 
933  visitedEdges.insert(e->GetEdge(edge)->m_id);
934  }
935  }
936 
937  // Add face nodes in manifold and full 3D case, but not for 2D.
938  if (m_mesh->m_spaceDim == 3)
939  {
940  vector<NodeSharedPtr> volNodes;
941 
942  if (e->GetConf().m_e == LibUtilities::eQuadrilateral)
943  {
944  volNodes.resize((nq - 2) * (nq - 2));
945  for (int j = 1; j < nq - 1; ++j)
946  {
947  for (int k = 1; k < nq - 1; ++k)
948  {
949  int v = j * nq + k;
950  volNodes[(j - 1) * (nq - 2) + (k - 1)] =
951  NodeSharedPtr(new Node(out[v]));
952  }
953  }
954  }
955  else
956  {
957  for (int j = 3 + 3 * (nq - 2); j < nquad; ++j)
958  {
959  volNodes.push_back(NodeSharedPtr(new Node(out[j])));
960  }
961  }
962 
963  e->SetVolumeNodes(volNodes);
964  }
965  }
966 
967  // Copy face nodes back into 3D element faces
968  if (m_mesh->m_expDim == 3)
969  {
970  set<pair<int, int> >::iterator it;
971  int elmt = 0;
972  for (it = m_mesh->m_spherigonSurfs.begin();
973  it != m_mesh->m_spherigonSurfs.end();
974  ++it, ++elmt)
975  {
976  FaceSharedPtr f =
977  m_mesh->m_element[m_mesh->m_expDim][it->first]->GetFace(
978  it->second);
979 
980  f->m_faceNodes = el[elmt]->GetVolumeNodes();
981  f->m_curveType = f->m_vertexList.size() == 3
984  }
985  }
986 
987  if (normalsGenerated)
988  {
989  m_mesh->m_vertexNormals.clear();
990  }
991 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Basic information about an element.
Definition: Element.h:58
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
map< string, ConfigOption > m_config
List of configuration values.
Principle Modified Functions .
Definition: BasisType.h:49
MeshSharedPtr m_mesh
Mesh object.
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:57
Principle Orthogonal Functions .
Definition: BasisType.h:47
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:266
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:293
double CrossProdMag(Node &a, Node &b)
Calculate the magnitude of the cross product .
double NekDouble
boost::shared_ptr< InputPly > InputPlySharedPtr
Definition: InputPly.h:68
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:196
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:137
void SuperBlend(std::vector< double > &r, std::vector< Node > &Q, Node &P, std::vector< double > &blend)
Calculate the blending function for spherigon implementation.
void GenerateNormals(std::vector< ElementSharedPtr > &el, MeshSharedPtr &mesh)
Generate a set of approximate vertex normals to a surface represented by line segments in 2D and a hy...
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
Definition: Face.h:378
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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:191
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:68
boost::shared_ptr< NodalTriExp > NodalTriExpSharedPtr
Definition: NodalTriExp.h:379
void Nektar::Utilities::ProcessSpherigon::SuperBlend ( std::vector< double > &  r,
std::vector< Node > &  Q,
Node P,
std::vector< double > &  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 173 of file ProcessSpherigon.cpp.

References TOL_BLEND.

Referenced by Process().

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

123 {
124  double inv_mag;
125 
126  c.m_x = a.m_y * b.m_z - a.m_z * b.m_y;
127  c.m_y = a.m_z * b.m_x - a.m_x * b.m_z;
128  c.m_z = a.m_x * b.m_y - a.m_y * b.m_x;
129 
130  inv_mag = 1.0 / sqrt(c.m_x * c.m_x + c.m_y * c.m_y + c.m_z * c.m_z);
131 
132  c.m_x = c.m_x * inv_mag;
133  c.m_y = c.m_y * inv_mag;
134  c.m_z = c.m_z * inv_mag;
135 }
NekDouble m_y
Y-coordinate.
Definition: Node.h:314
NekDouble m_x
X-coordinate.
Definition: Node.h:312
NekDouble m_z
Z-coordinate.
Definition: Node.h:316

Member Data Documentation

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

Definition at line 56 of file ProcessSpherigon.h.