45 namespace NekMeshUtils
62 if (!(loc[0] < octloc[0]) &&
63 !(loc[1] < octloc[1]) &&
64 !(loc[2] < octloc[2]))
68 else if (!(loc[0] < octloc[0]) &&
69 !(loc[1] < octloc[1]) &&
70 !(loc[2] > octloc[2]))
74 else if (!(loc[0] < octloc[0]) &&
75 !(loc[1] > octloc[1]) &&
76 !(loc[2] < octloc[2]))
80 else if (!(loc[0] < octloc[0]) &&
81 !(loc[1] > octloc[1]) &&
82 !(loc[2] > octloc[2]))
86 else if (!(loc[0] > octloc[0]) &&
87 !(loc[1] < octloc[1]) &&
88 !(loc[2] < octloc[2]))
92 else if (!(loc[0] > octloc[0]) &&
93 !(loc[1] < octloc[1]) &&
94 !(loc[2] > octloc[2]))
98 else if (!(loc[0] > octloc[0]) &&
99 !(loc[1] > octloc[1]) &&
100 !(loc[2] < octloc[2]))
104 else if (!(loc[0] > octloc[0]) &&
105 !(loc[1] > octloc[1]) &&
106 !(loc[2] > octloc[2]))
112 ASSERTL0(
false,
"Cannot locate quadrant");
115 n = n->GetChild(quad);
122 return n->GetDelta();
127 for (
int i = 0; i < m_octants.size(); i++)
134 vector<NodeSharedPtr> ns(8);
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)));
141 ns[1] = boost::shared_ptr<Node>(
new Node(0,
143 m_octants[i]->FX(
eDown),
144 m_octants[i]->FX(
eRight)));
146 ns[2] = boost::shared_ptr<Node>(
new Node(0,
148 m_octants[i]->FX(
eUp),
149 m_octants[i]->FX(
eRight)));
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)));
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)));
161 ns[5] = boost::shared_ptr<Node>(
new Node(0,
163 m_octants[i]->FX(
eDown),
164 m_octants[i]->FX(
eLeft)));
166 ns[6] = boost::shared_ptr<Node>(
new Node(0,
168 m_octants[i]->FX(
eUp),
169 m_octants[i]->FX(
eLeft)));
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)));
181 m->m_element[3].push_back(E);
190 cout << endl <<
"Octree system" << endl;
193 CompileCuravturePointList();
196 cout <<
"\tCurvature samples: " << m_cpList.size() << endl;
198 m_dim = max((boundingBox[1] - boundingBox[0]) / 2.0,
199 (boundingBox[3] - boundingBox[2]) / 2.0);
201 m_dim = max(m_dim, (boundingBox[5] - boundingBox[4]) / 2.0);
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;
210 0, m_centroid[0], m_centroid[1], m_centroid[2], m_dim, m_cpList);
215 m_masteroct->CompileLeaves(m_octants);
219 cout <<
"\tOctants: " << m_octants.size() << endl;
222 SmoothSurfaceOctants();
230 int elem = CountElemt();
231 cout <<
"\tPredicted mesh: " << elem << endl;
235 void Octree::SubDivide()
241 m_masteroct->Subdivide(m_masteroct, m_numoct);
245 cout <<
"\tSubdivide iteration: ";
259 m_masteroct->CompileLeaves(m_octants);
271 vector<vector<OctantSharedPtr> > dividelist;
275 vector<OctantSharedPtr> sublist;
276 for (
int i = 0; i < m_octants.size(); i++)
278 if (m_octants[i]->NeedDivide() &&
279 m_octants[i]->DX() / 4.0 > m_minDelta)
281 sublist.push_back(m_octants[i]);
282 inlist.insert(m_octants[i]->GetId());
287 dividelist.push_back(sublist);
294 vector<OctantSharedPtr> newsublist,
295 previouslist = dividelist.back();
296 for (
int i = 0; i < previouslist.size(); i++)
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++)
303 for (
int j = 0; j < it->second.size(); j++)
305 if (previouslist[i]->DX() < it->second[j]->DX())
308 inlist.find(it->second[j]->GetId());
309 if (s == inlist.end())
311 inlist.insert(it->second[j]->GetId());
312 newsublist.push_back(it->second[j]);
318 if (newsublist.size() == 0)
324 dividelist.push_back(newsublist);
328 vector<vector<OctantSharedPtr> >::reverse_iterator rit;
329 for (rit = dividelist.rbegin(); rit != dividelist.rend(); rit++)
331 vector<OctantSharedPtr> currentlist = *rit;
332 for (
int i = 0; i < currentlist.size(); i++)
334 currentlist[i]->Subdivide(currentlist[i], m_numoct);
345 bool Octree::VerifyNeigbours()
349 for (
int i = 0; i < m_octants.size(); i++)
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++)
357 if (it->second.size() == 0)
363 expectedfx = m_centroid[1] + m_dim;
366 expectedfx = m_centroid[1] - m_dim;
369 expectedfx = m_centroid[2] + m_dim;
372 expectedfx = m_centroid[2] - m_dim;
375 expectedfx = m_centroid[0] + m_dim;
378 expectedfx = m_centroid[0] - m_dim;
381 if (fabs(m_octants[i]->FX(it->first) - expectedfx) > 1E-6)
384 cout <<
"wall neigbour error" << endl;
385 cout << expectedfx <<
" " << m_octants[i]->FX(it->first)
386 <<
" " << it->first << endl;
389 else if (it->second.size() == 1)
391 if (!(m_octants[i]->DX() == it->second[0]->DX() ||
392 it->second[0]->DX() == 2.0 * m_octants[i]->DX()))
395 cout <<
" 1 neigbour error" << endl;
396 cout << m_octants[i]->DX() <<
" " << it->second[0]->DX()
400 else if (it->second.size() == 4)
402 if (!(m_octants[i]->DX() / 2.0 == it->second[0]->DX()))
405 cout <<
"4 neibour error" << endl;
406 cout << m_octants[i]->DX() <<
" " << it->second[0]->DX()
414 cout <<
"invalid neigbour config" << endl;
420 void Octree::SmoothSurfaceOctants()
430 for (
int i = 0; i < m_octants.size(); i++)
434 if (oct->IsDeltaKnown())
436 vector<OctantSharedPtr> check;
437 map<OctantFace, vector<OctantSharedPtr> > nList =
439 map<OctantFace, vector<OctantSharedPtr> >
::iterator it;
441 for (it = nList.begin(); it != nList.end(); it++)
443 for (
int j = 0; j < it->second.size(); j++)
445 if (it->second[j]->IsDeltaKnown() &&
446 it->second[j]->GetDelta() < oct->GetDelta() &&
447 ddx(oct, it->second[j]) > 0.1)
449 check.push_back(it->second[j]);
457 if (check.size() > 0)
459 NekDouble deltaSM = numeric_limits<double>::max();
460 for (
int j = 0; j < check.size(); j++)
464 if (0.099 * r + check[j]->GetDelta() < deltaSM)
466 deltaSM = 0.099 * r + check[j]->GetDelta();
469 oct->SetDelta(deltaSM);
477 void Octree::PropagateDomain()
488 for (
int i = 0; i < m_octants.size(); i++)
492 if (!oct->IsDeltaKnown())
494 vector<OctantSharedPtr> known;
495 map<OctantFace, vector<OctantSharedPtr> > nList =
497 map<OctantFace, vector<OctantSharedPtr> >
::iterator it;
499 for (it = nList.begin(); it != nList.end(); it++)
501 for (
int j = 0; j < it->second.size(); j++)
503 if (it->second[j]->IsDeltaKnown())
505 known.push_back(it->second[j]);
510 if (known.size() > 0)
512 vector<NekDouble> deltaPrime;
513 for (
int j = 0; j < known.size(); j++)
517 if (0.14 * r + known[j]->GetDelta() < m_maxDelta)
519 deltaPrime.push_back(0.14 * r +
520 known[j]->GetDelta());
524 deltaPrime.push_back(m_maxDelta);
527 NekDouble min = numeric_limits<double>::max();
528 for (
int j = 0; j < deltaPrime.size(); j++)
530 if (deltaPrime[j] < min)
544 vector<OctantSharedPtr> known;
545 map<OctantFace, vector<OctantSharedPtr> > nList =
547 map<OctantFace, vector<OctantSharedPtr> >
::iterator it;
549 for (it = nList.begin(); it != nList.end(); it++)
551 for (
int j = 0; j < it->second.size(); j++)
553 if (it->second[j]->GetLocation() !=
eUnknown)
555 known.push_back(it->second[j]);
560 if (known.size() > 0)
562 vector<OctantSharedPtr> isNotOnBound;
563 for (
int j = 0; j < known.size(); j++)
567 isNotOnBound.push_back(known[j]);
571 if (isNotOnBound.size() > 0)
573 oct->SetLocation(isNotOnBound[0]->GetLocation());
577 NekDouble dist = numeric_limits<double>::max();
581 for (
int j = 0; j < known.size(); j++)
583 if (oct->Distance(known[j]) < dist)
586 dist = oct->Distance(known[j]);
594 cp->GetCAD(surf, uv);
595 N = m_cad->GetSurf(surf)->N(uv);
597 octloc = oct->GetLoc();
598 cploc = cp->GetLoc();
600 vec[0] = octloc[0] - cploc[0];
601 vec[1] = octloc[1] - cploc[1];
602 vec[2] = octloc[2] - cploc[2];
605 vec[0] * N[0] + vec[1] * N[1] + vec[2] * N[2];
625 for (
int i = 0; i < m_octants.size(); i++)
627 ASSERTL0(m_octants[i]->IsDeltaKnown(),
628 "does not know delta after propergation");
632 void Octree::SmoothAllOctants()
641 for (
int i = 0; i < m_octants.size(); i++)
645 vector<OctantSharedPtr> check;
646 map<OctantFace, vector<OctantSharedPtr> > nList =
648 map<OctantFace, vector<OctantSharedPtr> >
::iterator it;
649 for (it = nList.begin(); it != nList.end(); it++)
651 for (
int j = 0; j < it->second.size(); j++)
653 if (it->second[j]->GetDelta() < oct->GetDelta() &&
654 ddx(oct, it->second[j]) > 0.2)
656 check.push_back(it->second[j]);
661 if (check.size() > 0)
663 NekDouble deltaSM = numeric_limits<double>::max();
664 for (
int j = 0; j < check.size(); j++)
668 if (0.199 * r + check[j]->GetDelta() < deltaSM)
670 deltaSM = 0.199 * r + check[j]->GetDelta();
673 oct->SetDelta(deltaSM);
681 int Octree::CountElemt()
690 for (
int i = 0; i < m_octants.size(); i++)
693 if (oct->GetLocation() ==
eInside)
695 total += 8.0 * oct->DX() * oct->DX() * oct->DX() /
696 (oct->GetDelta() * oct->GetDelta() * oct->GetDelta() /
702 if (oct->FX(
eBack) < boundingBox[1] &&
707 if (oct->FX(
eBack) > boundingBox[0])
709 min = oct->FX(
eBack);
713 min = boundingBox[0];
715 if (boundingBox[1] < oct->FX(
eForward))
717 max = boundingBox[1];
730 if (oct->FX(
eDown) < boundingBox[3] &&
731 oct->FX(
eUp) > boundingBox[2])
735 if (oct->FX(
eDown) > boundingBox[2])
737 min = oct->FX(
eDown);
741 min = boundingBox[2];
743 if (boundingBox[3] < oct->FX(
eUp))
745 max = boundingBox[3];
758 if (oct->FX(
eRight) < boundingBox[5] &&
759 oct->FX(
eLeft) > boundingBox[4])
763 if (oct->FX(
eRight) > boundingBox[4])
769 min = boundingBox[4];
771 if (boundingBox[5] < oct->FX(
eLeft))
773 max = boundingBox[5];
777 max = oct->FX(
eLeft);
785 total += vol / 2.0 / (oct->GetDelta() * oct->GetDelta() *
786 oct->GetDelta() / 6.0 / sqrt(2));
801 : x1(p1), x2(p2), R(r), delta(d)
808 for (
int i = 0; i < 3; i++)
810 Le[i] = p[i] - x1[i];
811 Re[i] = p[i] - x2[i];
812 s[i] = x2[i] - x1[i];
815 dev[0] = Le[1] * Re[2] - Re[1] * Le[2];
816 dev[1] = Le[0] * Re[2] - Re[0] * Le[2];
817 dev[2] = Le[0] * Re[1] - Re[0] * Le[1];
820 sqrt(dev[0] * dev[0] + dev[1] * dev[1] + dev[2] * dev[2]) /
821 sqrt(s[0] * s[0] + s[1] * s[1] + s[2] * s[2]);
823 NekDouble t = -1.0 * ((x1[0] - p[0]) * s[0] + (x1[1] - p[1]) * s[1] +
824 (x1[1] - p[1]) * s[1]) /
827 if (dist < R && !(t > 1) && !(t < 0))
839 return sqrt((x1[0] - x2[0]) * (x1[0] - x2[0]) +
840 (x1[1] - x2[1]) * (x1[1] - x2[1]) +
841 (x1[2] - x2[2]) * (x1[2] - x2[2]));
845 void Octree::CompileCuravturePointList()
847 for (
int i = 1; i <= m_cad->GetNumSurf(); i++)
863 NekDouble du = (bounds[1] - bounds[0]) / (40 - 1);
864 NekDouble dv = (bounds[3] - bounds[2]) / (40 - 1);
871 for (
int j = 0; j < 40; j++)
873 for (
int k = 0; k < 40; k++)
876 uv[0] = k * du + bounds[0];
877 uv[1] = j * dv + bounds[2];
882 samplepoints[k][j] = surf->P(uv);
886 for (
int j = 0; j < 40 - 1; j++)
888 for (
int k = 0; k < 40 - 1; k++)
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]));
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]));
914 int nu = ceil(DeltaU / m_minDelta) * 40;
915 int nv = ceil(DeltaV / m_minDelta) * 40;
917 for (
int j = 0; j < nu; j++)
919 for (
int k = 0; k < nv; k++)
922 uv[0] = (bounds[1] - bounds[0]) / (nu - 1) * j + bounds[0];
923 uv[1] = (bounds[3] - bounds[2]) / (nv - 1) * k + bounds[2];
938 bool minlimited =
false;
942 2.0 * (1.0 / C) * sqrt(m_eps * (2.0 - m_eps));
944 if (del > m_maxDelta)
948 if (del < m_minDelta)
959 surf->GetId(), uv, surf->P(uv), del, ideal);
961 m_cpList.push_back(newCPoint);
967 surf->GetId(), uv, surf->P(uv), del);
969 m_cpList.push_back(newCPoint);
976 surf->GetId(), uv, surf->P(uv));
978 m_cpList.push_back(newCPoint);
984 if (m_udsfile ==
"N")
990 vector<linesource> lsources;
992 fle.open(m_udsfile.c_str());
998 getline(fle, fileline);
999 stringstream s(fileline);
1009 x1[0] = boost::lexical_cast<
double>(word);
1010 s >> x1[1] >> x1[2] >> x2[0] >> x2[1] >> x2[2] >> r >> d;
1012 lsources.push_back(
linesource(x1, x2, r, d));
1016 for (
int j = 0; j < lsources.size(); j++)
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;
1026 for (
int i = 0; i < m_cpList.size(); i++)
1028 for (
int j = 0; j < lsources.size(); j++)
1030 if (lsources[j].withinRange(m_cpList[i]->GetLoc()))
1033 m_cpList[i]->SetDelta(lsources[j].delta);
1068 return fabs(i->GetDelta() - j->GetDelta()) / i->Distance(j);
#define ASSERTL0(condition, msg)
Basic information about an element.
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.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Array< OneD, NekDouble > x2
ElementFactory & GetElementFactory()
Represents a point in the domain.
boost::shared_ptr< CurvaturePoint > CurvaturePointSharedPtr
boost::shared_ptr< Octant > OctantSharedPtr
linesource(Array< OneD, NekDouble > p1, Array< OneD, NekDouble > p2, NekDouble r, NekDouble d)
boost::shared_ptr< CADSurf > CADSurfSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
boost::shared_ptr< Element > ElementSharedPtr
bool withinRange(Array< OneD, NekDouble > p)