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

class for octree More...

#include <Octree.h>

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

Public Member Functions

 Octree (MeshSharedPtr m)
 
void Process ()
 builds the octree based on curvature sampling and user defined spacing More...
 
NekDouble Query (Array< OneD, NekDouble > loc)
 once constructed queryies the octree based on x,y,z location to get a mesh spacing More...
 
NekDouble GetMinDelta ()
 returns the miminum spacing in the octree (for meshing purposes) More...
 
void SetParameters (NekDouble &min, NekDouble &max, NekDouble &ep)
 sets the parameters used for curvature sampling More...
 
void WriteOctree (std::string nm)
 populates the mesh m with a invalid hexahedral mesh based on the octree, used for visualisation More...
 
void Refinement (std::string nm)
 informs the octree there is a user defined spacing file More...
 

Private Member Functions

void SmoothAllOctants ()
 Smooths specification over all octants to a gradation criteria. More...
 
void CompileSourcePointList ()
 gets an optimum number of curvature sampling points and calculates the curavture at these points More...
 
void SubDivide ()
 Function which initiates and controls the subdivision process. More...
 
void SmoothSurfaceOctants ()
 Smooths specification over the surface encompasing octants to a gradation criteria. More...
 
void PropagateDomain ()
 takes the mesh specification from surface octants and progates that through the domain so all octants have a specification using gradiation crieteria More...
 
int CountElemt ()
 estimates the number of elements to be created in the mesh More...
 
NekDouble ddx (OctantSharedPtr i, OctantSharedPtr j)
 Calculates the difference in delta divided by the difference in location between two octants i and j. More...
 
bool VerifyNeigbours ()
 Looks over all leaf octants and checks that their neigbour assigments are valid. More...
 

Private Attributes

NekDouble m_minDelta
 minimum delta in the octree More...
 
NekDouble m_maxDelta
 maximum delta in the octree More...
 
NekDouble m_eps
 curavture sensivity paramter More...
 
Array< OneD, NekDoublem_centroid
 x,y,z location of the center of the octree More...
 
NekDouble m_dim
 physical size of the octree More...
 
std::vector< SPBaseSharedPtrm_SPList
 list of source points More...
 
std::vector< OctantSharedPtrm_octants
 list of leaf octants More...
 
OctantSharedPtr m_masteroct
 master octant for searching More...
 
int m_numoct
 number of octants made, used for id index More...
 
MeshSharedPtr m_mesh
 Mesh object. More...
 
std::string m_refinement
 
std::vector< linesourcem_lsources
 

Detailed Description

class for octree

This class contains the routines to generate and query a automatically generated set of mesh spacing parameters based on the CAD

Definition at line 107 of file Octree.h.

Constructor & Destructor Documentation

Nektar::NekMeshUtils::Octree::Octree ( MeshSharedPtr  m)
inline

Definition at line 111 of file Octree.h.

111  : m_mesh(m)
112  {
113  }
MeshSharedPtr m_mesh
Mesh object.
Definition: Octree.h:236

Member Function Documentation

void Nektar::NekMeshUtils::Octree::CompileSourcePointList ( )
private

gets an optimum number of curvature sampling points and calculates the curavture at these points

TODO add extra source points from the line souce to the octree

Definition at line 837 of file Octree.cpp.

References Nektar::LibUtilities::PrintProgressbar().

