Nektar++
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 
36 #include <string>
37 using namespace std;
38 
39 #include "../MeshElements.h"
40 #include "ProcessSpherigon.h"
41 
42 #include <LocalRegions/SegExp.h>
43 #include <LocalRegions/QuadExp.h>
44 #include <LocalRegions/TriExp.h>
49 
50 #define TOL_BLEND 1.0e-8
51 
52 namespace Nektar
53 {
54  namespace Utilities
55  {
56  ModuleKey ProcessSpherigon::className =
58  ModuleKey(eProcessModule, "spherigon"),
59  ProcessSpherigon::create);
60 
61  /**
62  * @class ProcessSpherigon
63  *
64  * This class implements the spherigon surface smoothing technique which
65  * is documented in
66  *
67  * "The SPHERIGON: A Simple Polygon Patch for Smoothing Quickly your
68  * Polygonal Meshes": P. Volino and N. Magnenat Thalmann, Computer
69  * Animation Proceedings (1998).
70  *
71  * This implementation works in both a 2D manifold setting (for
72  * triangles and quadrilaterals embedded in 3-space) and in a full 3D
73  * enviroment (prisms, tetrahedra and hexahedra).
74  *
75  * No additional information needs to be supplied directly to the module
76  * in order for it to work. However, 3D elements do rely on the mapping
77  * Mesh::spherigonSurfs to be populated by the relevant input modules.
78  *
79  * Additionally, since the algorithm assumes normals are supplied which
80  * are perpendicular to the true surface defined at each vertex. If
81  * these are specified in Mesh::vertexNormals by the input module,
82  * better smoothing results can be obtained. Otherwise, normals are
83  * estimated by taking the average of all edge/face normals which
84  * connect to the vertex.
85  */
86 
87  /**
88  * @brief Default constructor.
89  */
90  ProcessSpherigon::ProcessSpherigon(MeshSharedPtr m) : ProcessModule(m)
91  {
92  m_config["N"] = ConfigOption(false, "5",
93  "Number of points to add to face edges.");
94  m_config["surf"] = ConfigOption(false, "-1",
95  "Tag identifying surface to process.");
96  m_config["BothTriFacesOnPrism"] = ConfigOption(true, "-1",
97  "Curve both triangular faces of prism on boundary.");
98  m_config["usenormalfile"] = ConfigOption(false,"NoFile",
99  "Use alternative file for Spherigon definition");
100  m_config["scalefile"] = ConfigOption(false,"1.0",
101  "Apply scaling factor to coordinates in file ");
102  m_config["normalnoise"] = ConfigOption(false,"NotSpecified",
103  "Add randowm noise to normals of amplitude AMP in specified region. input string is Amp,xmin,xmax,ymin,ymax,zmin,zmax");
104  }
105 
106  /**
107  * @brief Destructor.
108  */
110  {
111 
112  }
113 
114  /*
115  * @brief Calcuate the normalised cross product \f$ \vec{c} =
116  * \vec{a}\times\vec{b} / \| \vec{a}\times\vec{b} \| \f$.
117  */
119  {
120  double inv_mag;
121 
122  c.m_x = a.m_y*b.m_z - a.m_z*b.m_y;
123  c.m_y = a.m_z*b.m_x - a.m_x*b.m_z;
124  c.m_z = a.m_x*b.m_y - a.m_y*b.m_x;
125 
126  inv_mag = 1.0/sqrt(c.m_x*c.m_x + c.m_y*c.m_y + c.m_z*c.m_z);
127 
128  c.m_x = c.m_x*inv_mag;
129  c.m_y = c.m_y*inv_mag;
130  c.m_z = c.m_z*inv_mag;
131  }
132 
133  /**
134  * @brief Calculate the magnitude of the cross product \f$
135  * \vec{a}\times\vec{b} \f$.
136  */
138  {
139  double tmp1 = a.m_y*b.m_z - a.m_z*b.m_y;
140  double tmp2 = a.m_z*b.m_x - a.m_x*b.m_z;
141  double tmp3 = a.m_x*b.m_y - a.m_y*b.m_x;
142  return sqrt(tmp1*tmp1 + tmp2*tmp2 + tmp3*tmp3);
143  }
144 
145  /**
146  * @brief Calculate the \f$ C^0 \f$ blending function for spherigon
147  * implementation.
148  *
149  * See equation (5) of the paper.
150  *
151  * @param r Barycentric coordinate.
152  */
153  double ProcessSpherigon::Blend(double r)
154  {
155  return r*r;
156  }
157 
158  /**
159  * @brief Calculate the \f$ C^1 \f$ blending function for spherigon
160  * implementation.
161  *
162  * See equation (10) of the paper.
163  *
164  * @param r Generalised barycentric coordinates of the point P.
165  * @param Q Vector of vertices denoting this triangle/quad.
166  * @param P Point in the triangle to apply blending to.
167  * @param blend The resulting blending components for each vertex.
168  */
170  vector<double> &r, vector<Node> &Q, Node &P, vector<double> &blend)
171  {
172  vector<double> tmp(r.size());
173  double totBlend = 0.0;
174  int i;
175  int nV = r.size();
176 
177  for (i = 0; i < nV; ++i)
178  {
179  blend[i] = 0.0;
180  tmp [i] = (Q[i]-P).abs2();
181  }
182 
183  for (i = 0; i < nV; ++i)
184  {
185  int ip = (i+1) % nV, im = (i-1+nV) % nV;
186 
187  if (r[i] > TOL_BLEND && r[i] < (1-TOL_BLEND))
188  {
189  blend[i] = r[i]*r[i]*(
190  r[im]*r[im]*tmp[im]/(tmp[im] + tmp[i]) +
191  r[ip]*r[ip]*tmp[ip]/(tmp[ip] + tmp[i]));
192  totBlend += blend[i];
193  }
194  }
195 
196  for (i = 0; i < nV; ++i)
197  {
198  blend[i] /= totBlend;
199  if (r[i] >= (1-TOL_BLEND))
200  {
201  blend[i] = 1.0;
202  }
203  if (r[i] <= TOL_BLEND)
204  {
205  blend[i] = 0.0;
206  }
207  }
208  }
209 
210  /**
211  * @brief Generate a set of approximate vertex normals to a surface
212  * represented by line segments in 2D and a hybrid
213  * triangular/quadrilateral mesh in 3D.
214  *
215  * This routine approximates the true vertex normals to a surface by
216  * averaging the normals of all edges/faces which connect to the
217  * vertex. It is better to use the exact surface normals which can be
218  * set in Mesh::vertexNormals, but where they are not supplied this
219  * routine calculates an approximation for the spherigon implementation.
220  *
221  * @param el Vector of elements denoting the surface mesh.
222  */
223  void ProcessSpherigon::GenerateNormals(std::vector<ElementSharedPtr> &el,
224  MeshSharedPtr &mesh)
225  {
227 
228 
229 
230  for (int i = 0; i < el.size(); ++i)
231  {
232  ElementSharedPtr e = el[i];
233 
234  // Ensure that element is a line, triangle or quad.
235  ASSERTL0(e->GetConf().m_e == LibUtilities::eSegment ||
236  e->GetConf().m_e == LibUtilities::eTriangle ||
237  e->GetConf().m_e == LibUtilities::eQuadrilateral,
238  "Spherigon expansions must be lines, triangles or "
239  "quadrilaterals.");
240 
241  // Calculate normal for this element.
242  int nV = e->GetVertexCount();
243  vector<NodeSharedPtr> node(nV);
244 
245  for (int j = 0; j < nV; ++j)
246  {
247  node[j] = e->GetVertex(j);
248  }
249 
250  Node n;
251 
252  if (mesh->m_spaceDim == 3)
253  {
254  // Create two tangent vectors and take unit cross product.
255  Node v1 = *(node[1]) - *(node[0]);
256  Node v2 = *(node[2]) - *(node[0]);
257  UnitCrossProd(v1, v2, n);
258  }
259  else
260  {
261  // Calculate gradient vector and invert.
262  Node dx = *(node[1]) - *(node[0]);
263  dx /= sqrt(dx.abs2());
264  n.m_x = -dx.m_y;
265  n.m_y = dx.m_x;
266  n.m_z = 0;
267  }
268 
269  // Insert face normal into vertex normal list or add to existing
270  // value.
271  for (int j = 0; j < nV; ++j)
272  {
273  nIt = mesh->m_vertexNormals.find(e->GetVertex(j)->m_id);
274  if (nIt == mesh->m_vertexNormals.end())
275  {
276  mesh->m_vertexNormals[e->GetVertex(j)->m_id] = n;
277  }
278  else
279  {
280  nIt->second += n;
281  }
282  }
283  }
284 
285  // Normalize resulting vectors.
286  for (nIt = mesh->m_vertexNormals.begin();
287  nIt != mesh->m_vertexNormals.end (); ++nIt)
288  {
289  Node &n = mesh->m_vertexNormals[nIt->first];
290  n /= sqrt(n.abs2());
291  }
292  }
293 
294  /**
295  * @brief Perform the spherigon smoothing technique on the mesh.
296  */
298  {
299  ASSERTL0(m_mesh->m_spaceDim == 3 || m_mesh->m_spaceDim == 2,
300  "Spherigon implementation only valid in 2D/3D.");
301 
303  boost::unordered_set<int> visitedEdges;
304 
305  // First construct vector of elements to process.
306  vector<ElementSharedPtr> el;
307 
308  if (m_mesh->m_verbose)
309  {
310  cout << "ProcessSpherigon: Smoothing mesh..." << endl;
311  }
312 
313  if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3)
314  {
315  // Manifold case - copy expansion dimension.
316  el = m_mesh->m_element[m_mesh->m_expDim];
317  }
318  else if (m_mesh->m_spaceDim == m_mesh->m_expDim)
319  {
320  // Full 2D or 3D case - iterate over stored edges/faces and
321  // create segments/triangles representing those edges/faces.
322  set<pair<int,int> >::iterator it;
323  vector<int> t;
324  t.push_back(0);
325 
326 
327  // Construct list of spherigon edges/faces from a tag.
328  string surfTag = m_config["surf"].as<string>();
329  bool prismTag = m_config["BothTriFacesOnPrism"].beenSet;
330 
331  if (surfTag != "")
332  {
333  vector<unsigned int> surfs;
334  ParseUtils::GenerateSeqVector(surfTag.c_str(), surfs);
335  sort(surfs.begin(), surfs.end());
336 
337  m_mesh->m_spherigonSurfs.clear();
338  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
339  {
340  ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
341  int nSurf = m_mesh->m_expDim == 3 ? el->GetFaceCount() :
342  el->GetEdgeCount();
343 
344  for (int j = 0; j < nSurf; ++j)
345  {
346  int bl = el->GetBoundaryLink(j);
347  if (bl == -1)
348  {
349  continue;
350  }
351 
352  ElementSharedPtr bEl = m_mesh->m_element[m_mesh->m_expDim-1][bl];
353  vector<int> tags = bEl->GetTagList();
354  vector<int> inter;
355 
356  sort(tags.begin(), tags.end());
357  set_intersection(surfs.begin(), surfs.end(),
358  tags .begin(), tags .end(),
359  back_inserter(inter));
360 
361  if (inter.size() == 1)
362  {
363  m_mesh->m_spherigonSurfs.insert(make_pair(i, j));
364 
365  // Curve other tri face on Prism. Note could be
366  // problem on pyramid when implemented.
367  if(nSurf == 5 && prismTag)
368  {
369  // add other end of prism on boundary for
370  // smoothing
371  int triFace = j == 1 ? 3 : 1;
372  m_mesh->m_spherigonSurfs.insert(
373  make_pair(i, triFace));
374  }
375  }
376  }
377  }
378  }
379 
380  if (m_mesh->m_spherigonSurfs.size() == 0)
381  {
382  cerr << "WARNING: Spherigon surfaces have not been defined "
383  << "-- ignoring smoothing." << endl;
384  return;
385  }
386 
387  if (m_mesh->m_expDim == 3)
388  {
389  for (it = m_mesh->m_spherigonSurfs.begin();
390  it != m_mesh->m_spherigonSurfs.end (); ++it)
391  {
392  FaceSharedPtr f = m_mesh->m_element[m_mesh->m_expDim][it->first]->
393  GetFace(it->second);
394  vector<NodeSharedPtr> nodes = f->m_vertexList;
395  LibUtilities::ShapeType eType = (LibUtilities::ShapeType)(nodes.size());
396  ElmtConfig conf(eType, 1, false, false);
397 
398  // Create 2D element.
400  CreateInstance(eType,conf,nodes,t);
401 
402  // Copy vertices/edges from face.
403  for (int i = 0; i < f->m_vertexList.size(); ++i)
404  {
405  elmt->SetVertex(i, f->m_vertexList[i]);
406  }
407  for (int i = 0; i < f->m_edgeList.size(); ++i)
408  {
409  elmt->SetEdge(i, f->m_edgeList[i]);
410  }
411 
412  el.push_back(elmt);
413  }
414  }
415  else
416  {
417  for (it = m_mesh->m_spherigonSurfs.begin();
418  it != m_mesh->m_spherigonSurfs.end (); ++it)
419  {
420  EdgeSharedPtr edge = m_mesh->m_element[m_mesh->m_expDim][it->first]->
421  GetEdge(it->second);
422  vector<NodeSharedPtr> nodes;
424  ElmtConfig conf(eType, 1, false, false);
425 
426  nodes.push_back(edge->m_n1);
427  nodes.push_back(edge->m_n2);
428 
429  // Create 2D element.
431  CreateInstance(eType,conf,nodes,t);
432 
433  // Copy vertices/edges from original element.
434  elmt->SetVertex(0, nodes[0]);
435  elmt->SetVertex(1, nodes[1]);
436  elmt->SetEdge(0, m_mesh->m_element[m_mesh->m_expDim][it->first]->
437  GetEdge(it->second));
438  el.push_back(elmt);
439  }
440  }
441  }
442  else
443  {
444  ASSERTL0(false, "Spherigon expansions must be 2/3 dimensional");
445  }
446 
447 
448  // See if vertex normals have been generated. If they have not,
449  // approximate them by summing normals of surrounding elements.
450  bool normalsGenerated = false;
451 
452  // Read Normal file if one exists
453  std::string normalfile = m_config["usenormalfile"].as<string>();
454  if(normalfile.compare("NoFile") != 0)
455  {
456  NekDouble scale = m_config["scalefile"].as<NekDouble>();
457 
458  if (m_mesh->m_verbose)
459  {
460  cout << "Inputing normal file: " << normalfile << " with scaling of " << scale << endl;
461  }
462 
463  ifstream inply;
464  InputPlySharedPtr plyfile;
465 
466  inply.open(normalfile.c_str());
467  ASSERTL0(inply,string("Could not open input ply file: ") +
468  normalfile);
469 
470  int j;
471  MeshSharedPtr m = boost::shared_ptr<Mesh>(new Mesh());
472  plyfile = boost::shared_ptr<InputPly>(new InputPly(m));
473  plyfile->ReadPly(inply,scale);
474  plyfile->ProcessVertices();
475 
476  MeshSharedPtr plymesh = plyfile->GetMesh();
477  GenerateNormals(plymesh->m_element[plymesh->m_expDim],plymesh);
478 
479  // finaly find nearest vertex and set normal to mesh surface file normal.
480  // probably should have a hex tree search ?
481  Array<OneD, NekDouble> len2(plymesh->m_vertexSet.size());
482  Node minx(0,0.0,0.0,0.0), tmp,tmpsav;
483  NodeSet::iterator it;
485  map<int,NodeSharedPtr> surfverts;
486 
487  // make a map of normal vertices to visit based on elements el
488  for (int i = 0; i < el.size(); ++i)
489  {
490  ElementSharedPtr e = el[i];
491  int nV = e->GetVertexCount();
492  for (int j = 0; j < nV; ++j)
493  {
494  int id = e->GetVertex(j)->m_id;
495  surfverts[id] = e->GetVertex(j);
496  }
497  }
498 
499  //loop over all element in ply mesh and determine
500  //xmin,xmax,ymin,ymax as search criterion
501 
502 
503  NekDouble mindiff,diff;
504  int cntmin;
505 
506  if (m_mesh->m_verbose)
507  {
508  cout << "\t Processing surface normals " << endl;
509  }
510  int cnt = 0;
511  map<int,int> locnorm;
512  for(vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt,++cnt)
513  {
514  mindiff = 1e12;
515 
516  for(j = 0, it = plymesh->m_vertexSet.begin(); it != plymesh->m_vertexSet.end(); ++it, ++j)
517  {
518  tmp = *(vIt->second)- *(*it);
519  diff = tmp.abs2();
520 
521  if(diff < mindiff)
522  {
523  mindiff = diff;
524  cntmin = (*it)->m_id;
525  tmpsav = tmp;
526  }
527  }
528  locnorm[cntmin] = vIt->first;
529 
530  ASSERTL1(cntmin < plymesh->m_vertexNormals.size(),"cntmin is out of range");
531  m_mesh->m_vertexNormals[vIt->first] = plymesh->m_vertexNormals[cntmin];
532 
533  }
534  if (m_mesh->m_verbose)
535  {
536  cout << "\t end of processing surface normals " << endl;
537  }
538  normalsGenerated = true;
539  }
540  else if (m_mesh->m_vertexNormals.size() == 0)
541  {
543  normalsGenerated = true;
544  }
545 
546 
547  // See if we should add noise to normals
548  std::string normalnoise = m_config["normalnoise"].as<string>();
549  if(normalnoise.compare("NotSpecified") != 0)
550  {
551  vector<NekDouble> values;
552  ASSERTL0(ParseUtils::GenerateUnOrderedVector(normalnoise.c_str(),values),"Failed to interpret normal noise string");
553 
554  int nvalues = values.size()/2;
555  NekDouble amp = values[0];
556 
557 
558  if (m_mesh->m_verbose)
559  {
560  cout << "\t adding noise to normals of amplitude "<< amp << " in range: ";
561  for(int i = 0; i < nvalues; ++i)
562  {
563  cout << values[2*i+1] <<"," << values[2*i+2] << " ";
564  }
565  cout << endl;
566  }
567 
569  map<int,NodeSharedPtr> surfverts;
570 
571  // make a map of normal vertices to visit based on elements el
572  for (int i = 0; i < el.size(); ++i)
573  {
574  ElementSharedPtr e = el[i];
575  int nV = e->GetVertexCount();
576  for (int j = 0; j < nV; ++j)
577  {
578  int id = e->GetVertex(j)->m_id;
579  surfverts[id] = e->GetVertex(j);
580  }
581  }
582 
583  for(vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt)
584  {
585  bool AddNoise = false;
586 
587  for(int i = 0; i < nvalues; ++i)
588  {
589  // check to see if point is in range
590  switch(nvalues)
591  {
592  case 1:
593  {
594  if(((vIt->second)->m_x > values[2*i+1])&&((vIt->second)->m_x < values[2*i+2]))
595  {
596  AddNoise = true;
597  }
598  }
599  break;
600  case 2:
601  {
602  if(((vIt->second)->m_x > values[2*i+1])&&
603  ((vIt->second)->m_x < values[2*i+2])&&
604  ((vIt->second)->m_y > values[2*i+3])&&
605  ((vIt->second)->m_y < values[2*i+4]))
606  {
607  AddNoise = true;
608  }
609  }
610  break;
611  case 3:
612  {
613  if(((vIt->second)->m_x > values[2*i+1])&&
614  ((vIt->second)->m_x < values[2*i+2])&&
615  ((vIt->second)->m_y > values[2*i+3])&&
616  ((vIt->second)->m_y < values[2*i+4])&&
617  ((vIt->second)->m_z > values[2*i+5])&&
618  ((vIt->second)->m_z < values[2*i+6]))
619 
620  {
621  AddNoise = true;
622  }
623  }
624  break;
625  }
626 
627  if(AddNoise)
628  {
629  // generate random unit vector;
630  Node rvec(0,rand(),rand(),rand());
631  rvec *= values[0]/sqrt(rvec.abs2());
632 
633  Node normal = m_mesh->m_vertexNormals[vIt->first];
634 
635  normal += rvec;
636  normal /= sqrt(normal.abs2());
637 
638  m_mesh->m_vertexNormals[vIt->first] = normal;
639  }
640  }
641  }
642  }
643 
644 
645  // Allocate storage for interior points.
646  int nq = m_config["N"].as<int>();
647  int nquad = m_mesh->m_spaceDim == 3 ? nq*nq : nq;
648  Array<OneD, NekDouble> x(nq*nq);
649  Array<OneD, NekDouble> y(nq*nq);
650  Array<OneD, NekDouble> z(nq*nq);
651 
652  Array<OneD, NekDouble> xc(nq*nq);
653  Array<OneD, NekDouble> yc(nq*nq);
654  Array<OneD, NekDouble> zc(nq*nq);
655 
656  ASSERTL0(nq > 2, "Number of points must be greater than 2.");
657 
669 
670  Array<OneD, NekDouble> xnodal(nq*(nq+1)/2), ynodal(nq*(nq+1)/2);
671  stdtri->GetNodalPoints(xnodal, ynodal);
672 
673  int edgeMap[3][4][2] = {
674  {{0, 1}, {-1, -1}, {-1, -1 }, {-1, -1}}, // seg
675  {{0, 1}, {nq-1, nq}, {nq*(nq-1), -nq}, {-1, -1}}, // tri
676  {{0, 1}, {nq-1, nq}, {nq*nq-1, -1 }, {nq*(nq-1), -nq}} // quad
677  };
678 
679  int vertMap[3][4][2] = {
680  {{0, 1}, {0, 0}, {0, 0}, {0, 0}}, // seg
681  {{0, 1}, {1, 2}, {2, 3}, {0, 0}}, // tri
682  {{0, 1}, {1, 2}, {2, 3}, {3, 0}}, // quad
683  };
684 
685  for (int i = 0; i < el.size(); ++i)
686  {
687  // Construct a Nektar++ element to obtain coordinate points
688  // inside the element. TODO: Add options for various
689  // nodal/tensor point distributions + number of points to add.
690  ElementSharedPtr e = el[i];
691 
696 
697  if (e->GetConf().m_e == LibUtilities::eSegment)
698  {
700  boost::dynamic_pointer_cast<SpatialDomains::SegGeom>(
701  e->GetGeom(m_mesh->m_spaceDim));
704  B2, geom);
705  seg->GetCoords(x,y,z);
706  nquad = nq;
707  }
708  else if (e->GetConf().m_e == LibUtilities::eTriangle)
709  {
711  boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
712  e->GetGeom(3));
716  B0, B1, LibUtilities::eNodalTriElec, geom);
717 
718  Array<OneD, NekDouble> coord(2);
719  tri->GetCoords(xc,yc,zc);
720  nquad = nq*(nq+1)/2;
721 
722  for (int j = 0; j < nquad; ++j)
723  {
724  coord[0] = xnodal[j];
725  coord[1] = ynodal[j];
726  x[j] = stdtri->PhysEvaluate(coord, xc);
727  }
728 
729  for (int j = 0; j < nquad; ++j)
730  {
731  coord[0] = xnodal[j];
732  coord[1] = ynodal[j];
733  y[j] = stdtri->PhysEvaluate(coord, yc);
734  }
735 
736  for (int j = 0; j < nquad; ++j)
737  {
738  coord[0] = xnodal[j];
739  coord[1] = ynodal[j];
740  z[j] = stdtri->PhysEvaluate(coord, zc);
741  }
742  }
743  else if (e->GetConf().m_e == LibUtilities::eQuadrilateral)
744  {
746  boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
747  e->GetGeom(3));
750  B2, B2, geom);
751  quad->GetCoords(x,y,z);
752  nquad = nq*nq;
753  }
754  else
755  {
756  ASSERTL0(false, "Unknown expansion type.");
757  }
758 
759  // Zero z-coordinate in 2D.
760  if (m_mesh->m_spaceDim == 2)
761  {
762  Vmath::Zero(nquad, z, 1);
763  }
764 
765  // Find vertex normals.
766  int nV = e->GetVertexCount();
767  vector<Node> v, vN;
768  for (int j = 0; j < nV; ++j)
769  {
770  v.push_back(*(e->GetVertex(j)));
771  ASSERTL1(m_mesh->m_vertexNormals.count(v[j].m_id) != 0,"Normal has not been defined");
772  vN.push_back(m_mesh->m_vertexNormals[v[j].m_id]);
773  }
774 
775  vector<Node> tmp (nV);
776  vector<double> r (nV);
777  vector<Node> K (nV);
778  vector<Node> Q (nV);
779  vector<Node> Qp (nV);
780  vector<double> blend(nV);
781  vector<Node> out (nquad);
782 
783  // Calculate segment length for 2D spherigon routine.
784  double segLength = sqrt((v[0] - v[1]).abs2());
785 
786  // Perform Spherigon method to smooth manifold.
787  for (int j = 0; j < nquad; ++j)
788  {
789  Node P(0, x[j], y[j], z[j]);
790  Node N(0,0,0,0);
791 
792  // Calculate generalised barycentric coordinates r[] and the
793  // Phong normal N = vN . r for this point of the element.
794  if (m_mesh->m_spaceDim == 2)
795  {
796  // In 2D the coordinates are given by a ratio of the
797  // segment length to the distance from one of the
798  // endpoints.
799  r[0] = sqrt((P - v[0]).abs2()) / segLength;
800  r[0] = max(min(1.0, r[0]), 0.0);
801  r[1] = 1.0 - r[0];
802 
803  // Calculate Phong normal.
804  N = vN[0]*r[0] + vN[1]*r[1];
805  }
806  else if (m_mesh->m_spaceDim == 3)
807  {
808  for (int k = 0; k < nV; ++k)
809  {
810  tmp[k] = P - v[k];
811  }
812 
813  // Calculate generalized barycentric coordinate system
814  // (see equation 6 of paper).
815  double weight = 0.0;
816  for (int k = 0; k < nV; ++k)
817  {
818  r[k] = 1.0;
819  for (int l = 0; l < nV-2; ++l)
820  {
821  r[k] *= CrossProdMag(tmp[(k+l+1) % nV],
822  tmp[(k+l+2) % nV]);
823  }
824  weight += r[k];
825  }
826 
827  // Calculate Phong normal (equation 1).
828  for (int k = 0; k < nV; ++k)
829  {
830  r[k] /= weight;
831  N += vN[k]*r[k];
832  }
833  }
834 
835  // Normalise Phong normal.
836  N /= sqrt(N.abs2());
837 
838  for (int k = 0; k < nV; ++k)
839  {
840  // Perform steps denoted in equations 2, 3, 8 for C1
841  // smoothing.
842  double tmp1;
843  K [k] = P+N*((v[k]-P).dot(N));
844  tmp1 = (v[k]-K[k]).dot(vN[k]) / (1.0 + N.dot(vN[k]));
845  Q [k] = K[k] + N*tmp1;
846  Qp[k] = v[k] - N*((v[k]-P).dot(N));
847  }
848 
849  // Apply C1 blending function to the surface. TODO: Add
850  // option to do (more efficient) C0 blending function.
851  SuperBlend(r, Qp, P, blend);
852  P.m_x = P.m_y = P.m_z = 0.0;
853 
854  // Apply blending (equation 4).
855  for (int k = 0; k < nV; ++k)
856  {
857  P += Q[k]*blend[k];
858  }
859 
860  out[j] = P;
861  }
862 
863  // Push nodes into lines - TODO: face interior nodes.
864  // offset = 0 (seg), 1 (tri) or 2 (quad)
865  int offset = (int)e->GetConf().m_e-2;
866 
867  for (int edge = 0; edge < e->GetEdgeCount(); ++edge)
868  {
869  eIt = visitedEdges.find(e->GetEdge(edge)->m_id);
870  if (eIt == visitedEdges.end())
871  {
872  bool reverseEdge = !(v[vertMap[offset][edge][0]] ==
873  *(e->GetEdge(edge)->m_n1));
874 
875  // Clear existing curvature.
876  e->GetEdge(edge)->m_edgeNodes.clear();
877 
878  if (e->GetConf().m_e != LibUtilities::eTriangle)
879  {
880  for (int j = 1; j < nq-1; ++j)
881  {
882  int v = edgeMap[offset][edge][0] +
883  j*edgeMap[offset][edge][1];
884  e->GetEdge(edge)->m_edgeNodes.push_back(
885  NodeSharedPtr(new Node(out[v])));
886  }
887  }
888  else
889  {
890  for (int j = 0; j < nq-2; ++j)
891  {
892  int v = 3 + edge*(nq-2) + j;
893  e->GetEdge(edge)->m_edgeNodes.push_back(
894  NodeSharedPtr(new Node(out[v])));
895  }
896  }
897 
898  if (reverseEdge)
899  {
900  reverse(e->GetEdge(edge)->m_edgeNodes.begin(),
901  e->GetEdge(edge)->m_edgeNodes.end());
902  }
903 
904  e->GetEdge(edge)->m_curveType =
906 
907  visitedEdges.insert(e->GetEdge(edge)->m_id);
908  }
909  }
910 
911  // Add face nodes in manifold and full 3D case, but not for 2D.
912  if (m_mesh->m_spaceDim == 3)
913  {
914  vector<NodeSharedPtr> volNodes;
915 
916  if (e->GetConf().m_e == LibUtilities::eQuadrilateral)
917  {
918  volNodes.resize((nq-2)*(nq-2));
919  for (int j = 1; j < nq-1; ++j)
920  {
921  for (int k = 1; k < nq-1; ++k)
922  {
923  int v = j*nq+k;
924  volNodes[(j-1)*(nq-2)+(k-1)] =
925  NodeSharedPtr(new Node(out[v]));
926  }
927  }
928  }
929  else
930  {
931  for (int j = 3+3*(nq-2); j < nquad; ++j)
932  {
933  volNodes.push_back(NodeSharedPtr(new Node(out[j])));
934  }
935  }
936 
937  e->SetVolumeNodes(volNodes);
938  }
939  }
940 
941  // Copy face nodes back into 3D element faces
942  if (m_mesh->m_expDim == 3)
943  {
944  set<pair<int,int> >::iterator it;
945  int elmt = 0;
946  for (it = m_mesh->m_spherigonSurfs.begin();
947  it != m_mesh->m_spherigonSurfs.end (); ++it, ++elmt)
948  {
949  FaceSharedPtr f = m_mesh->m_element[m_mesh->m_expDim][it->first]->
950  GetFace(it->second);
951 
952  f->m_faceNodes = el[elmt]->GetVolumeNodes();
953  f->m_curveType = f->m_vertexList.size() == 3 ?
956  }
957  }
958 
959  if (normalsGenerated)
960  {
961  m_mesh->m_vertexNormals.clear();
962  }
963  }
964  }
965 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
pair< ModuleType, string > ModuleKey
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void GenerateNormals(vector< ElementSharedPtr > &el, MeshSharedPtr &mesh)
Generate a set of approximate vertex normals to a surface represented by line segments in 2D and a hy...
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: MeshElements.h:318
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:46
MeshSharedPtr m_mesh
Mesh object.
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
Definition: MeshElements.h:550
NekDouble abs2() const
Definition: MeshElements.h:155
boost::shared_ptr< Node > NodeSharedPtr
Shared pointer to a Node.
Definition: MeshElements.h:195
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:78
boost::shared_ptr< StdNodalTriExp > StdNodalTriExpSharedPtr
virtual void Process()
Write mesh to output file.
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
NekDouble m_x
X-coordinate.
Definition: MeshElements.h:185
boost::shared_ptr< Element > ElementSharedPtr
Shared pointer to an element.
Definition: MeshElements.h:63
Principle Orthogonal Functions .
Definition: BasisType.h:47
void SuperBlend(vector< double > &r, vector< Node > &Q, Node &P, vector< double > &blend)
Calculate the blending function for spherigon implementation.
#define TOL_BLEND
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:255
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
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:284
NekDouble m_y
Y-coordinate.
Definition: MeshElements.h:187
double CrossProdMag(Node &a, Node &b)
Calculate the magnitude of the cross product .
double NekDouble
Basic information about an element.
Definition: MeshElements.h:583
Represents a point in the domain.
Definition: MeshElements.h:74
void UnitCrossProd(Node &a, Node &b, Node &c)
boost::shared_ptr< InputPly > InputPlySharedPtr
Definition: InputPly.h:65
Represents a command-line configuration option.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
NekDouble dot(const Node &pSrc) const
Definition: MeshElements.h:160
NekDouble m_z
Z-coordinate.
Definition: MeshElements.h:189
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
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:127
#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
ElementFactory & GetElementFactory()
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