Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Member Functions | Private Attributes | Friends | List of all members
Nektar::NekMeshUtils::FaceMesh Class Reference

class for surface meshes on individual surfaces (paramter plane meshes) More...

#include <FaceMesh.h>

Collaboration diagram for Nektar::NekMeshUtils::FaceMesh:
Collaboration graph
[legend]

Public Member Functions

 FaceMesh (const int id, MeshSharedPtr m, CADSurfSharedPtr cad, OctreeSharedPtr oct, const std::map< int, CurveMeshSharedPtr > &cmeshes)
 Default constructor. More...
 
 FaceMesh (const int id, MeshSharedPtr m, CADSurfSharedPtr cad, OctreeSharedPtr oct, const std::map< int, CurveMeshSharedPtr > &cmeshes, const NekDouble b)
 constructor for building with boundary layer More...
 
void Mesh ()
 mesh exectuation command More...
 

Private Member Functions

void Stretching ()
 Calculate the paramter plane streching factor. More...
 
void Smoothing ()
 performs node smoothing on face More...
 
void DiagonalSwap ()
 performs diagonal swapping of edges More...
 
void BuildLocalMesh ()
 build a local version of mesh elements More...
 
void OptimiseLocalMesh ()
 function which calls the optimisation routines More...
 
bool Validate ()
 Validate the surface mesh base on the octree and real dimensions of the edges. More...
 
void OrientateCurves ()
 Get the boundries of the surface and extracts the nodes from the curve meshes in the correct order. More...
 
void AddNewPoint (Array< OneD, NekDouble > uv)
 adds a new stiener point to the triangulation for meshing More...
 
void MakeBL ()
 adds a quad layer around any interior loops More...
 

Private Attributes

MeshSharedPtr m_mesh
 mesh pointer More...
 
CADSurfSharedPtr m_cadsurf
 CAD surface. More...
 
OctreeSharedPtr m_octree
 Octree object. More...
 
std::map< int, CurveMeshSharedPtrm_curvemeshes
 Map of the curve meshes which bound the surfaces. More...
 
std::vector< EdgeLoopm_edgeloops
 data structure containing the edges, their order and oreientation for the surface More...
 
int m_id
 id of the surface mesh More...
 
std::vector< std::vector
< NodeSharedPtr > > 
orderedLoops
 list of boundary nodes in their order loops More...
 
std::vector< NodeSharedPtrm_stienerpoints
 list of stiener points in the triangulation More...
 
NekDouble m_str
 pplane stretching More...
 
std::vector< std::vector
< NodeSharedPtr > > 
m_connec
 triangle connectiviities More...
 
NodeSet m_localNodes
 local set of nodes More...
 
EdgeSet m_localEdges
 local set of edges More...
 
std::vector< ElementSharedPtrm_localElements
 local list of elements More...
 
NekDouble m_bl
 boundary layer thickness More...
 
bool m_makebl
 should build boundary layer More...
 
std::vector< std::pair
< NodeSharedPtr, NodeSharedPtr > > 
blpairs
 list of node links between node on loop and its corresponding interior node in quads More...
 

Friends

class MemoryManager< FaceMesh >
 

Detailed Description

class for surface meshes on individual surfaces (paramter plane meshes)

Definition at line 53 of file FaceMesh.h.

Constructor & Destructor Documentation

Nektar::NekMeshUtils::FaceMesh::FaceMesh ( const int  id,
MeshSharedPtr  m,
CADSurfSharedPtr  cad,
OctreeSharedPtr  oct,
const std::map< int, CurveMeshSharedPtr > &  cmeshes 
)
inline

Default constructor.

Definition at line 61 of file FaceMesh.h.

References m_cadsurf, m_edgeloops, and m_makebl.

66  : m_mesh(m), m_cadsurf(cad), m_octree(oct), m_curvemeshes(cmeshes),
67  m_id(id)
68 
69  {
70  m_edgeloops = m_cadsurf->GetEdges();
71  m_makebl = false;
72  };
std::map< int, CurveMeshSharedPtr > m_curvemeshes
Map of the curve meshes which bound the surfaces.
Definition: FaceMesh.h:151
OctreeSharedPtr m_octree
Octree object.
Definition: FaceMesh.h:149
int m_id
id of the surface mesh
Definition: FaceMesh.h:156
std::vector< EdgeLoop > m_edgeloops
data structure containing the edges, their order and oreientation for the surface ...
Definition: FaceMesh.h:154
bool m_makebl
should build boundary layer
Definition: FaceMesh.h:174
CADSurfSharedPtr m_cadsurf
CAD surface.
Definition: FaceMesh.h:147
MeshSharedPtr m_mesh
mesh pointer
Definition: FaceMesh.h:145
Nektar::NekMeshUtils::FaceMesh::FaceMesh ( const int  id,
MeshSharedPtr  m,
CADSurfSharedPtr  cad,
OctreeSharedPtr  oct,
const std::map< int, CurveMeshSharedPtr > &  cmeshes,
const NekDouble  b 
)
inline

constructor for building with boundary layer

Definition at line 77 of file FaceMesh.h.

References m_cadsurf, m_edgeloops, and m_makebl.

83  : m_mesh(m), m_cadsurf(cad), m_octree(oct), m_curvemeshes(cmeshes),
84  m_id(id), m_bl(b)
85 
86  {
87  m_edgeloops = m_cadsurf->GetEdges();
88  m_makebl = true;
89  };
std::map< int, CurveMeshSharedPtr > m_curvemeshes
Map of the curve meshes which bound the surfaces.
Definition: FaceMesh.h:151
OctreeSharedPtr m_octree
Octree object.
Definition: FaceMesh.h:149
NekDouble m_bl
boundary layer thickness
Definition: FaceMesh.h:172
int m_id
id of the surface mesh
Definition: FaceMesh.h:156
std::vector< EdgeLoop > m_edgeloops
data structure containing the edges, their order and oreientation for the surface ...
Definition: FaceMesh.h:154
bool m_makebl
should build boundary layer
Definition: FaceMesh.h:174
CADSurfSharedPtr m_cadsurf
CAD surface.
Definition: FaceMesh.h:147
MeshSharedPtr m_mesh
mesh pointer
Definition: FaceMesh.h:145

Member Function Documentation

void Nektar::NekMeshUtils::FaceMesh::AddNewPoint ( Array< OneD, NekDouble uv)
private

adds a new stiener point to the triangulation for meshing

Definition at line 976 of file FaceMesh.cpp.