838 {
839 
840  if(m_mesh->m_cad->Is2D())
841  {
842  for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
843  {
844  CADCurveSharedPtr curve = m_mesh->m_cad->GetCurve(i);
845  Array<OneD, NekDouble> bds = curve->GetBounds();
846  int samples = 100;
847  NekDouble dt = (bds[1] - bds[0]) / (samples + 1);
848  for (int j = 1; j < samples - 1; j++) // dont want first and last point
849  {
850  NekDouble t = bds[0] + dt * j;
851  NekDouble C = curve->Curvature(t);
852 
853  Array<OneD, NekDouble> loc = curve->P(t);
854 
855  vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > ss =
856  curve->GetAdjSurf();
857  Array<OneD, NekDouble> uv = ss[0].first->locuv(loc);
858 
859  if (C != 0.0)
860  {
861  NekDouble del = 2.0 * (1.0 / C) * sqrt(m_eps * (2.0 - m_eps));
862 
863  if (del > m_maxDelta)
864  {
865  del = m_maxDelta;
866  }
867  if (del < m_minDelta)
868  {
869  del = m_minDelta;
870  }
871 
872  CPointSharedPtr newCPoint =
874  ss[0].first->GetId(), uv, loc, del);
875 
876  m_SPList.push_back(newCPoint);
877  }
878  else
879  {
880  BPointSharedPtr newBPoint =
882  ss[0].first->GetId(), uv, loc);
883 
884  m_SPList.push_back(newBPoint);
885  }
886  }
887  }
888  }
889 
890  for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
891  {
892  if (m_mesh->m_verbose)
893  {
894  LibUtilities::PrintProgressbar(i, m_mesh->m_cad->GetNumSurf(),
895  "\tCompiling source points");
896  }
897 
898  CADSurfSharedPtr surf = m_mesh->m_cad->GetSurf(i);
899  Array<OneD, NekDouble> bounds = surf->GetBounds();
900 
901  // to figure out the amount of curvature sampling to be conducted on
902  // each parameter plane the surface is first sampled with a 40x40 grid
903  // the real space lengths of this grid are analysed to find the largest
904  // strecthing in the u and v directions
905  // this stretching is then cosnidered with the mindelta user input
906  // to find a number of sampling points in each direction which
907  // enures that in the final octree each surface octant will have at
908  // least one sample point within its volume.
909  // the 40x40 grid is used to ensure each surface has a minimum of 40x40
910  // samples.
911  NekDouble du = (bounds[1] - bounds[0]) / (40 - 1);
912  NekDouble dv = (bounds[3] - bounds[2]) / (40 - 1);
913 
914  NekDouble DeltaU = 0.0;
915  NekDouble DeltaV = 0.0;
916 
917  Array<TwoD, Array<OneD, NekDouble> > samplepoints(40, 40);
918 
919  for (int j = 0; j < 40; j++)
920  {
921  for (int k = 0; k < 40; k++)
922  {
923  Array<OneD, NekDouble> uv(2);
924  uv[0] = k * du + bounds[0];
925  uv[1] = j * dv + bounds[2];
926  if (j == 40 - 1)
927  uv[1] = bounds[3];
928  if (k == 40 - 1)
929  uv[0] = bounds[1];
930  samplepoints[k][j] = surf->P(uv);
931  }
932  }
933 
934  for (int j = 0; j < 40 - 1; j++)
935  {
936  for (int k = 0; k < 40 - 1; k++)
937  {
938  NekDouble deltau = sqrt(
939  (samplepoints[k][j][0] - samplepoints[k + 1][j][0]) *
940  (samplepoints[k][j][0] - samplepoints[k + 1][j][0]) +
941  (samplepoints[k][j][1] - samplepoints[k + 1][j][1]) *
942  (samplepoints[k][j][1] - samplepoints[k + 1][j][1]) +
943  (samplepoints[k][j][2] - samplepoints[k + 1][j][2]) *
944  (samplepoints[k][j][2] - samplepoints[k + 1][j][2]));
945  NekDouble deltav = sqrt(
946  (samplepoints[k][j][0] - samplepoints[k][j + 1][0]) *
947  (samplepoints[k][j][0] - samplepoints[k][j + 1][0]) +
948  (samplepoints[k][j][1] - samplepoints[k][j + 1][1]) *
949  (samplepoints[k][j][1] - samplepoints[k][j + 1][1]) +
950  (samplepoints[k][j][2] - samplepoints[k][j + 1][2]) *
951  (samplepoints[k][j][2] - samplepoints[k][j + 1][2]));
952 
953  if (deltau > DeltaU)
954  DeltaU = deltau;
955  if (deltav > DeltaV)
956  DeltaV = deltav;
957  }
958  }
959 
960  // these are the acutal number of sample points in each parametric
961  // direction
962  int nu = ceil(DeltaU / m_minDelta) * 40 * 2;
963  int nv = ceil(DeltaV / m_minDelta) * 40 * 2;
964 
965  for (int j = 0; j < nu; j++)
966  {
967  for (int k = 0; k < nv; k++)
968  {
969  Array<OneD, NekDouble> uv(2);
970  uv[0] = (bounds[1] - bounds[0]) / (nu - 1) * j + bounds[0];
971  uv[1] = (bounds[3] - bounds[2]) / (nv - 1) * k + bounds[2];
972 
973  // this prevents round off error at the end of the surface
974  // may not be neseercary but works
975  if (j == nu - 1)
976  uv[0] = bounds[1];
977  if (k == nv - 1)
978  uv[1] = bounds[3];
979 
980  NekDouble C = surf->Curvature(uv);
981 
982  // create new point based on smallest R, flat surfaces have k=0
983  // but still need a point for element estimation
984  if (C != 0.0)
985  {
986  NekDouble del =
987  2.0 * (1.0 / C) * sqrt(m_eps * (2.0 - m_eps));
988 
989  if (del > m_maxDelta)
990  {
991  del = m_maxDelta;
992  }
993  if (del < m_minDelta)
994  {
995  del = m_minDelta;
996  }
997 
998  CPointSharedPtr newCPoint =
1000  surf->GetId(), uv, surf->P(uv), del);
1001 
1002  m_SPList.push_back(newCPoint);
1003  }
1004  else
1005  {
1006  BPointSharedPtr newBPoint =
1008  surf->GetId(), uv, surf->P(uv));
1009 
1010  m_SPList.push_back(newBPoint);
1011  }
1012  }
1013  }
1014  }
1015  if (m_mesh->m_verbose)
1016  {
1017  cout << endl;
1018  }
1019 
1020  if (m_refinement.size() > 0)
1021  {
1022  if (m_mesh->m_verbose)
1023  {
1024  cout << "\t\tModifying based on refinement lines" << endl;
1025  }
1026  // now deal with the user defined spacing
1027  vector<string> lines;
1028 
1029  boost::split(lines, m_refinement, boost::is_any_of(":"));
1030 
1031  for (int i = 0; i < lines.size(); i++)
1032  {
1033  vector<NekDouble> data;
1034  ParseUtils::GenerateUnOrderedVector(lines[i].c_str(), data);
1035 
1036  Array<OneD, NekDouble> x1(3), x2(3);
1037 
1038  x1[0] = data[0];
1039  x1[1] = data[1];
1040  x1[2] = data[2];
1041  x2[0] = data[3];
1042  x2[1] = data[4];
1043  x2[2] = data[5];
1044 
1045  m_lsources.push_back(linesource(x1, x2, data[6], data[7]));
1046  }
1047 
1048  // this takes any existing sourcepoints within the influence range
1049  // and modifies them
1050  /*for (int i = 0; i < m_SPList.size(); i++)
1051  {
1052  for (int j = 0; j < m_lsources.size(); j++)
1053  {
1054  if (m_lsources[j].withinRange(m_SPList[i]->GetLoc()))
1055  {
1056  if(m_SPList[i]->GetType() == ePBoundary)
1057  {
1058  BPointSharedPtr bp =
1059  boost::dynamic_pointer_cast<BPoint>
1060  (m_SPList[i]);
1061 
1062  m_SPList[i] = bp->ChangeType();
1063 
1064  }
1065  m_SPList[i]->SetDelta(m_lsources[j].delta);
1066  }
1067  }
1068  }*/
1069  /// TODO add extra source points from the line souce to the octree
1070  }
1071 }
std::vector< SPBaseSharedPtr > m_SPList
list of source points
Definition: Octree.h:228
int PrintProgressbar(const int position, const int goal, const string message, int lastprogress=-1)
Prints a progressbar.
Definition: Progressbar.hpp:69
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
NekDouble m_minDelta
minimum delta in the octree
Definition: Octree.h:218
std::string m_refinement
Definition: Octree.h:238
double NekDouble
std::vector< linesource > m_lsources
Definition: Octree.h:239
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:172
NekDouble m_eps
curavture sensivity paramter
Definition: Octree.h:222
MeshSharedPtr m_mesh
Mesh object.
Definition: Octree.h:236
boost::shared_ptr< BPoint > BPointSharedPtr
boost::shared_ptr< CPoint > CPointSharedPtr
static bool GenerateUnOrderedVector(const char *const str, std::vector< NekDouble > &vec)
Definition: ParseUtils.hpp:128
boost::shared_ptr< CADCurve > CADCurveSharedPtr
Definition: CADCurve.h:194
NekDouble m_maxDelta
maximum delta in the octree
Definition: Octree.h:220
int Nektar::NekMeshUtils::Octree::CountElemt ( )
private

