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::Octree Class Reference

class for octree More...

#include <Octree.h>

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

Public Member Functions

 Octree (CADSystemSharedPtr cad, const bool ver, const NekDouble min, const NekDouble max, const NekDouble eps, const string uds)
 Defualt constructor. More...
 
void Build ()
 executes octree building routines 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 GetOctreeMesh (MeshSharedPtr m)
 populates the mesh m with a invalid hexahedral mesh based on the octree, used for visualisation More...
 

Private Member Functions

void SmoothAllOctants ()
 Smooths specification over all octants to a gradation criteria. More...
 
void CompileCuravturePointList ()
 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...
 
CADSystemSharedPtr m_cad
 cad object More...
 
bool m_verbose
 verbose output 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
< CurvaturePointSharedPtr
m_cpList
 list of curvature sample 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...
 
string m_udsfile
 

Friends

class MemoryManager< Octree >
 

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 60 of file Octree.h.

Constructor & Destructor Documentation

Nektar::NekMeshUtils::Octree::Octree ( CADSystemSharedPtr  cad,
const bool  ver,
const NekDouble  min,
const NekDouble  max,
const NekDouble  eps,
const string  uds 
)
inline

Defualt constructor.

Parameters
cadCAD object
verbool verbose

Definition at line 71 of file Octree.h.

77  : m_minDelta(min), m_maxDelta(max), m_eps(eps), m_cad(cad),
78  m_verbose(ver), m_udsfile(uds)
79  {
80  }
CADSystemSharedPtr m_cad
cad object
Definition: Octree.h:167
NekDouble m_minDelta
minimum delta in the octree
Definition: Octree.h:161
bool m_verbose
verbose output
Definition: Octree.h:169
NekDouble m_eps
curavture sensivity paramter
Definition: Octree.h:165
NekDouble m_maxDelta
maximum delta in the octree
Definition: Octree.h:163

Member Function Documentation

void Nektar::NekMeshUtils::Octree::Build ( )

executes octree building routines

Definition at line 185 of file Octree.cpp.

186 {
187  Array<OneD, NekDouble> boundingBox = m_cad->GetBoundingBox();
188 
189  if (m_verbose)
190  cout << endl << "Octree system" << endl;
191 
192  // build curvature samples
194 
195  if (m_verbose)
196  cout << "\tCurvature samples: " << m_cpList.size() << endl;
197 
198  m_dim = max((boundingBox[1] - boundingBox[0]) / 2.0,
199  (boundingBox[3] - boundingBox[2]) / 2.0);
200 
201  m_dim = max(m_dim, (boundingBox[5] - boundingBox[4]) / 2.0);
202 
203  m_centroid = Array<OneD, NekDouble>(3);
204  m_centroid[0] = (boundingBox[1] + boundingBox[0]) / 2.0;
205  m_centroid[1] = (boundingBox[3] + boundingBox[2]) / 2.0;
206  m_centroid[2] = (boundingBox[5] + boundingBox[4]) / 2.0;
207 
208  // make master octant based on the bounding box of the domain
210  0, m_centroid[0], m_centroid[1], m_centroid[2], m_dim, m_cpList);
211 
212  SubDivide();
213 
214  m_octants.clear();
215  m_masteroct->CompileLeaves(m_octants);
216 
217  if (m_verbose)
218  {
219  cout << "\tOctants: " << m_octants.size() << endl;
220  }
221 
223 
224  PropagateDomain();
225 
227 
228  if (m_verbose)
229  {
230  int elem = CountElemt();
231  cout << "\tPredicted mesh: " << elem << endl;
232  }
233 }
CADSystemSharedPtr m_cad
cad object
Definition: Octree.h:167
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:235
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:177
bool m_verbose
verbose output
Definition: Octree.h:169
OctantSharedPtr m_masteroct
master octant for searching
Definition: Octree.h:179
void CompileCuravturePointList()
gets an optimum number of curvature sampling points and calculates the curavture at these points ...
Definition: Octree.cpp:845
void PropagateDomain()
takes the mesh specification from surface octants and progates that through the domain so all octants...
Definition: Octree.cpp:477
int CountElemt()
estimates the number of elements to be created in the mesh
Definition: Octree.cpp:681
NekDouble m_dim
physical size of the octree
Definition: Octree.h:173
Array< OneD, NekDouble > m_centroid
x,y,z location of the center of the octree
Definition: Octree.h:171
std::vector< CurvaturePointSharedPtr > m_cpList
list of curvature sample points
Definition: Octree.h:175
void SmoothSurfaceOctants()
Smooths specification over the surface encompasing octants to a gradation criteria.
Definition: Octree.cpp:420
void SmoothAllOctants()
Smooths specification over all octants to a gradation criteria.
Definition: Octree.cpp:632
void Nektar::NekMeshUtils::Octree::CompileCuravturePointList ( )
private

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