977 {
978  // adds a new point but checks that there are no other points nearby first
979  Array<OneD, NekDouble> np = m_cadsurf->P(uv);
980  NekDouble npDelta = m_octree->Query(np);
981 
982  NodeSharedPtr n = boost::shared_ptr<Node>(
983  new Node(m_mesh->m_numNodes++, np[0], np[1], np[2]));
984 
985  bool add = true;
986 
987  for (int i = 0; i < orderedLoops.size(); i++)
988  {
989  for (int j = 0; j < orderedLoops[i].size(); j++)
990  {
991  NekDouble r = orderedLoops[i][j]->Distance(n);
992 
993  if (r < npDelta / 2.0)
994  {
995  add = false;
996  break;
997  }
998  }
999  }
1000 
1001  if (add)
1002  {
1003  for (int i = 0; i < m_stienerpoints.size(); i++)
1004  {
1005  NekDouble r = m_stienerpoints[i]->Distance(n);
1006 
1007  if (r < npDelta / 2.0)
1008  {
1009  add = false;
1010  break;
1011  }
1012  }
1013  }
1014 
1015  if (add)
1016  {
1017  n->SetCADSurf(m_id, m_cadsurf, uv);
1018  m_stienerpoints.push_back(n);
1019  }
1020 }
std::vector< std::vector< NodeSharedPtr > > orderedLoops
list of boundary nodes in their order loops
Definition: FaceMesh.h:158
OctreeSharedPtr m_octree
Octree object.
Definition: FaceMesh.h:149
int m_id
id of the surface mesh
Definition: FaceMesh.h:156
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
double NekDouble
std::vector< NodeSharedPtr > m_stienerpoints
list of stiener points in the triangulation
Definition: FaceMesh.h:160
CADSurfSharedPtr m_cadsurf
CAD surface.
Definition: FaceMesh.h:147
MeshSharedPtr m_mesh
mesh pointer
Definition: FaceMesh.h:145
void Nektar::NekMeshUtils::FaceMesh::BuildLocalMesh ( )
private

build a local version of mesh elements

Definition at line 745 of file FaceMesh.cpp.

References Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::NekMeshUtils::GetElementFactory(), and Nektar::iterator.

746 {
747  /*************************
748  // build a local set of nodes edges and elemenets for optimstaion prior to
749  putting them into m_mesh
750  */
751  int tricomp = m_mesh->m_numcomp++;
752 
753  // first build quads is bl surface
754  if (m_makebl)
755  {
756  int quadcomp = m_mesh->m_numcomp++;
757  for (int i = 0; i < blpairs.size() - 1;
758  i++) // this wont work for more than one sym plane loop
759  {
760  ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
761  vector<NodeSharedPtr> ns;
762  ns.push_back(blpairs[i].first);
763  ns.push_back(blpairs[i].second);
764  ns.push_back(blpairs[i + 1].second);
765  ns.push_back(blpairs[i + 1].first);
766  vector<int> tags;
767  tags.push_back(quadcomp);
769  LibUtilities::eQuadrilateral, conf, ns, tags);
770  E->CADSurfId = m_id;
771  m_localElements.push_back(E);
772  }
773  {
774  ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
775  vector<NodeSharedPtr> ns;
776  ns.push_back(blpairs.back().first);
777  ns.push_back(blpairs.back().second);
778  ns.push_back(blpairs[0].second);
779  ns.push_back(blpairs[0].first);
780  vector<int> tags;
781  tags.push_back(quadcomp);
783  LibUtilities::eQuadrilateral, conf, ns, tags);
784  E->CADSurfId = m_id;
785  m_localElements.push_back(E);
786  }
787  for (int i = 0; i < m_localElements.size(); i++)
788  {
789  vector<NodeSharedPtr> nods = m_localElements[i]->GetVertexList();
790  for (int j = 0; j < nods.size(); j++)
791  {
792  // nodes are already unique some will insert some wont
793  m_localNodes.insert(nods[j]);
794  }
795  vector<EdgeSharedPtr> edgs = m_localElements[i]->GetEdgeList();
796  for (int j = 0; j < edgs.size(); j++)
797  {
798  // look for edge in m_mesh edgeset from curves
799  EdgeSet::iterator s = m_mesh->m_edgeSet.find(edgs[j]);
800  if (!(s == m_mesh->m_edgeSet.end()))
801  {
802  edgs[j] = *s;
803  m_localElements[i]->SetEdge(j, edgs[j]);
804  }
805 
806  pair<EdgeSet::iterator, bool> test =
807  m_localEdges.insert(edgs[j]);
808 
809  if (test.second)
810  {
811  (*test.first)
812  ->m_elLink.push_back(
813  pair<ElementSharedPtr, int>(m_localElements[i], j));
814  }
815  else
816  {
817  m_localElements[i]->SetEdge(j, *test.first);
818  (*test.first)
819  ->m_elLink.push_back(
820  pair<ElementSharedPtr, int>(m_localElements[i], j));
821  }
822  }
823  m_localElements[i]->SetId(i);
824  }
825  }
826 
827  for (int i = 0; i < m_connec.size(); i++)
828  {
829  ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);
830 
831  vector<int> tags;
832  tags.push_back(tricomp);
834  LibUtilities::eTriangle, conf, m_connec[i], tags);
835  E->CADSurfId = m_id;
836 
837  vector<NodeSharedPtr> nods = E->GetVertexList();
838  for (int j = 0; j < nods.size(); j++)
839  {
840  // nodes are already unique some will insert some wont
841  m_localNodes.insert(nods[j]);
842  }
843  vector<EdgeSharedPtr> edgs = E->GetEdgeList();
844  for (int j = 0; j < edgs.size(); j++)
845  {
846  // look for edge in m_mesh edgeset from curves
847  EdgeSet::iterator s = m_mesh->m_edgeSet.find(edgs[j]);
848  if (!(s == m_mesh->m_edgeSet.end()))
849  {
850  edgs[j] = *s;
851  E->SetEdge(j, edgs[j]);
852  }
853 
854  pair<EdgeSet::iterator, bool> test = m_localEdges.insert(edgs[j]);
855 
856  if (test.second)
857  {
858  (*test.first)
859  ->m_elLink.push_back(pair<ElementSharedPtr, int>(E, j));
860  }
861  else
862  {
863  E->SetEdge(j, *test.first);
864  (*test.first)
865  ->m_elLink.push_back(pair<ElementSharedPtr, int>(E, j));
866  }
867  }
868  E->SetId(m_localElements.size());
869  m_localElements.push_back(E);
870  }
871 }
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
std::vector< ElementSharedPtr > m_localElements
local list of elements
Definition: FaceMesh.h:170
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
std::vector< std::vector< NodeSharedPtr > > m_connec
triangle connectiviities
Definition: FaceMesh.h:164
std::vector< std::pair< NodeSharedPtr, NodeSharedPtr > > blpairs
list of node links between node on loop and its corresponding interior node in quads ...
Definition: FaceMesh.h:177
EdgeSet m_localEdges
local set of edges
Definition: FaceMesh.h:168
int m_id
id of the surface mesh
Definition: FaceMesh.h:156
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
bool m_makebl
should build boundary layer
Definition: FaceMesh.h:174
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
MeshSharedPtr m_mesh
mesh pointer
Definition: FaceMesh.h:145
NodeSet m_localNodes
local set of nodes
Definition: FaceMesh.h:166
void Nektar::NekMeshUtils::FaceMesh::DiagonalSwap ( )
private

performs diagonal swapping of edges

TODO fix this bit of code which figures out the node defect, doesnt work on quads or relfex angles

Definition at line 363 of file FaceMesh.cpp.

References ASSERTL0, Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::NekMeshUtils::GetElementFactory(), and Nektar::iterator.