estimates the number of elements to be created in the mesh

Definition at line 725 of file Octree.cpp.

References Nektar::NekMeshUtils::eBack, Nektar::NekMeshUtils::eDown, Nektar::NekMeshUtils::eForward, Nektar::NekMeshUtils::eInside, Nektar::NekMeshUtils::eLeft, Nektar::NekMeshUtils::eOnBoundary, Nektar::NekMeshUtils::eRight, and Nektar::NekMeshUtils::eUp.

726 {
727  // by considering the volume of a tet evaluate the number of elements in the
728  // mesh
729 
730  NekDouble total = 0.0;
731 
732  Array<OneD, NekDouble> boundingBox = m_mesh->m_cad->GetBoundingBox();
733 
734  for (int i = 0; i < m_octants.size(); i++)
735  {
736  OctantSharedPtr oct = m_octants[i];
737  if (oct->GetLocation() == eInside)
738  {
739  total += 8.0 * oct->DX() * oct->DX() * oct->DX() /
740  (oct->GetDelta() * oct->GetDelta() * oct->GetDelta() /
741  6.0 / sqrt(2));
742  }
743  else if (oct->GetLocation() == eOnBoundary)
744  {
745  NekDouble vol = 1.0;
746  if (oct->FX(eBack) < boundingBox[1] &&
747  oct->FX(eForward) > boundingBox[0])
748  {
749  // then there is some over lap in x
750  NekDouble min, max;
751  if (oct->FX(eBack) > boundingBox[0])
752  {
753  min = oct->FX(eBack);
754  }
755  else
756  {
757  min = boundingBox[0];
758  }
759  if (boundingBox[1] < oct->FX(eForward))
760  {
761  max = boundingBox[1];
762  }
763  else
764  {
765  max = oct->FX(eForward);
766  }
767  vol *= (max - min);
768  }
769  else
770  {
771  vol *= 0.0;
772  }
773 
774  if (oct->FX(eDown) < boundingBox[3] &&
775  oct->FX(eUp) > boundingBox[2])
776  {
777  // then there is some over lap in x
778  NekDouble min, max;
779  if (oct->FX(eDown) > boundingBox[2])
780  {
781  min = oct->FX(eDown);
782  }
783  else
784  {
785  min = boundingBox[2];
786  }
787  if (boundingBox[3] < oct->FX(eUp))
788  {
789  max = boundingBox[3];
790  }
791  else
792  {
793  max = oct->FX(eUp);
794  }
795  vol *= (max - min);
796  }
797  else
798  {
799  vol *= 0.0;
800  }
801 
802  if (oct->FX(eRight) < boundingBox[5] &&
803  oct->FX(eLeft) > boundingBox[4])
804  {
805  // then there is some over lap in x
806  NekDouble min, max;
807  if (oct->FX(eRight) > boundingBox[4])
808  {
809  min = oct->FX(eRight);
810  }
811  else
812  {
813  min = boundingBox[4];
814  }
815  if (boundingBox[5] < oct->FX(eLeft))
816  {
817  max = boundingBox[5];
818  }
819  else
820  {
821  max = oct->FX(eLeft);
822  }
823  vol *= (max - min);
824  }
825  else
826  {
827  vol *= 0.0;
828  }
829  total += vol / 2.0 / (oct->GetDelta() * oct->GetDelta() *
830  oct->GetDelta() / 6.0 / sqrt(2));
831  }
832  }
833 
834  return int(total);
835 }
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:230
double NekDouble
boost::shared_ptr< Octant > OctantSharedPtr
Definition: Octant.h:74
MeshSharedPtr m_mesh
Mesh object.
Definition: Octree.h:236
NekDouble Nektar::NekMeshUtils::Octree::ddx ( OctantSharedPtr  i,
OctantSharedPtr  j 
)
private

Calculates the difference in delta divided by the difference in location between two octants i and j.

Definition at line 1073 of file Octree.cpp.

1074 {
1075  return fabs(i->GetDelta() - j->GetDelta()) / i->Distance(j);
1076 }
NekDouble Nektar::NekMeshUtils::Octree::GetMinDelta ( )

returns the miminum spacing in the octree (for meshing purposes)

Returns
miminum delta in octree

Definition at line 191 of file Octree.cpp.