need to add curvature points with the false tag to make octree modification work

Definition at line 845 of file Octree.cpp.

References Nektar::NekMeshUtils::surf.

846 {
847  for (int i = 1; i <= m_cad->GetNumSurf(); i++)
848  {
849  CADSurfSharedPtr surf = m_cad->GetSurf(i);
850  Array<OneD, NekDouble> bounds = surf->GetBounds();
851 
852  // to figure out the amount of curvature sampling to be conducted on
853  // each parameter plane the surface is first sampled with a 40x40 grid
854  // the real space lengths of this grid are analysed to find the largest
855  // strecthing in the u and v directions
856  // this stretching this then cosnidered with the mindelta user input
857  // to find a number of sampling points in each direction which
858  // enures that in the final octree each surface octant will have at
859  // least
860  // one sample point within its volume.
861  // the 40x40 grid is used to ensure each surface has a minimum of 40x40
862  // samples.
863  NekDouble du = (bounds[1] - bounds[0]) / (40 - 1);
864  NekDouble dv = (bounds[3] - bounds[2]) / (40 - 1);
865 
866  NekDouble DeltaU = 0.0;
867  NekDouble DeltaV = 0.0;
868 
869  Array<TwoD, Array<OneD, NekDouble> > samplepoints(40, 40);
870 
871  for (int j = 0; j < 40; j++)
872  {
873  for (int k = 0; k < 40; k++)
874  {
875  Array<OneD, NekDouble> uv(2);
876  uv[0] = k * du + bounds[0];
877  uv[1] = j * dv + bounds[2];
878  if (j == 40 - 1)
879  uv[1] = bounds[3];
880  if (k == 40 - 1)
881  uv[0] = bounds[1];
882  samplepoints[k][j] = surf->P(uv);
883  }
884  }
885 
886  for (int j = 0; j < 40 - 1; j++)
887  {
888  for (int k = 0; k < 40 - 1; k++)
889  {
890  NekDouble deltau = sqrt(
891  (samplepoints[k][j][0] - samplepoints[k + 1][j][0]) *
892  (samplepoints[k][j][0] - samplepoints[k + 1][j][0]) +
893  (samplepoints[k][j][1] - samplepoints[k + 1][j][1]) *
894  (samplepoints[k][j][1] - samplepoints[k + 1][j][1]) +
895  (samplepoints[k][j][2] - samplepoints[k + 1][j][2]) *
896  (samplepoints[k][j][2] - samplepoints[k + 1][j][2]));
897  NekDouble deltav = sqrt(
898  (samplepoints[k][j][0] - samplepoints[k][j + 1][0]) *
899  (samplepoints[k][j][0] - samplepoints[k][j + 1][0]) +
900  (samplepoints[k][j][1] - samplepoints[k][j + 1][1]) *
901  (samplepoints[k][j][1] - samplepoints[k][j + 1][1]) +
902  (samplepoints[k][j][2] - samplepoints[k][j + 1][2]) *
903  (samplepoints[k][j][2] - samplepoints[k][j + 1][2]));
904 
905  if (deltau > DeltaU)
906  DeltaU = deltau;
907  if (deltav > DeltaV)
908  DeltaV = deltav;
909  }
910  }
911 
912  // these are the acutal number of sample points in each parametric
913  // direction
914  int nu = ceil(DeltaU / m_minDelta) * 40;
915  int nv = ceil(DeltaV / m_minDelta) * 40;
916 
917  for (int j = 0; j < nu; j++)
918  {
919  for (int k = 0; k < nv; k++)
920  {
921  Array<OneD, NekDouble> uv(2);
922  uv[0] = (bounds[1] - bounds[0]) / (nu - 1) * j + bounds[0];
923  uv[1] = (bounds[3] - bounds[2]) / (nv - 1) * k + bounds[2];
924 
925  // this prevents round off error at the end of the surface
926  // may not be neseercary but works
927  if (j == nu - 1)
928  uv[0] = bounds[1];
929  if (k == nv - 1)
930  uv[1] = bounds[3];
931 
932  NekDouble C = surf->Curvature(uv);
933 
934  // create new point based on smallest R, flat surfaces have k=0
935  // but still need a point for element estimation
936  if (C != 0.0)
937  {
938  bool minlimited = false;
939  NekDouble ideal;
940 
941  NekDouble del =
942  2.0 * (1.0 / C) * sqrt(m_eps * (2.0 - m_eps));
943 
944  if (del > m_maxDelta)
945  {
946  del = m_maxDelta;
947  }
948  if (del < m_minDelta)
949  {
950  ideal = del;
951  del = m_minDelta;
952  minlimited = true;
953  }
954 
955  if (minlimited)
956  {
957  CurvaturePointSharedPtr newCPoint =
959  surf->GetId(), uv, surf->P(uv), del, ideal);
960 
961  m_cpList.push_back(newCPoint);
962  }
963  else
964  {
965  CurvaturePointSharedPtr newCPoint =
967  surf->GetId(), uv, surf->P(uv), del);
968 
969  m_cpList.push_back(newCPoint);
970  }
971  }
972  else
973  {
974  CurvaturePointSharedPtr newCPoint =
976  surf->GetId(), uv, surf->P(uv));
977 
978  m_cpList.push_back(newCPoint);
979  }
980  }
981  }
982  }
983 
984  if (m_udsfile == "N")
985  {
986  return;
987  }
988 
989  // now deal with the user defined spacing
990  vector<linesource> lsources;
991  fstream fle;
992  fle.open(m_udsfile.c_str());
993 
994  string fileline;
995 
996  while (!fle.eof())
997  {
998  getline(fle, fileline);
999  stringstream s(fileline);
1000  string word;
1001  s >> word;
1002  if (word == "#")
1003  {
1004  continue;
1005  }
1006 
1007  Array<OneD, NekDouble> x1(3), x2(3);
1008  NekDouble r, d;
1009  x1[0] = boost::lexical_cast<double>(word);
1010  s >> x1[1] >> x1[2] >> x2[0] >> x2[1] >> x2[2] >> r >> d;
1011 
1012  lsources.push_back(linesource(x1, x2, r, d));
1013  }
1014  fle.close();
1015 
1016  for (int j = 0; j < lsources.size(); j++)
1017  {
1018  cout << lsources[j].x1[0] << " " << lsources[j].x1[1] << " "
1019  << lsources[j].x1[2] << endl;
1020  cout << lsources[j].x2[0] << " " << lsources[j].x2[1] << " "
1021  << lsources[j].x2[2] << endl;
1022  cout << lsources[j].Length() << endl;
1023  }
1024 
1025  int ct = 0;
1026  for (int i = 0; i < m_cpList.size(); i++)
1027  {
1028  for (int j = 0; j < lsources.size(); j++)
1029  {
1030  if (lsources[j].withinRange(m_cpList[i]->GetLoc()))
1031  {
1032  ct++;
1033  m_cpList[i]->SetDelta(lsources[j].delta);
1034  }
1035  }
1036  }
1037  cout << ct << endl;
1038 
1039  ///@TODO need to add curvature points with the false tag to make octree
1040  ///modification work
1041  // off surfaces
1042 
1043  /*for(int i = 0; i < lsources.size(); i++)
1044  {
1045  int nc; //number of point to add cicularly
1046  int nl; //number of point to add length
1047 
1048  nc = ceil(2.0*3.142*lsources[i].R / lsources[i].delta)*2;
1049  nl = ceil(lsources[i].Length() / lsources[i].delta)*2;
1050 
1051  NekDouble dr = lsources[i].Length() / nl;
1052  NekDouble dtheta = 2.0*3.142 / nc;
1053 
1054  for(int j = 0; j < nl; j++)
1055  {
1056  NekDouble len =
1057  for(int k = 0; k < nc; k++)
1058  {
1059 
1060  }
1061  }
1062 
1063  }*/
1064 }
CADSystemSharedPtr m_cad
cad object
Definition: Octree.h:167
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:161
boost::shared_ptr< CurvaturePoint > CurvaturePointSharedPtr
double NekDouble
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:217
NekDouble m_eps
curavture sensivity paramter
Definition: Octree.h:165
std::vector< CurvaturePointSharedPtr > m_cpList
list of curvature sample points
Definition: Octree.h:175
NekDouble m_maxDelta
maximum delta in the octree
Definition: Octree.h:163
int Nektar::NekMeshUtils::Octree::CountElemt ( )
private