364 {
365  /// TODO fix this bit of code which figures out the node defect, doesnt work
366  /// on quads or relfex angles
367  map<int, int> idealConnec;
368  map<int, int> actualConnec;
369  map<int, vector<EdgeSharedPtr> > nodetoedge;
370  // figure out ideal node count and actual node count
371  EdgeSet::iterator eit;
372  for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
373  {
374  nodetoedge[(*eit)->m_n1->m_id].push_back(*eit);
375  nodetoedge[(*eit)->m_n2->m_id].push_back(*eit);
376  }
377  NodeSet::iterator nit;
378  for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
379  {
380  if ((*nit)->GetNumCadCurve() == 0)
381  {
382  // node is interior
383  idealConnec[(*nit)->m_id] = 6;
384  }
385  else
386  {
387  // need to identify the two other nodes on the boundary to find
388  // interior angle
389  vector<NodeSharedPtr> ns;
390  vector<EdgeSharedPtr> e = nodetoedge[(*nit)->m_id];
391  for (int i = 0; i < e.size(); i++)
392  {
393  if (!e[i]->onCurve)
394  continue; // the linking nodes are not going to exist on
395  // interior edges
396 
397  if (e[i]->m_n1 == (*nit))
398  ns.push_back(e[i]->m_n2);
399  else
400  ns.push_back(e[i]->m_n1);
401  }
402  ASSERTL0(ns.size() == 2,
403  "failed to find 2 nodes in the angle system");
404 
405  idealConnec[(*nit)->m_id] =
406  ceil((*nit)->Angle(ns[0], ns[1]) / 3.142 * 3) + 1;
407  }
408  }
409  for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
410  {
411  actualConnec[(*nit)->m_id] = nodetoedge[(*nit)->m_id].size();
412  }
413 
414  // edgeswapping fun times
415  // perfrom edge swap based on node defect and then angle
416  for (int q = 0; q < 4; q++)
417  {
418  int edgesStart = m_localEdges.size();
419  EdgeSet edges = m_localEdges;
420  m_localEdges.clear();
421 
422  int swappedEdges = 0;
423 
425 
426  for (it = edges.begin(); it != edges.end(); it++)
427  {
428  EdgeSharedPtr e = *it;
429 
430  if (e->m_elLink.size() != 2)
431  {
432  m_localEdges.insert(e);
433  continue;
434  }
435  if (e->m_elLink[0].first->GetConf().m_e ==
437  e->m_elLink[1].first->GetConf().m_e ==
439  {
440  m_localEdges.insert(e);
441  continue;
442  }
443 
444  ElementSharedPtr tri1 = e->m_elLink[0].first;
445  ElementSharedPtr tri2 = e->m_elLink[1].first;
446 
447  NodeSharedPtr n1 = e->m_n1;
448  NodeSharedPtr n2 = e->m_n2;
449 
450  vector<NodeSharedPtr> nt = tri1->GetVertexList();
451 
452  // identify node a,b,c,d of the swapping
453  NodeSharedPtr A, B, C, D;
454  if (nt[0] != n1 && nt[0] != n2)
455  {
456  C = nt[0];
457  B = nt[1];
458  A = nt[2];
459  }
460  else if (nt[1] != n1 && nt[1] != n2)
461  {
462  C = nt[1];
463  B = nt[2];
464  A = nt[0];
465  }
466  else if (nt[2] != n1 && nt[2] != n2)
467  {
468  C = nt[2];
469  B = nt[0];
470  A = nt[1];
471  }
472  else
473  {
474  ASSERTL0(false, "failed to identify verticies in tri1");
475  }
476 
477  nt = tri2->GetVertexList();
478 
479  if (nt[0] != n1 && nt[0] != n2)
480  {
481  D = nt[0];
482  }
483  else if (nt[1] != n1 && nt[1] != n2)
484  {
485  D = nt[1];
486  }
487  else if (nt[2] != n1 && nt[2] != n2)
488  {
489  D = nt[2];
490  }
491  else
492  {
493  ASSERTL0(false, "failed to identify verticies in tri2");
494  }
495 
496  // determine signed area of alternate config
497  Array<OneD, NekDouble> ai, bi, ci, di;
498  ai = A->GetCADSurfInfo(m_id);
499  bi = B->GetCADSurfInfo(m_id);
500  ci = C->GetCADSurfInfo(m_id);
501  di = D->GetCADSurfInfo(m_id);
502 
503  NekDouble CDA, CBD;
504 
505  CDA = 0.5 * (-di[0] * ci[1] + ai[0] * ci[1] + ci[0] * di[1] -
506  ai[0] * di[1] - ci[0] * ai[1] + di[0] * ai[1]);
507 
508  CBD = 0.5 * (-bi[0] * ci[1] + di[0] * ci[1] + ci[0] * bi[1] -
509  di[0] * bi[1] - ci[0] * di[1] + bi[0] * di[1]);
510 
511  // if signed area of the swapping triangles is less than zero
512  // that configuration is invalid and swap cannot be performed
513  if (!(CDA > 0.001 && CBD > 0.001))
514  {
515  m_localEdges.insert(e);
516  continue;
517  }
518 
519  bool swap = false; // assume do not swap
520 
521  if (q < 2)
522  {
523  int nodedefectbefore = 0;
524  nodedefectbefore +=
525  abs(actualConnec[A->m_id] - idealConnec[A->m_id]);
526  nodedefectbefore +=
527  abs(actualConnec[B->m_id] - idealConnec[B->m_id]);
528  nodedefectbefore +=
529  abs(actualConnec[C->m_id] - idealConnec[C->m_id]);
530  nodedefectbefore +=
531  abs(actualConnec[D->m_id] - idealConnec[D->m_id]);
532 
533  int nodedefectafter = 0;
534  nodedefectafter +=
535  abs(actualConnec[A->m_id] - 1 - idealConnec[A->m_id]);
536  nodedefectafter +=
537  abs(actualConnec[B->m_id] - 1 - idealConnec[B->m_id]);
538  nodedefectafter +=
539  abs(actualConnec[C->m_id] + 1 - idealConnec[C->m_id]);
540  nodedefectafter +=
541  abs(actualConnec[D->m_id] + 1 - idealConnec[D->m_id]);
542 
543  if (nodedefectafter < nodedefectbefore)
544  {
545  swap = true;
546  }
547  }
548  else
549  {
550  NekDouble minanglebefore = C->Angle(A, B);
551  minanglebefore = min(minanglebefore, A->Angle(B, C));
552  minanglebefore = min(minanglebefore, B->Angle(A, C));
553  minanglebefore = min(minanglebefore, B->Angle(A, D));
554  minanglebefore = min(minanglebefore, A->Angle(B, D));
555  minanglebefore = min(minanglebefore, D->Angle(A, B));
556 
557  NekDouble minangleafter = C->Angle(B, D);
558  minangleafter = min(minangleafter, D->Angle(B, C));
559  minangleafter = min(minangleafter, B->Angle(C, D));
560  minangleafter = min(minangleafter, C->Angle(A, D));
561  minangleafter = min(minangleafter, A->Angle(C, D));
562  minangleafter = min(minangleafter, D->Angle(A, C));
563 
564  if (minangleafter > minanglebefore)
565  {
566  swap = true;
567  }
568  }
569 
570  if (swap)
571  {
572  actualConnec[A->m_id]--;
573  actualConnec[B->m_id]--;
574  actualConnec[C->m_id]++;
575  actualConnec[D->m_id]++;
576 
577  // make the 4 other edges
578  EdgeSharedPtr CA, AD, DB, BC, CAt, ADt, DBt, BCt;
579  CAt = boost::shared_ptr<Edge>(new Edge(C, A));
580  ADt = boost::shared_ptr<Edge>(new Edge(A, D));
581  DBt = boost::shared_ptr<Edge>(new Edge(D, B));
582  BCt = boost::shared_ptr<Edge>(new Edge(B, C));
583 
584  vector<EdgeSharedPtr> es = tri1->GetEdgeList();
585  for (int i = 0; i < 3; i++)
586  {
587  if (es[i] == CAt)
588  {
589  CA = es[i];
590  }
591  if (es[i] == BCt)
592  {
593  BC = es[i];
594  }
595  }
596  es = tri2->GetEdgeList();
597  for (int i = 0; i < 3; i++)
598  {
599  if (es[i] == DBt)
600  {
601  DB = es[i];
602  }
603  if (es[i] == ADt)
604  {
605  AD = es[i];
606  }
607  }
608 
609  // now sort out links for the 4 edges surrounding the patch
610  vector<pair<ElementSharedPtr, int> > links;
611 
612  links = CA->m_elLink;
613  CA->m_elLink.clear();
614  for (int i = 0; i < links.size(); i++)
615  {
616  if (links[i].first->GetId() == tri1->GetId())
617  continue;
618  CA->m_elLink.push_back(links[i]);
619  }
620 
621  links = BC->m_elLink;
622  BC->m_elLink.clear();
623  for (int i = 0; i < links.size(); i++)
624  {
625  if (links[i].first->GetId() == tri1->GetId())
626  continue;
627  BC->m_elLink.push_back(links[i]);
628  }
629 
630  links = AD->m_elLink;
631  AD->m_elLink.clear();
632  for (int i = 0; i < links.size(); i++)
633  {
634  if (links[i].first->GetId() == tri2->GetId())
635  continue;
636  AD->m_elLink.push_back(links[i]);
637  }
638 
639  links = DB->m_elLink;
640  DB->m_elLink.clear();
641  for (int i = 0; i < links.size(); i++)
642  {
643  if (links[i].first->GetId() == tri2->GetId())
644  continue;
645  DB->m_elLink.push_back(links[i]);
646  }
647 
648  EdgeSharedPtr newe = boost::shared_ptr<Edge>(new Edge(C, D));
649 
650  vector<NodeSharedPtr> t1, t2;
651  t1.push_back(B);
652  t1.push_back(D);
653  t1.push_back(C);
654  t2.push_back(A);
655  t2.push_back(C);
656  t2.push_back(D);
657 
658  ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);
659  vector<int> tags = tri1->GetTagList();
660 
661  int id1 = tri1->GetId();
662  int id2 = tri2->GetId();
663 
665  LibUtilities::eTriangle, conf, t1, tags);
666  tags = tri2->GetTagList();
668  LibUtilities::eTriangle, conf, t2, tags);
669 
670  ntri1->SetId(id1);
671  ntri2->SetId(id2);
672  ntri1->CADSurfId = m_id;
673  ntri2->CADSurfId = m_id;
674 
675  vector<EdgeSharedPtr> t1es = ntri1->GetEdgeList();
676  for (int i = 0; i < 3; i++)
677  {
678  if (t1es[i] == DB)
679  {
680  ntri1->SetEdge(i, DB);
681  DB->m_elLink.push_back(
682  pair<ElementSharedPtr, int>(ntri1, i));
683  }
684  else if (t1es[i] == BC)
685  {
686  ntri1->SetEdge(i, BC);
687  BC->m_elLink.push_back(
688  pair<ElementSharedPtr, int>(ntri1, i));
689  }
690  else if (t1es[i] == newe)
691  {
692  ntri1->SetEdge(i, newe);
693  newe->m_elLink.push_back(
694  pair<ElementSharedPtr, int>(ntri1, i));
695  }
696  else
697  {
698  ASSERTL0(false, "weird edge in new tri 1");
699  }
700  }
701  vector<EdgeSharedPtr> t2es = ntri2->GetEdgeList();
702  for (int i = 0; i < 3; i++)
703  {
704  if (t2es[i] == CA)
705  {
706  ntri2->SetEdge(i, CA);
707  CA->m_elLink.push_back(
708  pair<ElementSharedPtr, int>(ntri2, i));
709  }
710  else if (t2es[i] == AD)
711  {
712  ntri2->SetEdge(i, AD);
713  AD->m_elLink.push_back(
714  pair<ElementSharedPtr, int>(ntri2, i));
715  }
716  else if (t2es[i] == newe)
717  {
718  ntri2->SetEdge(i, newe);
719  newe->m_elLink.push_back(
720  pair<ElementSharedPtr, int>(ntri2, i));
721  }
722  else
723  {
724  ASSERTL0(false, "weird edge in new tri 2");
725  }
726  }
727 
728  m_localEdges.insert(newe);
729 
730  m_localElements[id1] = ntri1;
731  m_localElements[id2] = ntri2;
732 
733  swappedEdges++;
734  }
735  else
736  {
737  m_localEdges.insert(e);
738  }
739  }
740 
741  ASSERTL0(m_localEdges.size() == edgesStart, "mismatch edge count");
742  }
743 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
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
std::vector< ElementSharedPtr > m_localElements
local list of elements
Definition: FaceMesh.h:170
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
EdgeSet m_localEdges
local set of edges
Definition: FaceMesh.h:168
int m_id
id of the surface mesh
Definition: FaceMesh.h:156
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
double NekDouble
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:196
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
boost::unordered_set< EdgeSharedPtr, EdgeHash > EdgeSet
Definition: Edge.h:222
NodeSet m_localNodes
local set of nodes
Definition: FaceMesh.h:166
void Nektar::NekMeshUtils::FaceMesh::MakeBL ( )
private