192 {
193  NekDouble tmp = numeric_limits<double>::max();
194 
195  for (int i = 0; i < m_lsources.size(); i++)
196  {
197  tmp = min(m_lsources[i].delta, tmp);
198  }
199  return min(m_minDelta, tmp);
200 }
NekDouble m_minDelta
minimum delta in the octree
Definition: Octree.h:218
double NekDouble
std::vector< linesource > m_lsources
Definition: Octree.h:239
void Nektar::NekMeshUtils::Octree::Process ( )

builds the octree based on curvature sampling and user defined spacing

Definition at line 52 of file Octree.cpp.

53 {
54  Array<OneD, NekDouble> boundingBox = m_mesh->m_cad->GetBoundingBox();
55 
56  // build curvature samples
58 
59  if (m_mesh->m_verbose)
60  {
61  cout << "\tCurvature samples: " << m_SPList.size() << endl;
62  }
63 
64  // make master octant based on the bounding box of the domain
65  m_dim = max((boundingBox[1] - boundingBox[0]) / 2.0,
66  (boundingBox[3] - boundingBox[2]) / 2.0);
67 
68  m_dim = max(m_dim, (boundingBox[5] - boundingBox[4]) / 2.0);
69 
70  m_centroid = Array<OneD, NekDouble>(3);
71  m_centroid[0] = (boundingBox[1] + boundingBox[0]) / 2.0;
72  m_centroid[1] = (boundingBox[3] + boundingBox[2]) / 2.0;
73  m_centroid[2] = (boundingBox[5] + boundingBox[4]) / 2.0;
74 
76  0, m_centroid[0], m_centroid[1], m_centroid[2], m_dim, m_SPList);
77 
78  // begin recersive subdivision
79  SubDivide();
80 
81  m_octants.clear();
82  m_masteroct->CompileLeaves(m_octants);
83 
84  if (m_mesh->m_verbose)
85  {
86  cout << "\tOctants: " << m_octants.size() << endl;
87  }
88 
90 
92 
94 
95  if (m_mesh->m_verbose)
96  {
97  int elem = CountElemt();
98  cout << "\tPredicted mesh: " << elem << endl;
99  }
100 }
std::vector< SPBaseSharedPtr > m_SPList
list of source points
Definition: Octree.h:228
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void SubDivide()
Function which initiates and controls the subdivision process.
Definition: Octree.cpp:269
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:230
OctantSharedPtr m_masteroct
master octant for searching
Definition: Octree.h:232
void PropagateDomain()
takes the mesh specification from surface octants and progates that through the domain so all octants...
Definition: Octree.cpp:516
int CountElemt()
estimates the number of elements to be created in the mesh
Definition: Octree.cpp:725
void CompileSourcePointList()
gets an optimum number of curvature sampling points and calculates the curavture at these points ...
Definition: Octree.cpp:837
NekDouble m_dim
physical size of the octree
Definition: Octree.h:226
Array< OneD, NekDouble > m_centroid
x,y,z location of the center of the octree
Definition: Octree.h:224
MeshSharedPtr m_mesh
Mesh object.
Definition: Octree.h:236
void SmoothSurfaceOctants()
Smooths specification over the surface encompasing octants to a gradation criteria.
Definition: Octree.cpp:459
void SmoothAllOctants()
Smooths specification over all octants to a gradation criteria.
Definition: Octree.cpp:676
void Nektar::NekMeshUtils::Octree::PropagateDomain ( )
private

takes the mesh specification from surface octants and progates that through the domain so all octants have a specification using gradiation crieteria

Definition at line 516 of file Octree.cpp.

References ASSERTL0, Nektar::NekMeshUtils::eInside, Nektar::NekMeshUtils::eOnBoundary, Nektar::NekMeshUtils::eOutside, Nektar::SpatialDomains::eUnknown, and Nektar::iterator.

