Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessSpherigon.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessSpherigon.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Apply Spherigon surface smoothing technique to a 3D mesh.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
39 
40 #include <LocalRegions/SegExp.h>
41 #include <LocalRegions/QuadExp.h>
42 #include <LocalRegions/TriExp.h>
44 
46 
47 #include "ProcessSpherigon.h"
48 
49 #include <boost/geometry.hpp>
50 #include <boost/geometry/geometries/point.hpp>
51 #include <boost/geometry/geometries/box.hpp>
52 #include <boost/geometry/index/rtree.hpp>
53 
54 namespace bg = boost::geometry;
55 namespace bgi = boost::geometry::index;
56 
57 using namespace std;
58 using namespace Nektar::NekMeshUtils;
59 
60 #define TOL_BLEND 1.0e-8
61 
62 namespace Nektar
63 {
64 namespace Utilities
65 {
66 ModuleKey ProcessSpherigon::className =
68  ModuleKey(eProcessModule, "spherigon"), ProcessSpherigon::create);
69 
70 /**
71  * @class ProcessSpherigon
72  *
73  * This class implements the spherigon surface smoothing technique which
74  * is documented in
75  *
76  * "The SPHERIGON: A Simple Polygon Patch for Smoothing Quickly your
77  * Polygonal Meshes": P. Volino and N. Magnenat Thalmann, Computer
78  * Animation Proceedings (1998).
79  *
80  * This implementation works in both a 2D manifold setting (for
81  * triangles and quadrilaterals embedded in 3-space) and in a full 3D
82  * enviroment (prisms, tetrahedra and hexahedra).
83  *
84  * No additional information needs to be supplied directly to the module
85  * in order for it to work. However, 3D elements do rely on the mapping
86  * Mesh::spherigonSurfs to be populated by the relevant input modules.
87  *
88  * Additionally, since the algorithm assumes normals are supplied which
89  * are perpendicular to the true surface defined at each vertex. If
90  * these are specified in Mesh::vertexNormals by the input module,
91  * better smoothing results can be obtained. Otherwise, normals are
92  * estimated by taking the average of all edge/face normals which
93  * connect to the vertex.
94  */
95 
96 /**
97  * @brief Default constructor.
98  */
99 ProcessSpherigon::ProcessSpherigon(MeshSharedPtr m) : 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 }
118 
119 /**
120  * @brief Destructor.
121  */
123 {
124 }
125 
126 /*
127  * @brief Calcuate the normalised cross product \f$ \vec{c} =
128  * \vec{a}\times\vec{b} / \| \vec{a}\times\vec{b} \| \f$.
129  */
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 }
144 
145 /**
146  * @brief Calculate the magnitude of the cross product \f$
147  * \vec{a}\times\vec{b} \f$.
148  */
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 }
156 
157 /**
158  * @brief Calculate the \f$ C^0 \f$ blending function for spherigon
159  * implementation.
160  *
161  * See equation (5) of the paper.
162  *
163  * @param r Barycentric coordinate.
164  */
165 double ProcessSpherigon::Blend(double r)
166 {
167  return r * r;
168 }
169 
170 /**
171  * @brief Calculate the \f$ C^1 \f$ blending function for spherigon
172  * implementation.
173  *
174  * See equation (10) of the paper.
175  *
176  * @param r Generalised barycentric coordinates of the point P.
177  * @param Q Vector of vertices denoting this triangle/quad.
178  * @param P Point in the triangle to apply blending to.
179  * @param blend The resulting blending components for each vertex.
180  */
181 void ProcessSpherigon::SuperBlend(vector<double> &r,
182  vector<Node> &Q,
183  Node &P,
184  vector<double> &blend)
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 }
223 
225  map<int,NodeSharedPtr> &surfverts)
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 }
287 
288 /**
289  * @brief Generate a set of approximate vertex normals to a surface
290  * represented by line segments in 2D and a hybrid
291  * triangular/quadrilateral mesh in 3D.
292  *
293  * This routine approximates the true vertex normals to a surface by
294  * averaging the normals of all edges/faces which connect to the
295  * vertex. It is better to use the exact surface normals which can be
296  * set in Mesh::vertexNormals, but where they are not supplied this
297  * routine calculates an approximation for the spherigon implementation.
298  *
299  * @param el Vector of elements denoting the surface mesh.
300  */
301 void ProcessSpherigon::GenerateNormals(std::vector<ElementSharedPtr> &el,
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 }
372 
373 /**
374  * @brief Perform the spherigon smoothing technique on the mesh.
375  */
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 
725  nq,
729  nq,
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 
759  nq,
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 }
1041 }
1042 }
void UnitCrossProd(NekMeshUtils::Node &a, NekMeshUtils::Node &b, NekMeshUtils::Node &c)
NEKMESHUTILS_EXPORT NekDouble dot(const Node &pSrc) const
Definition: Node.h:159
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
int PrintProgressbar(const int position, const int goal, const string message, int lastprogress=-1)
Prints a progressbar.
Definition: Progressbar.hpp:69
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
STL namespace.
Converter for Ply files.
Definition: InputPly.h:47
pair< ModuleType, string > ModuleKey
NekDouble m_y
Y-coordinate.
Definition: Node.h:401
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
virtual void Process()
Write mesh to output file.
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
Principle Orthogonal Functions .
Definition: BasisType.h:47
NEKMESHUTILS_EXPORT NekDouble abs2() const
Definition: Node.h:154
#define TOL_BLEND
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:270
A 0-dimensional vertex.
Definition: Point.h:49
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
Principle Orthogonal Functions .
Definition: BasisType.h:46
Defines a specification for a set of points.
Definition: Points.h:58
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:135
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< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
Abstract base class for processing modules.
NekDouble Blend(NekDouble r)
Calculate the blending function for spherigon implementation.
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:148
NekDouble m_z
Z-coordinate.
Definition: Node.h:403
std::pair< ModuleType, std::string > ModuleKey
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
Describes the specification for a Basis.
Definition: Basis.h:50
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:70
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
boost::shared_ptr< NodalTriExp > NodalTriExpSharedPtr
Definition: NodalTriExp.h:381