estimates the number of elements to be created in the mesh

Definition at line 681 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.

682 {
683  // by considering the volume of a tet evaluate the number of elements in the
684  // mesh
685 
686  NekDouble total = 0.0;
687 
688  Array<OneD, NekDouble> boundingBox = m_cad->GetBoundingBox();
689 
690  for (int i = 0; i < m_octants.size(); i++)
691  {
692  OctantSharedPtr oct = m_octants[i];
693  if (oct->GetLocation() == eInside)
694  {
695  total += 8.0 * oct->DX() * oct->DX() * oct->DX() /
696  (oct->GetDelta() * oct->GetDelta() * oct->GetDelta() /
697  6.0 / sqrt(2));
698  }
699  else if (oct->GetLocation() == eOnBoundary)
700  {
701  NekDouble vol = 1.0;
702  if (oct->FX(eBack) < boundingBox[1] &&
703  oct->FX(eForward) > boundingBox[0])
704  {
705  // then there is some over lap in x
706  NekDouble min, max;
707  if (oct->FX(eBack) > boundingBox[0])
708  {
709  min = oct->FX(eBack);
710  }
711  else
712  {
713  min = boundingBox[0];
714  }
715  if (boundingBox[1] < oct->FX(eForward))
716  {
717  max = boundingBox[1];
718  }
719  else
720  {
721  max = oct->FX(eForward);
722  }
723  vol *= (max - min);
724  }
725  else
726  {
727  vol *= 0.0;
728  }
729 
730  if (oct->FX(eDown) < boundingBox[3] &&
731  oct->FX(eUp) > boundingBox[2])
732  {
733  // then there is some over lap in x
734  NekDouble min, max;
735  if (oct->FX(eDown) > boundingBox[2])
736  {
737  min = oct->FX(eDown);
738  }
739  else
740  {
741  min = boundingBox[2];
742  }
743  if (boundingBox[3] < oct->FX(eUp))
744  {
745  max = boundingBox[3];
746  }
747  else
748  {
749  max = oct->FX(eUp);
750  }
751  vol *= (max - min);
752  }
753  else
754  {
755  vol *= 0.0;
756  }
757 
758  if (oct->FX(eRight) < boundingBox[5] &&
759  oct->FX(eLeft) > boundingBox[4])
760  {
761  // then there is some over lap in x
762  NekDouble min, max;
763  if (oct->FX(eRight) > boundingBox[4])
764  {
765  min = oct->FX(eRight);
766  }
767  else
768  {
769  min = boundingBox[4];
770  }
771  if (boundingBox[5] < oct->FX(eLeft))
772  {
773  max = boundingBox[5];
774  }
775  else
776  {
777  max = oct->FX(eLeft);
778  }
779  vol *= (max - min);
780  }
781  else
782  {
783  vol *= 0.0;
784  }
785  total += vol / 2.0 / (oct->GetDelta() * oct->GetDelta() *
786  oct->GetDelta() / 6.0 / sqrt(2));
787  }
788  }
789 
790  return int(total);
791 }
CADSystemSharedPtr m_cad
cad object
Definition: Octree.h:167
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:177
double NekDouble
boost::shared_ptr< Octant > OctantSharedPtr
Definition: Octant.h:74
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 1066 of file Octree.cpp.