517 {
518  // until all m_octants know their delta specifcation and orientaion
519  // look over all m_octants and if their neighours know either their
520  // orientation
521  // or specifcation calculate one for this octant
522  int ct = 0;
523 
524  do
525  {
526  ct = 0;
527  for (int i = 0; i < m_octants.size(); i++)
528  {
529  OctantSharedPtr oct = m_octants[i];
530 
531  if (!oct->IsDeltaKnown())
532  { // if delta has not been asigned
533  vector<OctantSharedPtr> known;
534  map<OctantFace, vector<OctantSharedPtr> > nList =
535  oct->GetNeigbours();
536  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
537 
538  for (it = nList.begin(); it != nList.end(); it++)
539  {
540  for (int j = 0; j < it->second.size(); j++)
541  {
542  if (it->second[j]->IsDeltaKnown())
543  {
544  known.push_back(it->second[j]);
545  }
546  }
547  }
548 
549  if (known.size() > 0)
550  {
551  vector<NekDouble> deltaPrime;
552  for (int j = 0; j < known.size(); j++)
553  {
554  NekDouble r = oct->Distance(known[j]);
555 
556  if (0.199 * r + known[j]->GetDelta() < m_maxDelta)
557  {
558  deltaPrime.push_back(0.199 * r +
559  known[j]->GetDelta());
560  }
561  else
562  {
563  deltaPrime.push_back(m_maxDelta);
564  }
565  }
566  NekDouble min = numeric_limits<double>::max();
567  for (int j = 0; j < deltaPrime.size(); j++)
568  {
569  if (deltaPrime[j] < min)
570  {
571  min = deltaPrime[j];
572  }
573  }
574  oct->SetDelta(min);
575  ct += 1;
576  deltaPrime.clear();
577  }
578  known.clear();
579  }
580 
581  if (oct->GetLocation() == eUnknown)
582  { // if the node does not know its location
583  vector<OctantSharedPtr> known;
584  map<OctantFace, vector<OctantSharedPtr> > nList =
585  oct->GetNeigbours();
586  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
587 
588  for (it = nList.begin(); it != nList.end(); it++)
589  {
590  for (int j = 0; j < it->second.size(); j++)
591  {
592  if (it->second[j]->GetLocation() != eUnknown)
593  {
594  known.push_back(it->second[j]);
595  }
596  }
597  }
598 
599  if (known.size() > 0)
600  {
601  vector<OctantSharedPtr> isNotOnBound;
602  for (int j = 0; j < known.size(); j++)
603  {
604  if (known[j]->GetLocation() != eOnBoundary)
605  {
606  isNotOnBound.push_back(known[j]);
607  }
608  }
609 
610  if (isNotOnBound.size() > 0)
611  {
612  oct->SetLocation(isNotOnBound[0]->GetLocation());
613  }
614  else
615  {
616  NekDouble dist = numeric_limits<double>::max();
617 
618  OctantSharedPtr closest;
619 
620  bool f = false;
621  for (int j = 0; j < known.size(); j++)
622  {
623  if (oct->Distance(known[j]) < dist)
624  {
625  closest = known[j];
626  dist = oct->Distance(known[j]);
627  f = true;
628  }
629  }
630  ASSERTL0(f, "closest never set");
631 
632  SPBaseSharedPtr sp = closest->GetABoundPoint();
633 
634  Array<OneD, NekDouble> octloc, sploc, vec(3), uv, N;
635  int surf;
636  sp->GetCAD(surf, uv);
637  N = m_mesh->m_cad->GetSurf(surf)->N(uv);
638 
639  octloc = oct->GetLoc();
640  sploc = sp->GetLoc();
641 
642  vec[0] = octloc[0] - sploc[0];
643  vec[1] = octloc[1] - sploc[1];
644  vec[2] = octloc[2] - sploc[2];
645 
646  NekDouble dot =
647  vec[0] * N[0] + vec[1] * N[1] + vec[2] * N[2];
648 
649  if (dot <= 0.0)
650  {
651  oct->SetLocation(eOutside);
652  ct += 1;
653  }
654  else
655  {
656  oct->SetLocation(eInside);
657  ct += 1;
658  }
659  }
660  }
661  known.clear();
662  }
663  }
664 
665  } while (ct > 0);
666 
667  for (int i = 0; i < m_octants.size(); i++)
668  {
669  if (!m_octants[i]->IsDeltaKnown())
670  {
671  m_octants[i]->SetDelta(m_maxDelta);
672  }
673  }
674 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:230
boost::shared_ptr< SPBase > SPBaseSharedPtr
double NekDouble
boost::shared_ptr< Octant > OctantSharedPtr
Definition: Octant.h:74
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
MeshSharedPtr m_mesh
Mesh object.
Definition: Octree.h:236
NekDouble m_maxDelta
maximum delta in the octree
Definition: Octree.h:220
NekDouble Nektar::NekMeshUtils::Octree::Query ( Array< OneD, NekDouble loc)

once constructed queryies the octree based on x,y,z location to get a mesh spacing

Parameters
locarray of x,y,z
Returns
mesh spacing parameter

Definition at line 102 of file Octree.cpp.

References ASSERTL0.

103 {
104  // starting at master octant 0 move through succsesive m_octants which
105  // contain the point loc until a leaf is found
106  // first search through sourcepoints
107 
108  NekDouble tmp = numeric_limits<double>::max();
109 
110  for (int i = 0; i < m_lsources.size(); i++)
111  {
112  if (m_lsources[i].withinRange(loc))
113  {
114  tmp = min(m_lsources[i].delta, tmp);
115  }
116  }
117 
119  int quad;
120 
121  bool found = false;
122 
123  while (!found)
124  {
125  Array<OneD, NekDouble> octloc = n->GetLoc();
126 
127  if (!(loc[0] < octloc[0]) && // forward
128  !(loc[1] < octloc[1]) && // forward
129  !(loc[2] < octloc[2])) // forward
130  {
131  quad = 0;
132  }
133  else if (!(loc[0] < octloc[0]) && // forward
134  !(loc[1] < octloc[1]) && // forward
135  !(loc[2] > octloc[2])) // back
136  {
137  quad = 1;
138  }
139  else if (!(loc[0] < octloc[0]) && // forward
140  !(loc[1] > octloc[1]) && // back
141  !(loc[2] < octloc[2])) // forward
142  {
143  quad = 2;
144  }
145  else if (!(loc[0] < octloc[0]) && // forward
146  !(loc[1] > octloc[1]) && // back
147  !(loc[2] > octloc[2])) // back
148  {
149  quad = 3;
150  }
151  else if (!(loc[0] > octloc[0]) && // back
152  !(loc[1] < octloc[1]) && // forward
153  !(loc[2] < octloc[2])) // forward
154  {
155  quad = 4;
156  }
157  else if (!(loc[0] > octloc[0]) && // back
158  !(loc[1] < octloc[1]) && // forward
159  !(loc[2] > octloc[2])) // back
160  {
161  quad = 5;
162  }
163  else if (!(loc[0] > octloc[0]) && // back
164  !(loc[1] > octloc[1]) && // back
165  !(loc[2] < octloc[2])) // forward
166  {
167  quad = 6;
168  }
169  else if (!(loc[0] > octloc[0]) && // back
170  !(loc[1] > octloc[1]) && // back
171  !(loc[2] > octloc[2])) // back
172  {
173  quad = 7;
174  }
175  else
176  {
177  ASSERTL0(false, "Cannot locate quadrant");
178  }
179 
180  n = n->GetChild(quad);
181 
182  if (n->IsLeaf())
183  {
184  found = true;
185  }
186  }
187 
188  return min(n->GetDelta(), tmp);
189 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
OctantSharedPtr m_masteroct
master octant for searching
Definition: Octree.h:232
double NekDouble
boost::shared_ptr< Octant > OctantSharedPtr
Definition: Octant.h:74
std::vector< linesource > m_lsources
Definition: Octree.h:239
void Nektar::NekMeshUtils::Octree::Refinement ( std::string  nm)
inline

informs the octree there is a user defined spacing file

Parameters
nmname of the user defined spacing file

Definition at line 163 of file Octree.h.

References m_refinement.

164  {
165  m_refinement = nm;
166  }
std::string m_refinement
Definition: Octree.h:238
void Nektar::NekMeshUtils::Octree::SetParameters ( NekDouble min,
NekDouble max,
NekDouble ep 
)
inline

sets the parameters used for curvature sampling

Parameters
minminimum spacing to be found in the mesh
maxmaximum spacing to be found in the mesh
epscurvature sensivity relating radius of curvature to spacing

Definition at line 144 of file Octree.h.

References m_eps, m_maxDelta, and m_minDelta.

145  {
146  m_minDelta = min;
147  m_maxDelta = max;
148  m_eps = ep;
149  }
NekDouble m_minDelta
minimum delta in the octree
Definition: Octree.h:218
NekDouble m_eps
curavture sensivity paramter
Definition: Octree.h:222
NekDouble m_maxDelta
maximum delta in the octree
Definition: Octree.h:220
void Nektar::NekMeshUtils::Octree::SmoothAllOctants ( )
private

Smooths specification over all octants to a gradation criteria.

Definition at line 676 of file Octree.cpp.

References Nektar::iterator.

677 {
678  // until no more changes occur smooth the mesh specification between all
679  // m_octants not particualrly strictly
680  int ct = 0;
681 
682  do
683  {
684  ct = 0;
685  for (int i = 0; i < m_octants.size(); i++)
686  {
687  OctantSharedPtr oct = m_octants[i];
688 
689  vector<OctantSharedPtr> check;
690  map<OctantFace, vector<OctantSharedPtr> > nList =
691  oct->GetNeigbours();
692  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
693  for (it = nList.begin(); it != nList.end(); it++)
694  {
695  for (int j = 0; j < it->second.size(); j++)
696  {
697  if (it->second[j]->GetDelta() < oct->GetDelta() &&
698  ddx(oct, it->second[j]) > 0.2)
699  {
700  check.push_back(it->second[j]);
701  }
702  }
703  }
704 
705  if (check.size() > 0)
706  {
707  NekDouble deltaSM = numeric_limits<double>::max();
708  for (int j = 0; j < check.size(); j++)
709  {
710  NekDouble r = oct->Distance(check[j]);
711 
712  if (0.199 * r + check[j]->GetDelta() < deltaSM)
713  {
714  deltaSM = 0.199 * r + check[j]->GetDelta();
715  }
716  }
717  oct->SetDelta(deltaSM);
718  ct += 1;
719  }
720  }
721 
722  } while (ct > 0);
723 }
NekDouble ddx(OctantSharedPtr i, OctantSharedPtr j)
Calculates the difference in delta divided by the difference in location between two octants i and j...
Definition: Octree.cpp:1073
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:230
double NekDouble
boost::shared_ptr< Octant > OctantSharedPtr
Definition: Octant.h:74
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void Nektar::NekMeshUtils::Octree::SmoothSurfaceOctants ( )
private

Smooths specification over the surface encompasing octants to a gradation criteria.

Definition at line 459 of file Octree.cpp.

References Nektar::iterator.

460 {
461  // for all the m_octants which are surface containing and know their delta
462  // specification, look over all neighbours and ensure the specification
463  // between them is smooth
464  int ct = 0;
465 
466  do
467  {
468  ct = 0;
469  for (int i = 0; i < m_octants.size(); i++)
470  {
471  OctantSharedPtr oct = m_octants[i];
472 
473  if (oct->IsDeltaKnown())
474  {
475  vector<OctantSharedPtr> check;
476  map<OctantFace, vector<OctantSharedPtr> > nList =
477  oct->GetNeigbours();
478  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
479 
480  for (it = nList.begin(); it != nList.end(); it++)
481  {
482  for (int j = 0; j < it->second.size(); j++)
483  {
484  if (it->second[j]->IsDeltaKnown() &&
485  it->second[j]->GetDelta() < oct->GetDelta() &&
486  ddx(oct, it->second[j]) > 0.2)
487  {
488  check.push_back(it->second[j]);
489  }
490  }
491  }
492 
493  // for each neighbour listed in check_id, figure out the
494  // smoothed
495  // delta, and asign the miminum of these to nodes[i].GetDelta()
496  if (check.size() > 0)
497  {
498  NekDouble deltaSM = numeric_limits<double>::max();
499  for (int j = 0; j < check.size(); j++)
500  {
501  NekDouble r = oct->Distance(check[j]);
502 
503  if (0.199 * r + check[j]->GetDelta() < deltaSM)
504  {
505  deltaSM = 0.199 * r + check[j]->GetDelta();
506  }
507  }
508  oct->SetDelta(deltaSM);
509  ct += 1;
510  }
511  }
512  }
513  } while (ct > 0);
514 }
NekDouble ddx(OctantSharedPtr i, OctantSharedPtr j)
Calculates the difference in delta divided by the difference in location between two octants i and j...
Definition: Octree.cpp:1073
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:230
double NekDouble
boost::shared_ptr< Octant > OctantSharedPtr
Definition: Octant.h:74
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void Nektar::NekMeshUtils::Octree::SubDivide ( )
private

Function which initiates and controls the subdivision process.

Definition at line 269 of file Octree.cpp.

References Nektar::iterator.

270 {
271  bool repeat;
272  int ct = 1;
273 
274  m_numoct = 1;
275  m_masteroct->Subdivide(m_masteroct, m_numoct);
276 
277  if (m_mesh->m_verbose)
278  {
279  cout << "\tSubdivide iteration: ";
280  }
281 
282  do
283  {
284  if (m_mesh->m_verbose)
285  {
286  cout << "\r ";
287  cout << "\r";
288  cout << "\tSubdivide iteration: " << ct;
289  cout.flush();
290  }
291  ct++;
292  repeat = false;
293  m_octants.clear();
294  // grab a list of the leaves curently in the octree
295  m_masteroct->CompileLeaves(m_octants);
296 
297  VerifyNeigbours();
298 
299  // neeed to create a divide list, in the first list will be m_octants
300  // which need to
301  // subdivide based on curvature,
302  // in the next list will be ocants which need to subdivide to make sure
303  // level criteria is statisified for the previous list and so on
304  // the list will keep building till no more lists are required.
305  // the list will then be iterated through backwards to subdivide all the
306  // sublists in turn.
307  vector<vector<OctantSharedPtr> > dividelist;
308  set<int> inlist;
309  // build initial list
310  {
311  vector<OctantSharedPtr> sublist;
312  for (int i = 0; i < m_octants.size(); i++)
313  {
314  if (m_octants[i]->NeedDivide() &&
315  m_octants[i]->DX() / 4.0 > m_minDelta)
316  {
317  sublist.push_back(m_octants[i]);
318  inlist.insert(m_octants[i]->GetId());
319  repeat = true; // if there is a subdivision this whole
320  // process needs to be repeated
321  }
322  }
323  dividelist.push_back(sublist);
324  }
325  // then loop over building sublists until no more are required
326  int ct2 = 0;
327  while (true)
328  {
329  ct2++;
330  vector<OctantSharedPtr> newsublist,
331  previouslist = dividelist.back();
332  for (int i = 0; i < previouslist.size(); i++)
333  {
334  map<OctantFace, vector<OctantSharedPtr> > nlist =
335  previouslist[i]->GetNeigbours();
336  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
337  for (it = nlist.begin(); it != nlist.end(); it++)
338  {
339  for (int j = 0; j < it->second.size(); j++)
340  {
341  if (previouslist[i]->DX() < it->second[j]->DX())
342  {
344  inlist.find(it->second[j]->GetId());
345  if (s == inlist.end())
346  {
347  inlist.insert(it->second[j]->GetId());
348  newsublist.push_back(it->second[j]);
349  }
350  }
351  }
352  }
353  }
354  if (newsublist.size() == 0)
355  {
356  break;
357  }
358  else
359  {
360  dividelist.push_back(newsublist);
361  }
362  }
363 
364  vector<vector<OctantSharedPtr> >::reverse_iterator rit;
365  for (rit = dividelist.rbegin(); rit != dividelist.rend(); rit++)
366  {
367  vector<OctantSharedPtr> currentlist = *rit;
368  for (int i = 0; i < currentlist.size(); i++)
369  {
370  currentlist[i]->Subdivide(currentlist[i], m_numoct);
371  }
372  }
373  } while (repeat);
374 
375  if (m_mesh->m_verbose)
376  {
377  cout << endl;
378  }
379 }
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:230
NekDouble m_minDelta
minimum delta in the octree
Definition: Octree.h:218
bool VerifyNeigbours()
Looks over all leaf octants and checks that their neigbour assigments are valid.
Definition: Octree.cpp:381
OctantSharedPtr m_masteroct
master octant for searching
Definition: Octree.h:232
int m_numoct
number of octants made, used for id index
Definition: Octree.h:234
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
MeshSharedPtr m_mesh
Mesh object.
Definition: Octree.h:236
bool Nektar::NekMeshUtils::Octree::VerifyNeigbours ( )
private

Looks over all leaf octants and checks that their neigbour assigments are valid.

Definition at line 381 of file Octree.cpp.

References Nektar::NekMeshUtils::eBack, Nektar::NekMeshUtils::eDown, Nektar::NekMeshUtils::eForward, Nektar::NekMeshUtils::eLeft, Nektar::NekMeshUtils::eRight, Nektar::NekMeshUtils::eUp, and Nektar::iterator.

382 {
383  // check all octant links to their neighbours
384  // at all times in the subdivision the set of neigbours must
385  // conform to a set of criteria such as smoothness in size
386  // this checks that
387  bool error = false;
388  for (int i = 0; i < m_octants.size(); i++)
389  {
390  bool valid = true;
391  map<OctantFace, vector<OctantSharedPtr> > nlist =
392  m_octants[i]->GetNeigbours();
393  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
394  for (it = nlist.begin(); it != nlist.end(); it++)
395  {
396  if (it->second.size() == 0)
397  {
398  NekDouble expectedfx = 0.0;
399  switch (it->first)
400  {
401  case eUp:
402  expectedfx = m_centroid[1] + m_dim;
403  break;
404  case eDown:
405  expectedfx = m_centroid[1] - m_dim;
406  break;
407  case eLeft:
408  expectedfx = m_centroid[2] + m_dim;
409  break;
410  case eRight:
411  expectedfx = m_centroid[2] - m_dim;
412  break;
413  case eForward:
414  expectedfx = m_centroid[0] + m_dim;
415  break;
416  case eBack:
417  expectedfx = m_centroid[0] - m_dim;
418  break;
419  }
420  if (fabs(m_octants[i]->FX(it->first) - expectedfx) > 1E-6)
421  {
422  valid = false;
423  cout << "wall neigbour error" << endl;
424  cout << expectedfx << " " << m_octants[i]->FX(it->first)
425  << " " << it->first << endl;
426  }
427  }
428  else if (it->second.size() == 1)
429  {
430  if (!(m_octants[i]->DX() == it->second[0]->DX() ||
431  it->second[0]->DX() == 2.0 * m_octants[i]->DX()))
432  {
433  valid = false;
434  cout << " 1 neigbour error" << endl;
435  cout << m_octants[i]->DX() << " " << it->second[0]->DX()
436  << endl;
437  }
438  }
439  else if (it->second.size() == 4)
440  {
441  if (!(m_octants[i]->DX() / 2.0 == it->second[0]->DX()))
442  {
443  valid = false;
444  cout << "4 neibour error" << endl;
445  cout << m_octants[i]->DX() << " " << it->second[0]->DX()
446  << endl;
447  }
448  }
449  }
450  if (!valid)
451  {
452  error = true;
453  cout << "invalid neigbour config" << endl;
454  }
455  }
456  return !error;
457 }
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:230
double NekDouble
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
NekDouble m_dim
physical size of the octree
Definition: Octree.h:226
Array< OneD, NekDouble > m_centroid
x,y,z location of the center of the octree
Definition: Octree.h:224
void Nektar::NekMeshUtils::Octree::WriteOctree ( std::string  nm)

populates the mesh m with a invalid hexahedral mesh based on the octree, used for visualisation

Parameters
nmname of the mesh file to be made

Definition at line 202 of file Octree.cpp.

References Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::NekMeshUtils::eBack, Nektar::NekMeshUtils::eDown, Nektar::NekMeshUtils::eForward, Nektar::LibUtilities::eHexahedron, Nektar::NekMeshUtils::eLeft, Nektar::FieldUtils::eOutputModule, Nektar::NekMeshUtils::eRight, Nektar::NekMeshUtils::eUp, Nektar::NekMeshUtils::GetElementFactory(), Nektar::FieldUtils::GetModuleFactory(), CG_Iterations::Mesh, and class_topology::Node.

203 {
204  MeshSharedPtr oct = boost::shared_ptr<Mesh>(new Mesh());
205  oct->m_expDim = 3;
206  oct->m_spaceDim = 3;
207  oct->m_nummode = 2;
208 
209  for (int i = 0; i < m_octants.size(); i++)
210  {
211  /*if(m_octants[i]->GetLocation() != eOnBoundary)
212  {
213  continue;
214  }*/
215 
216  vector<NodeSharedPtr> ns(8);
217 
218  ns[0] = boost::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eBack),
219  m_octants[i]->FX(eDown),
220  m_octants[i]->FX(eRight)));
221 
222  ns[1] = boost::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eForward),
223  m_octants[i]->FX(eDown),
224  m_octants[i]->FX(eRight)));
225 
226  ns[2] = boost::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eForward),
227  m_octants[i]->FX(eUp),
228  m_octants[i]->FX(eRight)));
229 
230  ns[3] = boost::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eBack),
231  m_octants[i]->FX(eUp),
232  m_octants[i]->FX(eRight)));
233 
234  ns[4] = boost::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eBack),
235  m_octants[i]->FX(eDown),
236  m_octants[i]->FX(eLeft)));
237 
238  ns[5] = boost::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eForward),
239  m_octants[i]->FX(eDown),
240  m_octants[i]->FX(eLeft)));
241 
242  ns[6] = boost::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eForward),
243  m_octants[i]->FX(eUp),
244  m_octants[i]->FX(eLeft)));
245 
246  ns[7] = boost::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eBack),
247  m_octants[i]->FX(eUp),
248  m_octants[i]->FX(eLeft)));
249 
250  vector<int> tags;
251  tags.push_back(0);
252  ElmtConfig conf(LibUtilities::eHexahedron, 1, false, false);
254  LibUtilities::eHexahedron, conf, ns, tags);
255  oct->m_element[3].push_back(E);
256  }
257 
258  ModuleSharedPtr mod =
260  mod->RegisterConfig("outfile", nm);
261  mod->ProcessVertices();
262  mod->ProcessEdges();
263  mod->ProcessFaces();
264  mod->ProcessElements();
265  mod->ProcessComposites();
266  mod->Process();
267 }
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< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:230
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
boost::shared_ptr< Module > ModuleSharedPtr
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
std::pair< ModuleType, std::string > ModuleKey