adds a quad layer around any interior loops

Definition at line 147 of file FaceMesh.cpp.

References ASSERTL0.

148 {
149  for (int i = 1; i < orderedLoops.size(); i++) // dont do this to first loop
150  {
151  for (int j = 0; j < orderedLoops[i].size(); j++)
152  {
153  // for each of the nodes make a new node which exists off the
154  // surface
155  vector<pair<int, CADSurfSharedPtr> > surfs =
156  orderedLoops[i][j]->GetCADSurfs();
157  ASSERTL0(surfs.size() > 1,
158  "point does not have enough surfs to make quad");
159 
160  Array<OneD, NekDouble> AN(3, 0.0);
161  // make a averaged normal ignoring surface mid
162  for (int s = 0; s < surfs.size(); s++)
163  {
164  if (surfs[s].first == m_id)
165  continue; // does not contribute to norm
166 
167  Array<OneD, NekDouble> uv =
168  orderedLoops[i][j]->GetCADSurfInfo(surfs[s].first);
169  Array<OneD, NekDouble> N = surfs[s].second->N(uv);
170  for (int k = 0; k < 3; k++)
171  AN[k] += N[k];
172  }
173 
174  // renormalise normal
175  NekDouble mag = sqrt(AN[0] * AN[0] + AN[1] * AN[1] + AN[2] * AN[2]);
176  for (int k = 0; k < 3; k++)
177  AN[k] /= mag;
178 
179  Array<OneD, NekDouble> loc = orderedLoops[i][j]->GetLoc();
180  Array<OneD, NekDouble> tp(3);
181  for (int k = 0; k < 3; k++)
182  tp[k] = m_bl * AN[k] + loc[k];
183  // project tp onto to surface to get new point
184  Array<OneD, NekDouble> uv(2);
185  m_cadsurf->ProjectTo(tp, uv);
186  NodeSharedPtr nn = boost::shared_ptr<Node>(
187  new Node(m_mesh->m_numNodes++, tp[0], tp[1], tp[2]));
188  nn->SetCADSurf(m_id, m_cadsurf, uv);
189 
190  blpairs.push_back(
191  pair<NodeSharedPtr, NodeSharedPtr>(orderedLoops[i][j], nn));
192  // place the new node into ordered loops for the surface
193  // triangulation to work With
194  orderedLoops[i][j] = nn;
195  }
196  }
197 }
std::vector< std::vector< NodeSharedPtr > > orderedLoops
list of boundary nodes in their order loops
Definition: FaceMesh.h:158
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
NekDouble m_bl
boundary layer thickness
Definition: FaceMesh.h:172
std::vector< std::pair< NodeSharedPtr, NodeSharedPtr > > blpairs
list of node links between node on loop and its corresponding interior node in quads ...
Definition: FaceMesh.h:177
int m_id
id of the surface mesh
Definition: FaceMesh.h:156
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
double NekDouble
CADSurfSharedPtr m_cadsurf
CAD surface.
Definition: FaceMesh.h:147
MeshSharedPtr m_mesh
mesh pointer
Definition: FaceMesh.h:145
void Nektar::NekMeshUtils::FaceMesh::Mesh ( )