1067 {
1068  return fabs(i->GetDelta() - j->GetDelta()) / i->Distance(j);
1069 }
NekDouble Nektar::NekMeshUtils::Octree::GetMinDelta ( )
inline

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

Returns
miminum delta in octree

Definition at line 101 of file Octree.h.

References m_minDelta.

102  {
103  return m_minDelta;
104  }
NekDouble m_minDelta
minimum delta in the octree
Definition: Octree.h:161
void Nektar::NekMeshUtils::Octree::GetOctreeMesh ( MeshSharedPtr  m)

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

Definition at line 125 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::NekMeshUtils::eRight, Nektar::NekMeshUtils::eUp, and Nektar::NekMeshUtils::GetElementFactory().

126 {
127  for (int i = 0; i < m_octants.size(); i++)
128  {
129  /*if(m_octants[i]->GetLocation() != eOnBoundary)
130  {
131  continue;
132  }*/
133 
134  vector<NodeSharedPtr> ns(8);
135 
136  ns[0] = boost::shared_ptr<Node>(new Node(0,
137  m_octants[i]->FX(eBack),
138  m_octants[i]->FX(eDown),
139  m_octants[i]->FX(eRight)));
140 
141  ns[1] = boost::shared_ptr<Node>(new Node(0,
142  m_octants[i]->FX(eForward),
143  m_octants[i]->FX(eDown),
144  m_octants[i]->FX(eRight)));
145 
146  ns[2] = boost::shared_ptr<Node>(new Node(0,
147  m_octants[i]->FX(eForward),
148  m_octants[i]->FX(eUp),
149  m_octants[i]->FX(eRight)));
150 
151  ns[3] = boost::shared_ptr<Node>(new Node(0,
152  m_octants[i]->FX(eBack),
153  m_octants[i]->FX(eUp),
154  m_octants[i]->FX(eRight)));
155 
156  ns[4] = boost::shared_ptr<Node>(new Node(0,
157  m_octants[i]->FX(eBack),
158  m_octants[i]->FX(eDown),
159  m_octants[i]->FX(eLeft)));
160 
161  ns[5] = boost::shared_ptr<Node>(new Node(0,
162  m_octants[i]->FX(eForward),
163  m_octants[i]->FX(eDown),
164  m_octants[i]->FX(eLeft)));
165 
166  ns[6] = boost::shared_ptr<Node>(new Node(0,
167  m_octants[i]->FX(eForward),
168  m_octants[i]->FX(eUp),
169  m_octants[i]->FX(eLeft)));
170 
171  ns[7] = boost::shared_ptr<Node>(new Node(0,
172  m_octants[i]->FX(eBack),
173  m_octants[i]->FX(eUp),
174  m_octants[i]->FX(eLeft)));
175 
176  vector<int> tags;
177  tags.push_back(0);
178  ElmtConfig conf(LibUtilities::eHexahedron, 1, false, false);
180  LibUtilities::eHexahedron, conf, ns, tags);
181  m->m_element[3].push_back(E);
182  }
183 }
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:177
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
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 477 of file Octree.cpp.

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