Member Data Documentation

Array<OneD, NekDouble> Nektar::NekMeshUtils::Octree::m_centroid
private

x,y,z location of the center of the octree

Definition at line 224 of file Octree.h.

NekDouble Nektar::NekMeshUtils::Octree::m_dim
private

physical size of the octree

Definition at line 226 of file Octree.h.

NekDouble Nektar::NekMeshUtils::Octree::m_eps
private

curavture sensivity paramter

Definition at line 222 of file Octree.h.

Referenced by SetParameters().

std::vector<linesource> Nektar::NekMeshUtils::Octree::m_lsources
private

Definition at line 239 of file Octree.h.

OctantSharedPtr Nektar::NekMeshUtils::Octree::m_masteroct
private

master octant for searching

Definition at line 232 of file Octree.h.

NekDouble Nektar::NekMeshUtils::Octree::m_maxDelta
private

maximum delta in the octree

Definition at line 220 of file Octree.h.

Referenced by SetParameters().

MeshSharedPtr Nektar::NekMeshUtils::Octree::m_mesh
private

Mesh object.

Definition at line 236 of file Octree.h.

NekDouble Nektar::NekMeshUtils::Octree::m_minDelta
private

minimum delta in the octree

Definition at line 218 of file Octree.h.

Referenced by SetParameters().

int Nektar::NekMeshUtils::Octree::m_numoct
private

number of octants made, used for id index

Definition at line 234 of file Octree.h.

std::vector<OctantSharedPtr> Nektar::NekMeshUtils::Octree::m_octants
private

list of leaf octants

Definition at line 230 of file Octree.h.

std::string Nektar::NekMeshUtils::Octree::m_refinement
private

Definition at line 238 of file Octree.h.

Referenced by Refinement().

std::vector<SPBaseSharedPtr> Nektar::NekMeshUtils::Octree::m_SPList
private

list of source points

Definition at line 228 of file Octree.h.