Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 
38 
39 #include <LocalRegions/SegExp.h>
40 #include <LocalRegions/QuadExp.h>
41 #include <LocalRegions/TriExp.h>
43 
45 
46 #include "ProcessSpherigon.h"
47 
48 using namespace std;
49 using namespace Nektar::NekMeshUtils;
50 
51 #define TOL_BLEND 1.0e-8
52 
53 namespace Nektar
54 {
55 namespace Utilities
56 {
57 
58 ModuleKey ProcessSpherigon::className =
60  ModuleKey(eProcessModule, "spherigon"), ProcessSpherigon::create);
61 
62 /**
63  * @class ProcessSpherigon
64  *
65  * This class implements the spherigon surface smoothing technique which
66  * is documented in
67  *
68  * "The SPHERIGON: A Simple Polygon Patch for Smoothing Quickly your
69  * Polygonal Meshes": P. Volino and N. Magnenat Thalmann, Computer
70  * Animation Proceedings (1998).
71  *
72  * This implementation works in both a 2D manifold setting (for
73  * triangles and quadrilaterals embedded in 3-space) and in a full 3D
74  * enviroment (prisms, tetrahedra and hexahedra).
75  *
76  * No additional information needs to be supplied directly to the module
77  * in order for it to work. However, 3D elements do rely on the mapping
78  * Mesh::spherigonSurfs to be populated by the relevant input modules.
79  *
80  * Additionally, since the algorithm assumes normals are supplied which
81  * are perpendicular to the true surface defined at each vertex. If
82  * these are specified in Mesh::vertexNormals by the input module,
83  * better smoothing results can be obtained. Otherwise, normals are
84  * estimated by taking the average of all edge/face normals which
85  * connect to the vertex.
86  */
87 
88 /**
89  * @brief Default constructor.
90  */
91 ProcessSpherigon::ProcessSpherigon(MeshSharedPtr m) : 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 }
110 
111 /**
112  * @brief Destructor.
113  */
115 {
116 }
117 
118 /*
119  * @brief Calcuate the normalised cross product \f$ \vec{c} =
120  * \vec{a}\times\vec{b} / \| \vec{a}\times\vec{b} \| \f$.
121  */
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 }
136 
137 /**
138  * @brief Calculate the magnitude of the cross product \f$
139  * \vec{a}\times\vec{b} \f$.
140  */
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 }
148 
149 /**
150  * @brief Calculate the \f$ C^0 \f$ blending function for spherigon
151  * implementation.
152  *
153  * See equation (5) of the paper.
154  *
155  * @param r Barycentric coordinate.
156  */
157 double ProcessSpherigon::Blend(double r)
158 {
159  return r * r;
160 }
161 
162 /**
163  * @brief Calculate the \f$ C^1 \f$ blending function for spherigon
164  * implementation.
165  *
166  * See equation (10) of the paper.
167  *
168  * @param r Generalised barycentric coordinates of the point P.
169  * @param Q Vector of vertices denoting this triangle/quad.
170  * @param P Point in the triangle to apply blending to.
171  * @param blend The resulting blending components for each vertex.
172  */
173 void ProcessSpherigon::SuperBlend(vector<double> &r,
174  vector<Node> &Q,
175  Node &P,
176  vector<double> &blend)
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 }
215 
216 /**
217  * @brief Generate a set of approximate vertex normals to a surface
218  * represented by line segments in 2D and a hybrid
219  * triangular/quadrilateral mesh in 3D.
220  *
221  * This routine approximates the true vertex normals to a surface by
222  * averaging the normals of all edges/faces which connect to the
223  * vertex. It is better to use the exact surface normals which can be
224  * set in Mesh::vertexNormals, but where they are not supplied this
225  * routine calculates an approximation for the spherigon implementation.
226  *
227  * @param el Vector of elements denoting the surface mesh.
228  */
229 void ProcessSpherigon::GenerateNormals(std::vector<ElementSharedPtr> &el,
230  MeshSharedPtr &mesh)
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 }
298 
299 /**
300  * @brief Perform the spherigon smoothing technique on the mesh.
301  */
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 
676  nq,
680  nq,
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 
710  nq,
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 }
992 }
993 }
NEKMESHUTILS_EXPORT NekDouble dot(const Node &pSrc) const
Definition: Node.h:153
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Basic information about an element.
Definition: Element.h:58
pair< ModuleType, string > ModuleKey
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
STL namespace.
Converter for Ply files.
Definition: InputPly.h:47
MeshSharedPtr m_mesh
Mesh object.
NekDouble m_y
Y-coordinate.
Definition: Node.h:314
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:57
Principle Orthogonal Functions .
Definition: BasisType.h:47
NEKMESHUTILS_EXPORT NekDouble abs2() const
Definition: Node.h:148
#define TOL_BLEND
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
Defines a specification for a set of points.
Definition: Points.h:58
boost::shared_ptr< QuadExp > QuadExpSharedPtr
Definition: QuadExp.h:293
double CrossProdMag(Node &a, Node &b)
Calculate the magnitude of the cross product .
double NekDouble
void UnitCrossProd(Node &a, Node &b, Node &c)
boost::shared_ptr< InputPly > InputPlySharedPtr
Definition: InputPly.h:68
Represents a command-line configuration option.
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
NekDouble m_x
X-coordinate.
Definition: Node.h:312
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
NekDouble m_z
Z-coordinate.
Definition: Node.h:316
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
Describes the specification for a Basis.
Definition: Basis.h:50
double Blend(double r)
Calculate the blending function for spherigon implementation.
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:68
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
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:379