mesh exectuation command

Definition at line 49 of file FaceMesh.cpp.

References ASSERTL0, and Nektar::iterator.

50 {
51  Stretching();
52 
54 
55  if (m_makebl)
56  {
57  MakeBL();
58  }
59 
60  int numPoints = 0;
61  for (int i = 0; i < orderedLoops.size(); i++)
62  {
63  numPoints += orderedLoops[i].size();
64  }
65 
66  stringstream ss;
67  ss << "3 points required for triangulation, " << numPoints << " in loop"
68  << endl;
69  ss << "curves: ";
70  for (int i = 0; i < m_edgeloops.size(); i++)
71  {
72  for (int j = 0; j < m_edgeloops[i].edges.size(); j++)
73  {
74  ss << m_edgeloops[i].edges[j]->GetId() << " ";
75  }
76  }
77 
78  ASSERTL0(numPoints > 2, ss.str());
79 
80  // create interface to triangle thirdparty library
81  TriangleInterfaceSharedPtr pplanemesh =
83 
84  vector<Array<OneD, NekDouble> > centers;
85  for (int i = 0; i < m_edgeloops.size(); i++)
86  {
87  centers.push_back(m_edgeloops[i].center);
88  }
89 
90  pplanemesh->Assign(orderedLoops, centers, m_id, m_str);
91 
92  pplanemesh->Mesh();
93 
94  pplanemesh->Extract(m_connec);
95 
96  bool repeat = true;
97  int meshcounter = 1;
98 
99  while (repeat)
100  {
101  repeat = Validate();
102  if (!repeat)
103  {
104  break;
105  }
106  m_connec.clear();
107  pplanemesh->AssignStiener(m_stienerpoints);
108  pplanemesh->Mesh();
109  pplanemesh->Extract(m_connec);
110  meshcounter++;
111  }
112 
113  BuildLocalMesh();
114 
116 
117  // clear local element links
118  EdgeSet::iterator eit;
119  for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
120  {
121  (*eit)->m_elLink.clear();
122  }
123 
124  // make new elements and add to list from list of nodes and connectivity
125  // from triangle
126  for (int i = 0; i < m_localElements.size(); i++)
127  {
128  vector<EdgeSharedPtr> e = m_localElements[i]->GetEdgeList();
129  for (int j = 0; j < e.size(); j++)
130  {
131  e[j]->m_elLink.clear();
132  }
133  m_mesh->m_element[m_mesh->m_expDim].push_back(m_localElements[i]);
134  }
135 
136  cout << "\r "
137  " ";
138  cout << scientific << "\r\t\tFace " << m_id << endl
139  << "\t\t\tNodes: " << m_localNodes.size() << endl
140  << "\t\t\tEdges: " << m_localEdges.size() << endl
141  << "\t\t\tTriangles: " << m_localElements.size() << endl
142  << "\t\t\tLoops: " << m_edgeloops.size() << endl
143  << "\t\t\tSTR: " << m_str << endl
144  << endl;
145 }
std::vector< std::vector< NodeSharedPtr > > orderedLoops
list of boundary nodes in their order loops
Definition: FaceMesh.h:158
bool Validate()
Validate the surface mesh base on the octree and real dimensions of the edges.
Definition: FaceMesh.cpp:921
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
boost::shared_ptr< TriangleInterface > TriangleInterfaceSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Stretching()
Calculate the paramter plane streching factor.
Definition: FaceMesh.cpp:873
std::vector< ElementSharedPtr > m_localElements
local list of elements
Definition: FaceMesh.h:170
std::vector< std::vector< NodeSharedPtr > > m_connec
triangle connectiviities
Definition: FaceMesh.h:164
EdgeSet m_localEdges
local set of edges
Definition: FaceMesh.h:168
int m_id
id of the surface mesh
Definition: FaceMesh.h:156
void OrientateCurves()
Get the boundries of the surface and extracts the nodes from the curve meshes in the correct order...
Definition: FaceMesh.cpp:1022
std::vector< NodeSharedPtr > m_stienerpoints
list of stiener points in the triangulation
Definition: FaceMesh.h:160
std::vector< EdgeLoop > m_edgeloops
data structure containing the edges, their order and oreientation for the surface ...
Definition: FaceMesh.h:154
void MakeBL()
adds a quad layer around any interior loops
Definition: FaceMesh.cpp:147
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
bool m_makebl
should build boundary layer
Definition: FaceMesh.h:174
NekDouble m_str
pplane stretching
Definition: FaceMesh.h:162
void BuildLocalMesh()
build a local version of mesh elements
Definition: FaceMesh.cpp:745
void OptimiseLocalMesh()
function which calls the optimisation routines
Definition: FaceMesh.cpp:199
MeshSharedPtr m_mesh
mesh pointer
Definition: FaceMesh.h:145
NodeSet m_localNodes
local set of nodes
Definition: FaceMesh.h:166
void Nektar::NekMeshUtils::FaceMesh::OptimiseLocalMesh ( )
private

function which calls the optimisation routines

Definition at line 199 of file FaceMesh.cpp.

200 {
201  DiagonalSwap();
202 
203  Smoothing();
204 
205  DiagonalSwap();
206 
207  Smoothing();
208 }
void DiagonalSwap()
performs diagonal swapping of edges
Definition: FaceMesh.cpp:363
void Smoothing()
performs node smoothing on face
Definition: FaceMesh.cpp:210
void Nektar::NekMeshUtils::FaceMesh::OrientateCurves ( )
private

Get the boundries of the surface and extracts the nodes from the curve meshes in the correct order.

Definition at line 1022 of file FaceMesh.cpp.

References ASSERTL0.