478 {
479  // until all m_octants know their delta specifcation and orientaion
480  // look over all m_octants and if their neighours know either their
481  // orientation
482  // or specifcation calculate one for this octant
483  int ct = 0;
484 
485  do
486  {
487  ct = 0;
488  for (int i = 0; i < m_octants.size(); i++)
489  {
490  OctantSharedPtr oct = m_octants[i];
491 
492  if (!oct->IsDeltaKnown())
493  { // if delta has not been asigned
494  vector<OctantSharedPtr> known;
495  map<OctantFace, vector<OctantSharedPtr> > nList =
496  oct->GetNeigbours();
497  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
498 
499  for (it = nList.begin(); it != nList.end(); it++)
500  {
501  for (int j = 0; j < it->second.size(); j++)
502  {
503  if (it->second[j]->IsDeltaKnown())
504  {
505  known.push_back(it->second[j]);
506  }
507  }
508  }
509 
510  if (known.size() > 0)
511  {
512  vector<NekDouble> deltaPrime;
513  for (int j = 0; j < known.size(); j++)
514  {
515  NekDouble r = oct->Distance(known[j]);
516 
517  if (0.14 * r + known[j]->GetDelta() < m_maxDelta)
518  {
519  deltaPrime.push_back(0.14 * r +
520  known[j]->GetDelta());
521  }
522  else
523  {
524  deltaPrime.push_back(m_maxDelta);
525  }
526  }
527  NekDouble min = numeric_limits<double>::max();
528  for (int j = 0; j < deltaPrime.size(); j++)
529  {
530  if (deltaPrime[j] < min)
531  {
532  min = deltaPrime[j];
533  }
534  }
535  oct->SetDelta(min);
536  ct += 1;
537  deltaPrime.clear();
538  }
539  known.clear();
540  }
541 
542  if (oct->GetLocation() == eUnknown)
543  { // if the node does not know its location
544  vector<OctantSharedPtr> known;
545  map<OctantFace, vector<OctantSharedPtr> > nList =
546  oct->GetNeigbours();
547  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
548 
549  for (it = nList.begin(); it != nList.end(); it++)
550  {
551  for (int j = 0; j < it->second.size(); j++)
552  {
553  if (it->second[j]->GetLocation() != eUnknown)
554  {
555  known.push_back(it->second[j]);
556  }
557  }
558  }
559 
560  if (known.size() > 0)
561  {
562  vector<OctantSharedPtr> isNotOnBound;
563  for (int j = 0; j < known.size(); j++)
564  {
565  if (known[j]->GetLocation() != eOnBoundary)
566  {
567  isNotOnBound.push_back(known[j]);
568  }
569  }
570 
571  if (isNotOnBound.size() > 0)
572  {
573  oct->SetLocation(isNotOnBound[0]->GetLocation());
574  }
575  else
576  {
577  NekDouble dist = numeric_limits<double>::max();
578 
579  OctantSharedPtr closest;
580 
581  for (int j = 0; j < known.size(); j++)
582  {
583  if (oct->Distance(known[j]) < dist)
584  {
585  closest = known[j];
586  dist = oct->Distance(known[j]);
587  }
588  }
589 
590  CurvaturePointSharedPtr cp = closest->GetCPPoint();
591 
592  Array<OneD, NekDouble> octloc, cploc, vec(3), uv, N;
593  int surf;
594  cp->GetCAD(surf, uv);
595  N = m_cad->GetSurf(surf)->N(uv);
596 
597  octloc = oct->GetLoc();
598  cploc = cp->GetLoc();
599 
600  vec[0] = octloc[0] - cploc[0];
601  vec[1] = octloc[1] - cploc[1];
602  vec[2] = octloc[2] - cploc[2];
603 
604  NekDouble dot =
605  vec[0] * N[0] + vec[1] * N[1] + vec[2] * N[2];
606 
607  if (dot <= 0.0)
608  {
609  oct->SetLocation(eOutside);
610  ct += 1;
611  }
612  else
613  {
614  oct->SetLocation(eInside);
615  ct += 1;
616  }
617  }
618  }
619  known.clear();
620  }
621  }
622 
623  } while (ct > 0);
624 
625  for (int i = 0; i < m_octants.size(); i++)
626  {
627  ASSERTL0(m_octants[i]->IsDeltaKnown(),
628  "does not know delta after propergation");
629  }
630 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
CADSystemSharedPtr m_cad
cad object
Definition: Octree.h:167
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:177
boost::shared_ptr< CurvaturePoint > CurvaturePointSharedPtr
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
NekDouble m_maxDelta
maximum delta in the octree
Definition: Octree.h:163
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 48 of file Octree.cpp.

References ASSERTL0.

49 {
50  // starting at master octant 0 move through succsesive m_octants which
51  // contain
52  // the point loc until a leaf is found
54  int quad;
55 
56  bool found = false;
57 
58  while (!found)
59  {
60  Array<OneD, NekDouble> octloc = n->GetLoc();
61 
62  if (!(loc[0] < octloc[0]) && // forward
63  !(loc[1] < octloc[1]) && // forward
64  !(loc[2] < octloc[2])) // forward
65  {
66  quad = 0;
67  }
68  else if (!(loc[0] < octloc[0]) && // forward
69  !(loc[1] < octloc[1]) && // forward
70  !(loc[2] > octloc[2])) // back
71  {
72  quad = 1;
73  }
74  else if (!(loc[0] < octloc[0]) && // forward
75  !(loc[1] > octloc[1]) && // back
76  !(loc[2] < octloc[2])) // forward
77  {
78  quad = 2;
79  }
80  else if (!(loc[0] < octloc[0]) && // forward
81  !(loc[1] > octloc[1]) && // back
82  !(loc[2] > octloc[2])) // back
83  {
84  quad = 3;
85  }
86  else if (!(loc[0] > octloc[0]) && // back
87  !(loc[1] < octloc[1]) && // forward
88  !(loc[2] < octloc[2])) // forward
89  {
90  quad = 4;
91  }
92  else if (!(loc[0] > octloc[0]) && // back
93  !(loc[1] < octloc[1]) && // forward
94  !(loc[2] > octloc[2])) // back
95  {
96  quad = 5;
97  }
98  else if (!(loc[0] > octloc[0]) && // back
99  !(loc[1] > octloc[1]) && // back
100  !(loc[2] < octloc[2])) // forward
101  {
102  quad = 6;
103  }
104  else if (!(loc[0] > octloc[0]) && // back
105  !(loc[1] > octloc[1]) && // back
106  !(loc[2] > octloc[2])) // back
107  {
108  quad = 7;
109  }
110  else
111  {
112  ASSERTL0(false, "Cannot locate quadrant");
113  }
114 
115  n = n->GetChild(quad);
116 
117  if (n->IsLeaf())
118  {
119  found = true;
120  }
121  }
122  return n->GetDelta();
123 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
OctantSharedPtr m_masteroct
master octant for searching
Definition: Octree.h:179
boost::shared_ptr< Octant > OctantSharedPtr
Definition: Octant.h:74
void Nektar::NekMeshUtils::Octree::SmoothAllOctants ( )
private

Smooths specification over all octants to a gradation criteria.

Definition at line 632 of file Octree.cpp.

References Nektar::iterator.

633 {
634  // until no more changes occur smooth the mesh specification between all
635  // m_octants not particualrly strictly
636  int ct = 0;
637 
638  do
639  {
640  ct = 0;
641  for (int i = 0; i < m_octants.size(); i++)
642  {
643  OctantSharedPtr oct = m_octants[i];
644 
645  vector<OctantSharedPtr> check;
646  map<OctantFace, vector<OctantSharedPtr> > nList =
647  oct->GetNeigbours();
648  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
649  for (it = nList.begin(); it != nList.end(); it++)
650  {
651  for (int j = 0; j < it->second.size(); j++)
652  {
653  if (it->second[j]->GetDelta() < oct->GetDelta() &&
654  ddx(oct, it->second[j]) > 0.2)
655  {
656  check.push_back(it->second[j]);
657  }
658  }
659  }
660 
661  if (check.size() > 0)
662  {
663  NekDouble deltaSM = numeric_limits<double>::max();
664  for (int j = 0; j < check.size(); j++)
665  {
666  NekDouble r = oct->Distance(check[j]);
667 
668  if (0.199 * r + check[j]->GetDelta() < deltaSM)
669  {
670  deltaSM = 0.199 * r + check[j]->GetDelta();
671  }
672  }
673  oct->SetDelta(deltaSM);
674  ct += 1;
675  }
676  }
677 
678  } while (ct > 0);
679 }
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:1066
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:177
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 420 of file Octree.cpp.

References Nektar::iterator.

421 {
422  // for all the m_octants which are surface containing and know their delta
423  // specification, look over all neighbours and ensure the specification
424  // between them is smooth
425  int ct = 0;
426 
427  do
428  {
429  ct = 0;
430  for (int i = 0; i < m_octants.size(); i++)
431  {
432  OctantSharedPtr oct = m_octants[i];
433 
434  if (oct->IsDeltaKnown())
435  {
436  vector<OctantSharedPtr> check;
437  map<OctantFace, vector<OctantSharedPtr> > nList =
438  oct->GetNeigbours();
439  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
440 
441  for (it = nList.begin(); it != nList.end(); it++)
442  {
443  for (int j = 0; j < it->second.size(); j++)
444  {
445  if (it->second[j]->IsDeltaKnown() &&
446  it->second[j]->GetDelta() < oct->GetDelta() &&
447  ddx(oct, it->second[j]) > 0.1)
448  {
449  check.push_back(it->second[j]);
450  }
451  }
452  }
453 
454  // for each neighbour listed in check_id, figure out the
455  // smoothed
456  // delta, and asign the miminum of these to nodes[i].GetDelta()
457  if (check.size() > 0)
458  {
459  NekDouble deltaSM = numeric_limits<double>::max();
460  for (int j = 0; j < check.size(); j++)
461  {
462  NekDouble r = oct->Distance(check[j]);
463 
464  if (0.099 * r + check[j]->GetDelta() < deltaSM)
465  {
466  deltaSM = 0.099 * r + check[j]->GetDelta();
467  }
468  }
469  oct->SetDelta(deltaSM);
470  ct += 1;
471  }
472  }
473  }
474  } while (ct > 0);
475 }
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:1066
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:177
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 235 of file Octree.cpp.

References Nektar::iterator.

236 {
237  bool repeat;
238  int ct = 1;
239 
240  m_numoct = 1;
241  m_masteroct->Subdivide(m_masteroct, m_numoct);
242 
243  if (m_verbose)
244  {
245  cout << "\tSubdivide iteration: ";
246  }
247 
248  do
249  {
250  ct++;
251  if (m_verbose)
252  {
253  cout << ct << " ";
254  cout.flush();
255  }
256  repeat = false;
257  m_octants.clear();
258  // grab a list of the leaves curently in the octree
259  m_masteroct->CompileLeaves(m_octants);
260 
261  VerifyNeigbours();
262 
263  // neeed to create a divide list, in the first list will be m_octants
264  // which need to
265  // subdivide based on curvature,
266  // in the next list will be ocants which need to subdivide to make sure
267  // level criteria is statisified for the previous list and so on
268  // the list will keep building till no more lists are required.
269  // the list will then be iterated through backwards to subdivide all the
270  // sublists in turn.
271  vector<vector<OctantSharedPtr> > dividelist;
272  set<int> inlist;
273  // build initial list
274  {
275  vector<OctantSharedPtr> sublist;
276  for (int i = 0; i < m_octants.size(); i++)
277  {
278  if (m_octants[i]->NeedDivide() &&
279  m_octants[i]->DX() / 4.0 > m_minDelta)
280  {
281  sublist.push_back(m_octants[i]);
282  inlist.insert(m_octants[i]->GetId());
283  repeat = true; // if there is a subdivision this whole
284  // process needs to be repeated
285  }
286  }
287  dividelist.push_back(sublist);
288  }
289  // then loop over building sublists until no more are required
290  int ct2 = 0;
291  while (true)
292  {
293  ct2++;
294  vector<OctantSharedPtr> newsublist,
295  previouslist = dividelist.back();
296  for (int i = 0; i < previouslist.size(); i++)
297  {
298  map<OctantFace, vector<OctantSharedPtr> > nlist =
299  previouslist[i]->GetNeigbours();
300  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
301  for (it = nlist.begin(); it != nlist.end(); it++)
302  {
303  for (int j = 0; j < it->second.size(); j++)
304  {
305  if (previouslist[i]->DX() < it->second[j]->DX())
306  {
308  inlist.find(it->second[j]->GetId());
309  if (s == inlist.end())
310  {
311  inlist.insert(it->second[j]->GetId());
312  newsublist.push_back(it->second[j]);
313  }
314  }
315  }
316  }
317  }
318  if (newsublist.size() == 0)
319  {
320  break;
321  }
322  else
323  {
324  dividelist.push_back(newsublist);
325  }
326  }
327 
328  vector<vector<OctantSharedPtr> >::reverse_iterator rit;
329  for (rit = dividelist.rbegin(); rit != dividelist.rend(); rit++)
330  {
331  vector<OctantSharedPtr> currentlist = *rit;
332  for (int i = 0; i < currentlist.size(); i++)
333  {
334  currentlist[i]->Subdivide(currentlist[i], m_numoct);
335  }
336  }
337  } while (repeat);
338 
339  if (m_verbose)
340  {
341  cout << endl;
342  }
343 }
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:177
NekDouble m_minDelta
minimum delta in the octree
Definition: Octree.h:161
bool m_verbose
verbose output
Definition: Octree.h:169
bool VerifyNeigbours()
Looks over all leaf octants and checks that their neigbour assigments are valid.
Definition: Octree.cpp:345
OctantSharedPtr m_masteroct
master octant for searching
Definition: Octree.h:179
int m_numoct
number of octants made, used for id index
Definition: Octree.h:181
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
bool Nektar::NekMeshUtils::Octree::VerifyNeigbours ( )
private

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

Definition at line 345 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.

346 {
347  // check all neibours
348  bool error = false;
349  for (int i = 0; i < m_octants.size(); i++)
350  {
351  bool valid = true;
352  map<OctantFace, vector<OctantSharedPtr> > nlist =
353  m_octants[i]->GetNeigbours();
354  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
355  for (it = nlist.begin(); it != nlist.end(); it++)
356  {
357  if (it->second.size() == 0)
358  {
359  NekDouble expectedfx;
360  switch (it->first)
361  {
362  case eUp:
363  expectedfx = m_centroid[1] + m_dim;
364  break;
365  case eDown:
366  expectedfx = m_centroid[1] - m_dim;
367  break;
368  case eLeft:
369  expectedfx = m_centroid[2] + m_dim;
370  break;
371  case eRight:
372  expectedfx = m_centroid[2] - m_dim;
373  break;
374  case eForward:
375  expectedfx = m_centroid[0] + m_dim;
376  break;
377  case eBack:
378  expectedfx = m_centroid[0] - m_dim;
379  break;
380  }
381  if (fabs(m_octants[i]->FX(it->first) - expectedfx) > 1E-6)
382  {
383  valid = false;
384  cout << "wall neigbour error" << endl;
385  cout << expectedfx << " " << m_octants[i]->FX(it->first)
386  << " " << it->first << endl;
387  }
388  }
389  else if (it->second.size() == 1)
390  {
391  if (!(m_octants[i]->DX() == it->second[0]->DX() ||
392  it->second[0]->DX() == 2.0 * m_octants[i]->DX()))
393  {
394  valid = false;
395  cout << " 1 neigbour error" << endl;
396  cout << m_octants[i]->DX() << " " << it->second[0]->DX()
397  << endl;
398  }
399  }
400  else if (it->second.size() == 4)
401  {
402  if (!(m_octants[i]->DX() / 2.0 == it->second[0]->DX()))
403  {
404  valid = false;
405  cout << "4 neibour error" << endl;
406  cout << m_octants[i]->DX() << " " << it->second[0]->DX()
407  << endl;
408  }
409  }
410  }
411  if (!valid)
412  {
413  error = true;
414  cout << "invalid neigbour config" << endl;
415  }
416  }
417  return !error;
418 }
std::vector< OctantSharedPtr > m_octants
list of leaf octants
Definition: Octree.h:177
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:173
Array< OneD, NekDouble > m_centroid
x,y,z location of the center of the octree
Definition: Octree.h:171

Friends And Related Function Documentation

friend class MemoryManager< Octree >
friend

Definition at line 63 of file Octree.h.

Member Data Documentation

CADSystemSharedPtr Nektar::NekMeshUtils::Octree::m_cad
private

cad object

Definition at line 167 of file Octree.h.

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

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

Definition at line 171 of file Octree.h.

std::vector<CurvaturePointSharedPtr> Nektar::NekMeshUtils::Octree::m_cpList
private

list of curvature sample points

Definition at line 175 of file Octree.h.

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

physical size of the octree

Definition at line 173 of file Octree.h.

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

curavture sensivity paramter

Definition at line 165 of file Octree.h.

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

master octant for searching

Definition at line 179 of file Octree.h.

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

maximum delta in the octree

Definition at line 163 of file Octree.h.

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

minimum delta in the octree

Definition at line 161 of file Octree.h.

Referenced by GetMinDelta().

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

number of octants made, used for id index

Definition at line 181 of file Octree.h.

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

list of leaf octants

Definition at line 177 of file Octree.h.

string Nektar::NekMeshUtils::Octree::m_udsfile
private

Definition at line 183 of file Octree.h.

bool Nektar::NekMeshUtils::Octree::m_verbose
private

verbose output

Definition at line 169 of file Octree.h.