1023 {
1024  // create list of bounding loop nodes
1025  for (int i = 0; i < m_edgeloops.size(); i++)
1026  {
1027  vector<NodeSharedPtr> cE;
1028  for (int j = 0; j < m_edgeloops[i].edges.size(); j++)
1029  {
1030  int cid = m_edgeloops[i].edges[j]->GetId();
1031  vector<NodeSharedPtr> edgePoints =
1032  m_curvemeshes[cid]->GetMeshPoints();
1033 
1034  int numPoints = m_curvemeshes[cid]->GetNumPoints();
1035 
1036  if (m_edgeloops[i].edgeo[j] == 0)
1037  {
1038  for (int k = 0; k < numPoints - 1; k++)
1039  {
1040  cE.push_back(edgePoints[k]);
1041  }
1042  }
1043  else
1044  {
1045  for (int k = numPoints - 1; k > 0; k--)
1046  {
1047  cE.push_back(edgePoints[k]);
1048  }
1049  }
1050  }
1051  orderedLoops.push_back(cE);
1052  }
1053 
1054  // loops made need to orientate on which is biggest and define holes
1055  for (int i = 0; i < orderedLoops.size(); i++)
1056  {
1057  NekDouble area = 0.0;
1058  for (int j = 0; j < orderedLoops[i].size() - 1; j++)
1059  {
1060  Array<OneD, NekDouble> n1info, n2info;
1061  n1info = orderedLoops[i][j]->GetCADSurfInfo(m_id);
1062  n2info = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);
1063 
1064  area += -n2info[1] * (n2info[0] - n1info[0]) +
1065  n1info[0] * (n2info[1] - n1info[1]);
1066  }
1067  area *= 0.5;
1068  m_edgeloops[i].area = area;
1069  }
1070 
1071  int ct = 0;
1072 
1073  do
1074  {
1075  ct = 0;
1076  for (int i = 0; i < m_edgeloops.size() - 1; i++)
1077  {
1078  if (fabs(m_edgeloops[i].area) < fabs(m_edgeloops[i + 1].area))
1079  {
1080  // swap
1081  vector<NodeSharedPtr> orderedlooptmp = orderedLoops[i];
1082  EdgeLoop edgeLoopstmp = m_edgeloops[i];
1083 
1084  orderedLoops[i] = orderedLoops[i + 1];
1085  m_edgeloops[i] = m_edgeloops[i + 1];
1086 
1087  orderedLoops[i + 1] = orderedlooptmp;
1088  m_edgeloops[i + 1] = edgeLoopstmp;
1089 
1090  ct += 1;
1091  }
1092  }
1093 
1094  } while (ct > 0);
1095 
1096  for (int i = 0; i < orderedLoops.size(); i++)
1097  {
1098  NodeSharedPtr n1, n2;
1099 
1100  n1 = orderedLoops[i][0];
1101  n2 = orderedLoops[i][1];
1102 
1103  Array<OneD, NekDouble> n1info, n2info;
1104  n1info = n1->GetCADSurfInfo(m_id);
1105  n2info = n2->GetCADSurfInfo(m_id);
1106 
1107  Array<OneD, NekDouble> N(2);
1108  NekDouble mag = sqrt((n1info[0] - n2info[0]) * (n1info[0] - n2info[0]) +
1109  (n1info[1] - n2info[1]) * (n1info[1] - n2info[1]));
1110  ASSERTL0(mag > 1e-30, "infinity");
1111  N[0] = -1.0 * (n2info[1] - n1info[1]) / mag;
1112  N[1] = (n2info[0] - n1info[0]) / mag;
1113 
1114  Array<OneD, NekDouble> P(2);
1115  P[0] = (n1info[0] + n2info[0]) / 2.0 + 1e-6 * N[0];
1116  P[1] = (n1info[1] + n2info[1]) / 2.0 + 1e-6 * N[1];
1117 
1118  // now test to see if p is inside or outside the shape
1119  // vector to the right
1120  int intercepts = 0;
1121  for (int j = 0; j < orderedLoops[i].size() - 1; j++)
1122  {
1123  Array<OneD, NekDouble> nt1, nt2;
1124  nt1 = orderedLoops[i][j]->GetCADSurfInfo(m_id);
1125  nt2 = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);
1126 
1127  if (fabs(nt2[1] - nt1[1]) < 1e-30)
1128  continue;
1129 
1130  NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
1131  NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
1132 
1133  if (!(lam < 0) && !(lam > 1) && S > 0)
1134  {
1135  intercepts++;
1136  }
1137  }
1138  {
1139  Array<OneD, NekDouble> nt1, nt2;
1140  nt1 = orderedLoops[i].back()->GetCADSurfInfo(m_id);
1141  nt2 = orderedLoops[i][0]->GetCADSurfInfo(m_id);
1142 
1143  if (fabs(nt2[1] - nt1[1]) < 1e-30)
1144  continue;
1145 
1146  NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
1147  NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
1148 
1149  if (!(lam < 0) && !(lam > 1) && S > 0)
1150  {
1151  intercepts++;
1152  }
1153  }
1154  if (intercepts % 2 == 0)
1155  {
1156  P[0] = (n1info[0] + n2info[0]) / 2.0 - 1e-6 * N[0];
1157  P[1] = (n1info[1] + n2info[1]) / 2.0 - 1e-6 * N[1];
1158  intercepts = 0;
1159  for (int j = 0; j < orderedLoops[i].size() - 1; j++)
1160  {
1161  Array<OneD, NekDouble> nt1, nt2;
1162  nt1 = orderedLoops[i][j]->GetCADSurfInfo(m_id);
1163  nt2 = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);
1164 
1165  if (fabs(nt2[1] - nt1[1]) < 1e-30)
1166  continue;
1167 
1168  NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
1169  NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
1170 
1171  if (!(lam < 0) && !(lam > 1) && S > 0)
1172  {
1173  intercepts++;
1174  }
1175  }
1176  {
1177  Array<OneD, NekDouble> nt1, nt2;
1178  nt1 = orderedLoops[i].back()->GetCADSurfInfo(m_id);
1179  nt2 = orderedLoops[i][0]->GetCADSurfInfo(m_id);
1180 
1181  if (fabs(nt2[1] - nt1[1]) < 1e-30)
1182  continue;
1183 
1184  NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
1185  NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
1186 
1187  if (!(lam < 0) && !(lam > 1) && S > 0)
1188  {
1189  intercepts++;
1190  }
1191  }
1192  if (intercepts % 2 == 0)
1193  {
1194  cerr << "still failed to find point inside loop" << endl;
1195  }
1196  }
1197 
1198  m_edgeloops[i].center = P;
1199  }
1200 
1201  if (m_edgeloops[0].area < 0) // reverse the first uvLoop
1202  {
1203  vector<NodeSharedPtr> tmp = orderedLoops[0];
1204  reverse(tmp.begin(), tmp.end());
1205  orderedLoops[0] = tmp;
1206  }
1207 
1208  for (int i = 1; i < orderedLoops.size(); i++)
1209  {
1210  if (m_edgeloops[i].area > 0) // reverse the loop
1211  {
1212  vector<NodeSharedPtr> tmp = orderedLoops[i];
1213  reverse(tmp.begin(), tmp.end());
1214  orderedLoops[i] = tmp;
1215  }
1216  }
1217 }
std::vector< std::vector< NodeSharedPtr > > orderedLoops
list of boundary nodes in their order loops
Definition: FaceMesh.h:158
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
std::map< int, CurveMeshSharedPtr > m_curvemeshes
Map of the curve meshes which bound the surfaces.
Definition: FaceMesh.h:151
int m_id
id of the surface mesh
Definition: FaceMesh.h:156
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
double NekDouble
std::vector< EdgeLoop > m_edgeloops
data structure containing the edges, their order and oreientation for the surface ...
Definition: FaceMesh.h:154
void Nektar::NekMeshUtils::FaceMesh::Smoothing ( )
private

performs node smoothing on face

Definition at line 210 of file FaceMesh.cpp.

References ASSERTL0, Nektar::LibUtilities::eQuadrilateral, and Nektar::iterator.

211 {
212  EdgeSet::iterator eit;
213  NodeSet::iterator nit;
214 
215  map<int, vector<EdgeSharedPtr> > connectingedges;
216 
217  map<int, vector<ElementSharedPtr> > connectingelements;
218 
219  for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
220  {
221  connectingedges[(*eit)->m_n1->m_id].push_back(*eit);
222  connectingedges[(*eit)->m_n2->m_id].push_back(*eit);
223  }
224 
225  for (int i = 0; i < m_localElements.size(); i++)
226  {
227  vector<NodeSharedPtr> v = m_localElements[i]->GetVertexList();
228  for (int j = 0; j < 3; j++)
229  {
230  connectingelements[v[j]->m_id].push_back(m_localElements[i]);
231  }
232  }
233 
234  // perform 4 runs of elastic relaxation based on the octree
235  for (int q = 0; q < 4; q++)
236  {
237  for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
238  {
239  if ((*nit)->GetNumCadCurve() > 0) // node is on curve so skip
240  continue;
241 
242  vector<NodeSharedPtr> connodes; // this can be real nodes or dummy
243  // nodes depending on the system
244 
245  vector<EdgeSharedPtr> edges = connectingedges[(*nit)->m_id];
246  vector<ElementSharedPtr> els = connectingelements[(*nit)->m_id];
247 
248  bool perfrom = true;
249  for (int i = 0; i < els.size(); i++)
250  {
251  if (els[i]->GetConf().m_e == LibUtilities::eQuadrilateral)
252  {
253  perfrom = false;
254  break;
255  }
256  }
257 
258  if (!perfrom)
259  continue;
260 
261  vector<NodeSharedPtr> nodesystem;
262  vector<NekDouble> lamp;
263 
264  for (int i = 0; i < edges.size(); i++)
265  {
266  vector<NekDouble> lambda;
267 
268  NodeSharedPtr J;
269  if (*nit == edges[i]->m_n1)
270  J = edges[i]->m_n2;
271  else if (*nit == edges[i]->m_n2)
272  J = edges[i]->m_n1;
273  else
274  ASSERTL0(false, "could not find node");
275 
276  Array<OneD, NekDouble> ui = (*nit)->GetCADSurfInfo(m_id);
277  Array<OneD, NekDouble> uj = J->GetCADSurfInfo(m_id);
278 
279  for (int j = 0; j < els.size(); j++)
280  {
281  vector<NodeSharedPtr> v = els[j]->GetVertexList();
282  if (v[0] == J || v[1] == J || v[2] == J)
283  continue; // elememt is adjacent to J therefore no
284  // intersection on IJ
285 
286  // need to find other edge
287  EdgeSharedPtr AtoB;
288  bool found = false;
289  vector<EdgeSharedPtr> es = els[j]->GetEdgeList();
290  for (int k = 0; k < es.size(); k++)
291  {
292  if (!(es[k]->m_n1 == *nit || es[k]->m_n2 == *nit))
293  {
294  found = true;
295  AtoB = es[k];
296  break;
297  }
298  }
299  ASSERTL0(found, "failed to find edge to test");
300 
301  Array<OneD, NekDouble> A = AtoB->m_n1->GetCADSurfInfo(m_id);
302  Array<OneD, NekDouble> B = AtoB->m_n2->GetCADSurfInfo(m_id);
303 
304  NekDouble lam = ((A[0] - uj[0]) * (B[1] - A[1]) -
305  (A[1] - uj[1]) * (B[0] - A[0])) /
306  ((ui[0] - uj[0]) * (B[1] - A[1]) -
307  (ui[1] - uj[1]) * (B[0] - A[0]));
308 
309  if (!(lam < 0) && !(lam > 1))
310  lambda.push_back(lam);
311  }
312 
313  if (lambda.size() > 0)
314  {
315  sort(lambda.begin(), lambda.end());
316  // make a new dummy node based on the system
317  Array<OneD, NekDouble> ud(2);
318  ud[0] = uj[0] + lambda[0] * (ui[0] - uj[0]);
319  ud[1] = uj[1] + lambda[0] * (ui[1] - uj[1]);
320  Array<OneD, NekDouble> locd = m_cadsurf->P(ud);
321  NodeSharedPtr dn = boost::shared_ptr<Node>(
322  new Node(0, locd[0], locd[1], locd[2]));
323  dn->SetCADSurf(m_id, m_cadsurf, ud);
324 
325  nodesystem.push_back(dn);
326  lamp.push_back(lambda[0]);
327  }
328  else
329  {
330  nodesystem.push_back(J);
331  lamp.push_back(1.0);
332  }
333  }
334 
335  Array<OneD, NekDouble> ui(2);
336  ui[0] = 0.0;
337  ui[1] = 0.0;
338 
339  for (int i = 0; i < nodesystem.size(); i++)
340  {
341  Array<OneD, NekDouble> uj = nodesystem[i]->GetCADSurfInfo(m_id);
342  ui[0] += uj[0] / nodesystem.size();
343  ui[1] += uj[1] / nodesystem.size();
344  }
345 
346  Array<OneD, NekDouble> bounds = m_cadsurf->GetBounds();
347 
348  Array<OneD, NekDouble> uvn(2);
349 
350  uvn[0] = ui[0];
351  uvn[1] = ui[1];
352 
353  if (!(uvn[0] < bounds[0] || uvn[0] > bounds[1] ||
354  uvn[1] < bounds[2] || uvn[1] > bounds[3]))
355  {
356  Array<OneD, NekDouble> l2 = m_cadsurf->P(uvn);
357  (*nit)->Move(l2, m_id, uvn);
358  }
359  }
360  }
361 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
std::vector< ElementSharedPtr > m_localElements
local list of elements
Definition: FaceMesh.h:170
EdgeSet m_localEdges
local set of edges
Definition: FaceMesh.h:168
int m_id
id of the surface mesh
Definition: FaceMesh.h:156
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
double NekDouble
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
CADSurfSharedPtr m_cadsurf
CAD surface.
Definition: FaceMesh.h:147
NodeSet m_localNodes
local set of nodes
Definition: FaceMesh.h:166
void Nektar::NekMeshUtils::FaceMesh::Stretching ( )
private

Calculate the paramter plane streching factor.

Definition at line 873 of file FaceMesh.cpp.

874 {
875  // define a sampling and calculate the aspect ratio of the paramter plane
876  m_str = 0.0;
877  Array<OneD, NekDouble> bnds = m_cadsurf->GetBounds();
878 
879  NekDouble dxu = int(bnds[1] - bnds[0] < bnds[3] - bnds[2]
880  ? 40
881  : (bnds[1] - bnds[0]) / (bnds[3] - bnds[2]) * 40);
882  NekDouble dxv = int(bnds[3] - bnds[2] < bnds[1] - bnds[0]
883  ? 40
884  : (bnds[3] - bnds[2]) / (bnds[1] - bnds[0]) * 40);
885 
886  NekDouble du = (bnds[1] - bnds[0]) / dxu;
887  NekDouble dv = (bnds[3] - bnds[2]) / dxv;
888 
889  int ct = 0;
890 
891  for (int i = 0; i < dxu; i++)
892  {
893  for (int j = 0; j < dxv; j++)
894  {
895  Array<OneD, NekDouble> uv(2);
896  uv[0] = bnds[0] + i * du;
897  uv[1] = bnds[2] + j * dv;
898  if (i == dxu - 1)
899  uv[0] = bnds[1];
900  if (j == dxv - 1)
901  uv[1] = bnds[3];
902  Array<OneD, NekDouble> r = m_cadsurf->D1(uv);
903 
904  NekDouble ru = sqrt(r[3] * r[3] + r[4] * r[4] + r[5] * r[5]);
905  NekDouble rv = sqrt(r[6] * r[6] + r[7] * r[7] + r[8] * r[8]);
906 
907  ru *= du;
908  rv *= dv;
909 
910  if (rv < 1E-8)
911  continue;
912 
913  m_str += ru / rv;
914  ct++;
915  }
916  }
917 
918  m_str /= ct;
919 }
double NekDouble
NekDouble m_str
pplane stretching
Definition: FaceMesh.h:162
CADSurfSharedPtr m_cadsurf
CAD surface.
Definition: FaceMesh.h:147
bool Nektar::NekMeshUtils::FaceMesh::Validate ( )
private

Validate the surface mesh base on the octree and real dimensions of the edges.

Definition at line 921 of file FaceMesh.cpp.

922 {
923  // check all edges in the current mesh for length against the octree
924  // if the octree is not conformed to add a new point inside the triangle
925  // if no new points are added meshing can stop
926  int pointBefore = m_stienerpoints.size();
927  for (int i = 0; i < m_connec.size(); i++)
928  {
929  Array<OneD, NekDouble> triDelta(3);
930 
931  Array<OneD, NekDouble> r(3);
932 
933  r[0] = m_connec[i][0]->Distance(m_connec[i][1]);
934  r[1] = m_connec[i][1]->Distance(m_connec[i][2]);
935  r[2] = m_connec[i][2]->Distance(m_connec[i][0]);
936 
937  triDelta[0] = m_octree->Query(m_connec[i][0]->GetLoc());
938  triDelta[1] = m_octree->Query(m_connec[i][1]->GetLoc());
939  triDelta[2] = m_octree->Query(m_connec[i][2]->GetLoc());
940 
941  int numValid = 0;
942 
943  if (r[0] < triDelta[0] && r[2] < triDelta[0])
944  numValid++;
945 
946  if (r[1] < triDelta[1] && r[0] < triDelta[1])
947  numValid++;
948 
949  if (r[2] < triDelta[2] && r[1] < triDelta[2])
950  numValid++;
951 
952  if (numValid != 3)
953  {
954  Array<OneD, NekDouble> ainfo, binfo, cinfo;
955  ainfo = m_connec[i][0]->GetCADSurfInfo(m_id);
956  binfo = m_connec[i][1]->GetCADSurfInfo(m_id);
957  cinfo = m_connec[i][2]->GetCADSurfInfo(m_id);
958 
959  Array<OneD, NekDouble> uvc(2);
960  uvc[0] = (ainfo[0] + binfo[0] + cinfo[0]) / 3.0;
961  uvc[1] = (ainfo[1] + binfo[1] + cinfo[1]) / 3.0;
962  AddNewPoint(uvc);
963  }
964  }
965 
966  if (m_stienerpoints.size() == pointBefore)
967  {
968  return false;
969  }
970  else
971  {
972  return true;
973  }
974 }
OctreeSharedPtr m_octree
Octree object.
Definition: FaceMesh.h:149
std::vector< std::vector< NodeSharedPtr > > m_connec
triangle connectiviities
Definition: FaceMesh.h:164
int m_id
id of the surface mesh
Definition: FaceMesh.h:156
std::vector< NodeSharedPtr > m_stienerpoints
list of stiener points in the triangulation
Definition: FaceMesh.h:160
void AddNewPoint(Array< OneD, NekDouble > uv)
adds a new stiener point to the triangulation for meshing
Definition: FaceMesh.cpp:976

Friends And Related Function Documentation

friend class MemoryManager< FaceMesh >
friend

Definition at line 56 of file FaceMesh.h.

Member Data Documentation

std::vector<std::pair<NodeSharedPtr, NodeSharedPtr> > Nektar::NekMeshUtils::FaceMesh::blpairs
private

list of node links between node on loop and its corresponding interior node in quads

Definition at line 177 of file FaceMesh.h.

NekDouble Nektar::NekMeshUtils::FaceMesh::m_bl
private

boundary layer thickness

Definition at line 172 of file FaceMesh.h.

CADSurfSharedPtr Nektar::NekMeshUtils::FaceMesh::m_cadsurf
private

CAD surface.

Definition at line 147 of file FaceMesh.h.

Referenced by FaceMesh().

std::vector<std::vector<NodeSharedPtr> > Nektar::NekMeshUtils::FaceMesh::m_connec
private

triangle connectiviities

Definition at line 164 of file FaceMesh.h.

std::map<int, CurveMeshSharedPtr> Nektar::NekMeshUtils::FaceMesh::m_curvemeshes
private

Map of the curve meshes which bound the surfaces.

Definition at line 151 of file FaceMesh.h.

std::vector<EdgeLoop> Nektar::NekMeshUtils::FaceMesh::m_edgeloops
private

data structure containing the edges, their order and oreientation for the surface

Definition at line 154 of file FaceMesh.h.

Referenced by FaceMesh().

int Nektar::NekMeshUtils::FaceMesh::m_id
private

id of the surface mesh

Definition at line 156 of file FaceMesh.h.

EdgeSet Nektar::NekMeshUtils::FaceMesh::m_localEdges
private

local set of edges

Definition at line 168 of file FaceMesh.h.

std::vector<ElementSharedPtr> Nektar::NekMeshUtils::FaceMesh::m_localElements
private

local list of elements

Definition at line 170 of file FaceMesh.h.

NodeSet Nektar::NekMeshUtils::FaceMesh::m_localNodes
private

local set of nodes

Definition at line 166 of file FaceMesh.h.

bool Nektar::NekMeshUtils::FaceMesh::m_makebl
private

should build boundary layer

Definition at line 174 of file FaceMesh.h.

Referenced by FaceMesh().

MeshSharedPtr Nektar::NekMeshUtils::FaceMesh::m_mesh
private

mesh pointer

Definition at line 145 of file FaceMesh.h.

OctreeSharedPtr Nektar::NekMeshUtils::FaceMesh::m_octree
private

Octree object.

Definition at line 149 of file FaceMesh.h.

std::vector<NodeSharedPtr> Nektar::NekMeshUtils::FaceMesh::m_stienerpoints
private

list of stiener points in the triangulation

Definition at line 160 of file FaceMesh.h.

NekDouble Nektar::NekMeshUtils::FaceMesh::m_str
private

pplane stretching

Definition at line 162 of file FaceMesh.h.

std::vector<std::vector<NodeSharedPtr> > Nektar::NekMeshUtils::FaceMesh::orderedLoops
private

list of boundary nodes in their order loops

Definition at line 158 of file FaceMesh.h.