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

#include <BLMesh.h>

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

Classes

struct  blInfo
 

Public Types

typedef boost::shared_ptr< blInfoblInfoSharedPtr
 

Public Member Functions

 BLMesh (MeshSharedPtr m, std::vector< unsigned int > bls, NekDouble b, int l, NekDouble p, int id)
 default constructor More...
 
void Mesh ()
 Execute boundary layer meshing. More...
 
std::vector< unsigned int > GetSymSurfs ()
 
std::vector< unsigned int > GetBLSurfs ()
 
std::map< NodeSharedPtr,
NodeSharedPtr
GetSymNodes ()
 
std::vector< ElementSharedPtrGetPseudoSurface ()
 

Private Member Functions

void Setup ()
 
void GrowLayers ()
 
void Shrink ()
 
void BuildElements ()
 
bool TestIntersectionEl (ElementSharedPtr e1, ElementSharedPtr e2)
 
bool IsPrismValid (ElementSharedPtr el)
 
NekDouble Proximity (NodeSharedPtr n, ElementSharedPtr el)
 
NekDouble Visability (std::vector< ElementSharedPtr > tris, Array< OneD, NekDouble > N)
 
Array< OneD, NekDoubleGetNormal (std::vector< ElementSharedPtr > tris)
 

Private Attributes

MeshSharedPtr m_mesh
 mesh object containing surface mesh More...
 
std::vector< unsigned int > m_blsurfs
 List of surfaces onto which boundary layers are placed. More...
 
NekDouble m_bl
 thickness of the boundary layer More...
 
NekDouble m_prog
 
int m_layer
 
int m_id
 
Array< OneD, NekDoublem_layerT
 
std::vector< unsigned int > m_symSurfs
 list of surfaces to be remeshed due to the boundary layer More...
 
std::map< NodeSharedPtr,
blInfoSharedPtr
m_blData
 data structure used to store and develop bl information More...
 
std::map< NodeSharedPtr,
std::vector< blInfoSharedPtr > > 
m_nToNInfo
 
std::map< ElementSharedPtr,
ElementSharedPtr
m_priToTri
 
std::vector< ElementSharedPtrm_psuedoSurface
 
NekMatrix< NekDoublem_deriv [3]
 

Friends

class MemoryManager< BLMesh >
 

Detailed Description

Definition at line 49 of file BLMesh.h.

Member Typedef Documentation

Definition at line 105 of file BLMesh.h.

Constructor & Destructor Documentation

Nektar::NekMeshUtils::BLMesh::BLMesh ( MeshSharedPtr  m,
std::vector< unsigned int >  bls,
NekDouble  b,
int  l,
NekDouble  p,
int  id 
)
inline

default constructor

Definition at line 57 of file BLMesh.h.

59  : m_mesh(m), m_blsurfs(bls), m_bl(b), m_prog(p), m_layer(l), m_id(id)
60  {
61 
62  };
MeshSharedPtr m_mesh
mesh object containing surface mesh
Definition: BLMesh.h:121
NekDouble m_bl
thickness of the boundary layer
Definition: BLMesh.h:125
std::vector< unsigned int > m_blsurfs
List of surfaces onto which boundary layers are placed.
Definition: BLMesh.h:123

Member Function Documentation

void Nektar::NekMeshUtils::BLMesh::BuildElements ( )
private

Definition at line 753 of file BLMesh.cpp.

References Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::StdRegions::eBackwards, Nektar::StdRegions::eForwards, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::eTriangle, Nektar::StdRegions::find(), Nektar::NekMeshUtils::GetElementFactory(), and Nektar::iterator.

754 {
755  // make prisms
756  map<CADOrientation::Orientation, vector<int> > baseTri;
757  map<CADOrientation::Orientation, vector<int> > topTri;
758 
759  vector<int> tmp;
760  // back-base
761  tmp.push_back(0);
762  tmp.push_back(4);
763  tmp.push_back(1);
764  baseTri[CADOrientation::eBackwards] = tmp;
765  tmp.clear();
766  // for-base
767  tmp.push_back(0);
768  tmp.push_back(1);
769  tmp.push_back(4);
770  baseTri[CADOrientation::eForwards] = tmp;
771  // back-top
772  tmp.clear();
773  tmp.push_back(3);
774  tmp.push_back(5);
775  tmp.push_back(2);
776  topTri[CADOrientation::eBackwards] = tmp;
777  // for-top
778  tmp.clear();
779  tmp.push_back(3);
780  tmp.push_back(2);
781  tmp.push_back(5);
782  topTri[CADOrientation::eForwards] = tmp;
783 
784  ElmtConfig pconf(LibUtilities::ePrism, 1, false, false);
785  ElmtConfig tconf(LibUtilities::eTriangle, 1, false, false);
786 
787  for (int i = 0; i < m_mesh->m_element[2].size(); i++)
788  {
789  ElementSharedPtr el = m_mesh->m_element[2][i];
791  find(m_blsurfs.begin(), m_blsurfs.end(), el->m_parentCAD->GetId());
792 
793  if (f == m_blsurfs.end())
794  {
795  // if this triangle is not in bl surfs continue
796  continue;
797  }
798 
799  vector<NodeSharedPtr> tn(3); // nodes for pseduo surface
800  vector<NodeSharedPtr> pn(6); // all prism nodes
801  vector<NodeSharedPtr> n = el->GetVertexList();
802  CADOrientation::Orientation o = el->m_parentCAD->Orientation();
803 
804  for (int j = 0; j < 3; j++)
805  {
806  pn[baseTri[o][j]] = n[j];
807  pn[topTri[o][j]] = m_blData[n[j]]->pNode;
808  tn[j] = m_blData[n[j]]->pNode;
809  }
810 
811  vector<int> tags;
812  tags.push_back(m_id);
814  LibUtilities::ePrism, pconf, pn, tags);
815  E->SetId(i);
816 
817  m_mesh->m_element[3].push_back(E);
818 
819  // tag of this element doesnt matter so can just be 1
821  LibUtilities::eTriangle, tconf, tn, tags);
822  m_psuedoSurface.push_back(T);
823 
824  T->m_parentCAD = el->m_parentCAD;
825 
826  m_priToTri[E] = el;
827  }
828 }
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
std::vector< ElementSharedPtr > m_psuedoSurface
Definition: BLMesh.h:137
MeshSharedPtr m_mesh
mesh object containing surface mesh
Definition: BLMesh.h:121
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
std::map< NodeSharedPtr, blInfoSharedPtr > m_blData
data structure used to store and develop bl information
Definition: BLMesh.h:133
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:316
std::map< ElementSharedPtr, ElementSharedPtr > m_priToTri
Definition: BLMesh.h:136
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
std::vector< unsigned int > m_blsurfs
List of surfaces onto which boundary layers are placed.
Definition: BLMesh.h:123
std::vector<unsigned int> Nektar::NekMeshUtils::BLMesh::GetBLSurfs ( )
inline

Definition at line 73 of file BLMesh.h.

References m_blsurfs.

74  {
75  return m_blsurfs;
76  }
std::vector< unsigned int > m_blsurfs
List of surfaces onto which boundary layers are placed.
Definition: BLMesh.h:123
Array< OneD, NekDouble > Nektar::NekMeshUtils::BLMesh::GetNormal ( std::vector< ElementSharedPtr tris)
private

Definition at line 844 of file BLMesh.cpp.

845 {
846  // compile list of normals
847  vector<Array<OneD, NekDouble> > N;
848  for (int i = 0; i < tris.size(); i++)
849  {
850  N.push_back(tris[i]->Normal(true));
851  }
852 
853  vector<NekDouble> w(N.size());
854  Array<OneD, NekDouble> Np(3, 0.0);
855 
856  for (int i = 0; i < N.size(); i++)
857  {
858  w[i] = 1.0 / N.size();
859  }
860 
861  for (int i = 0; i < N.size(); i++)
862  {
863  Np[0] += w[i] * N[i][0];
864  Np[1] += w[i] * N[i][1];
865  Np[2] += w[i] * N[i][2];
866  }
867  NekDouble mag = sqrt(Np[0] * Np[0] + Np[1] * Np[1] + Np[2] * Np[2]);
868  Np[0] /= mag;
869  Np[1] /= mag;
870  Np[2] /= mag;
871 
872  Array<OneD, NekDouble> Ninital = Np;
873 
874  NekDouble dot = 0.0;
875  int ct = 0;
876  vector<NekDouble> a(N.size());
877  while (fabs(dot - 1) > 1e-6)
878  {
879  ct++;
880  Array<OneD, NekDouble> Nplast(3);
881  Nplast[0] = Np[0];
882  Nplast[1] = Np[1];
883  Nplast[2] = Np[2];
884 
885  NekDouble aSum = 0.0;
886  for (int i = 0; i < N.size(); i++)
887  {
888  NekDouble dot2 =
889  Np[0] * N[i][0] + Np[1] * N[i][1] + Np[2] * N[i][2];
890  if (fabs(dot2 - 1) < 1e-9)
891  {
892  a[i] = dot2 / fabs(dot2) * 1e-9;
893  }
894  else
895  {
896  a[i] = acos(dot2);
897  }
898 
899  aSum += a[i];
900  }
901 
902  NekDouble wSum = 0.0;
903  for (int i = 0; i < N.size(); i++)
904  {
905  w[i] = w[i] * a[i] / aSum;
906  wSum += w[i];
907  }
908 
909  for (int i = 0; i < N.size(); i++)
910  {
911  w[i] = w[i] / wSum;
912  }
913 
914  Array<OneD, NekDouble> NpN(3, 0.0);
915  for (int i = 0; i < N.size(); i++)
916  {
917  NpN[0] += w[i] * N[i][0];
918  NpN[1] += w[i] * N[i][1];
919  NpN[2] += w[i] * N[i][2];
920  }
921  mag = sqrt(NpN[0] * NpN[0] + NpN[1] * NpN[1] + NpN[2] * NpN[2]);
922  NpN[0] /= mag;
923  NpN[1] /= mag;
924  NpN[2] /= mag;
925 
926  Np[0] = 0.8 * NpN[0] + (1.0 - 0.8) * Np[0];
927  Np[1] = 0.8 * NpN[1] + (1.0 - 0.8) * Np[1];
928  Np[2] = 0.8 * NpN[2] + (1.0 - 0.8) * Np[2];
929  mag = sqrt(Np[0] * Np[0] + Np[1] * Np[1] + Np[2] * Np[2]);
930  Np[0] /= mag;
931  Np[1] /= mag;
932  Np[2] /= mag;
933 
934  dot = Np[0] * Nplast[0] + Np[1] * Nplast[1] + Np[2] * Nplast[2];
935 
936  if (ct > 100000)
937  {
938  cout << "run out of iterations" << endl;
939  Np = Ninital;
940  break;
941  }
942  }
943 
944  return Np;
945 }
double NekDouble
std::vector<ElementSharedPtr> Nektar::NekMeshUtils::BLMesh::GetPseudoSurface ( )
inline

Definition at line 80 of file BLMesh.h.

References m_psuedoSurface.

81  {
82  return m_psuedoSurface;
83  }
std::vector< ElementSharedPtr > m_psuedoSurface
Definition: BLMesh.h:137
map< NodeSharedPtr, NodeSharedPtr > Nektar::NekMeshUtils::BLMesh::GetSymNodes ( )

Definition at line 159 of file BLMesh.cpp.

References Nektar::iterator.

160 {
161  map<NodeSharedPtr, NodeSharedPtr> ret;
162 
164  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
165  {
166  if (!bit->second->onSym)
167  {
168  continue;
169  }
170  CADSurfSharedPtr s = m_mesh->m_cad->GetSurf(bit->second->symsurf);
171  Array<OneD, NekDouble> loc = bit->second->pNode->GetLoc();
172  Array<OneD, NekDouble> uv(2);
173  uv = s->locuv(loc);
174  bit->second->pNode->SetCADSurf(bit->second->symsurf, s, uv);
175  ret[bit->first] = bit->second->pNode;
176  }
177  return ret;
178 }
MeshSharedPtr m_mesh
mesh object containing surface mesh
Definition: BLMesh.h:121
std::map< NodeSharedPtr, blInfoSharedPtr > m_blData
data structure used to store and develop bl information
Definition: BLMesh.h:133
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:172
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::vector<unsigned int> Nektar::NekMeshUtils::BLMesh::GetSymSurfs ( )
inline

Definition at line 69 of file BLMesh.h.

References m_symSurfs.

70  {
71  return m_symSurfs;
72  }
std::vector< unsigned int > m_symSurfs
list of surfaces to be remeshed due to the boundary layer
Definition: BLMesh.h:131
void Nektar::NekMeshUtils::BLMesh::GrowLayers ( )
private

Definition at line 197 of file BLMesh.cpp.

References Nektar::StdRegions::find(), Nektar::NekMeshUtils::GetBox(), Nektar::NekMeshUtils::Infont(), and Nektar::iterator.

198 {
200 
201  // setup up a tree which is formed of boxes of each surface plus some
202  // extra room (ideal bl thick)
203 
204  // in each iteration a tree is made for all the triangles in each surface
205  // when considering to stop a bounary layer growing, it first
206  // looks at the top tree to find surfaces which are canditates
207  // it then searches the subtrees for each triangle canditate
208  // it then does a distance calcation.
209  // if a boundary layer should be close to that from another surface, it
210  // should stop
211 
212  map<int, vector<ElementSharedPtr> > psElements;
213  for (int i = 0; i < m_mesh->m_element[2].size(); i++)
214  {
215  ElementSharedPtr el = m_mesh->m_element[2][i];
217  find(m_blsurfs.begin(), m_blsurfs.end(), el->m_parentCAD->GetId());
218 
220  m_symSurfs.begin(), m_symSurfs.end(), el->m_parentCAD->GetId());
221 
222  if (f == m_blsurfs.end() && s == m_symSurfs.end())
223  {
224  psElements[el->m_parentCAD->GetId()].push_back(el);
225  }
226  }
227  for (int i = 0; i < m_psuedoSurface.size(); i++)
228  {
229  psElements[m_psuedoSurface[i]->m_parentCAD->GetId()].push_back(
230  m_psuedoSurface[i]);
231  }
232 
233  bgi::rtree<boxI, bgi::quadratic<16> > TopTree;
234  map<int, bgi::rtree<boxI, bgi::quadratic<16> > > SubTrees;
235 
236  // ofstream file;
237  // file.open("pts.3D");
238  // file << "x y z value" << endl;
239 
240  for (int l = 1; l < m_layer; l++)
241  {
242  NekDouble delta = (m_layerT[l] - m_layerT[l - 1]);
243  TopTree.clear();
244  SubTrees.clear();
245  map<int, vector<ElementSharedPtr> >::iterator it;
246  for (it = psElements.begin(); it != psElements.end(); it++)
247  {
248  TopTree.insert(make_pair(GetBox(it->second, m_bl), it->first));
249  vector<boxI> toInsert;
250  for (int i = 0; i < it->second.size(); i++)
251  {
252  toInsert.push_back(make_pair(GetBox(it->second[i], m_bl), i));
253  }
254  SubTrees[it->first].insert(toInsert.begin(), toInsert.end());
255  }
256 
257  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
258  {
259  if (bit->second->stopped)
260  {
261  continue;
262  }
263 
264  vector<boxI> results;
265  TopTree.query(bgi::intersects(point(bit->second->pNode->m_x,
266  bit->second->pNode->m_y,
267  bit->second->pNode->m_z)),
268  back_inserter(results));
269  set<int> surfs;
270  for (int i = 0; i < results.size(); i++)
271  {
273  bit->second->surfs.find(results[i].second);
274  if (f == bit->second->surfs.end())
275  {
276  // hit
277  surfs.insert(results[i].second);
278  }
279  }
280 
281  set<int>::iterator iit;
282  bool hit = false;
283  for (iit = surfs.begin(); iit != surfs.end(); iit++)
284  {
285  results.clear();
286  SubTrees[*iit].query(
287  bgi::intersects(GetBox(bit->second->pNode, m_bl)),
288  back_inserter(results));
289  for (int i = 0; i < results.size(); i++)
290  {
291  if (Infont(bit->second->pNode,
292  psElements[*iit][results[i].second]))
293  {
294  NekDouble prox =
295  Proximity(bit->second->pNode,
296  psElements[*iit][results[i].second]);
297  if (prox < delta * 2.5)
298  {
299  hit = true;
300  // cout << "hit" << endl;
301  bit->second->stopped = true;
302  /*file << bit->first->m_x << " " << bit->first->m_y
303  << " " << bit->first->m_z << " " << l << endl;
304  file << bit->second->pNode->m_x << " " <<
305  bit->second->pNode->m_y << " " <<
306  bit->second->pNode->m_z << " " << l << endl;
307  m_mesh->m_element[2].clear();
308  m_mesh->m_expDim--;
309  m_mesh->m_element[2].push_back(psElements[*iit][results[i].second]);*/
310  break;
311  // return;
312  }
313  }
314  }
315  if (hit)
316  break;
317  }
318  }
319 
320  // if after proximity scanning all is okay, advance the layer
321  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
322  {
323  if (bit->second->stopped)
324  {
325  continue;
326  }
327 
328  // test the smoothness
329  bool shouldStop = false;
330  vector<blInfoSharedPtr> ne = m_nToNInfo[bit->first];
331  for (int i = 0; i < ne.size(); i++)
332  {
333  if (ne[i]->bl < bit->second->bl)
334  {
335  shouldStop = true;
336  break;
337  }
338  }
339  if (shouldStop)
340  {
341  bit->second->stopped = true;
342  continue;
343  }
344 
345  bit->second->AlignNode(m_layerT[l]);
346  bit->second->bl = l;
347  }
348  }
349  // file.close();
350 }
std::vector< unsigned int > m_symSurfs
list of surfaces to be remeshed due to the boundary layer
Definition: BLMesh.h:131
std::vector< ElementSharedPtr > m_psuedoSurface
Definition: BLMesh.h:137
MeshSharedPtr m_mesh
mesh object containing surface mesh
Definition: BLMesh.h:121
Array< OneD, NekDouble > m_layerT
Definition: BLMesh.h:129
NekDouble m_bl
thickness of the boundary layer
Definition: BLMesh.h:125
bool Infont(NodeSharedPtr n, ElementSharedPtr el)
Definition: BLMesh.cpp:180
std::map< NodeSharedPtr, blInfoSharedPtr > m_blData
data structure used to store and develop bl information
Definition: BLMesh.h:133
NekDouble Proximity(NodeSharedPtr n, ElementSharedPtr el)
Definition: BLMesh.cpp:362
double NekDouble
bg::model::point< double, 3, bg::cs::cartesian > point
Definition: BLMesh.cpp:54
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::map< NodeSharedPtr, std::vector< blInfoSharedPtr > > m_nToNInfo
Definition: BLMesh.h:135
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:316
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
box GetBox(ElementSharedPtr el, NekDouble ov)
Definition: BLMesh.cpp:64
std::vector< unsigned int > m_blsurfs
List of surfaces onto which boundary layers are placed.
Definition: BLMesh.h:123
bool Nektar::NekMeshUtils::BLMesh::IsPrismValid ( ElementSharedPtr  el)
private

Definition at line 698 of file BLMesh.cpp.

References Nektar::eFULL.

699 {
700  NekDouble mn = numeric_limits<double>::max();
701  NekDouble mx = -1.0 * numeric_limits<double>::max();
702  vector<NodeSharedPtr> ns = el->GetVertexList();
703  NekVector<NekDouble> X(6), Y(6), Z(6);
704  for (int j = 0; j < ns.size(); j++)
705  {
706  X(j) = ns[j]->m_x;
707  Y(j) = ns[j]->m_y;
708  Z(j) = ns[j]->m_z;
709  }
710  NekVector<NekDouble> x1(6), y1(6), z1(6), x2(6), y2(6), z2(6), x3(6), y3(6),
711  z3(6);
712 
713  x1 = m_deriv[0] * X;
714  y1 = m_deriv[0] * Y;
715  z1 = m_deriv[0] * Z;
716  x2 = m_deriv[1] * X;
717  y2 = m_deriv[1] * Y;
718  z2 = m_deriv[1] * Z;
719  x3 = m_deriv[2] * X;
720  y3 = m_deriv[2] * Y;
721  z3 = m_deriv[2] * Z;
722 
723  for (int j = 0; j < 6; j++)
724  {
725  DNekMat dxdz(3, 3, 1.0, eFULL);
726  dxdz(0, 0) = x1(j);
727  dxdz(0, 1) = x2(j);
728  dxdz(0, 2) = x3(j);
729  dxdz(1, 0) = y1(j);
730  dxdz(1, 1) = y2(j);
731  dxdz(1, 2) = y3(j);
732  dxdz(2, 0) = z1(j);
733  dxdz(2, 1) = z2(j);
734  dxdz(2, 2) = z3(j);
735 
736  NekDouble jacDet =
737  dxdz(0, 0) * (dxdz(1, 1) * dxdz(2, 2) - dxdz(2, 1) * dxdz(1, 2)) -
738  dxdz(0, 1) * (dxdz(1, 0) * dxdz(2, 2) - dxdz(2, 0) * dxdz(1, 2)) +
739  dxdz(0, 2) * (dxdz(1, 0) * dxdz(2, 1) - dxdz(2, 0) * dxdz(1, 1));
740  mn = min(mn, jacDet);
741  mx = max(mx, jacDet);
742  }
743 
744  /*SpatialDomains::GeometrySharedPtr geom = el->GetGeom(3);
745  SpatialDomains::GeomFactorsSharedPtr gfac = geom->GetGeomFactors();
746 
747  cout << mn << " " << mx << " " << (mn > 0) << " " << gfac->IsValid() <<
748  endl;*/
749 
750  return mn > 0;
751 }
NekMatrix< NekDouble > m_deriv[3]
Definition: BLMesh.h:138
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:52
double NekDouble
void Nektar::NekMeshUtils::BLMesh::Mesh ( )

Execute boundary layer meshing.

Definition at line 121 of file BLMesh.cpp.

References Nektar::iterator.

122 {
123  Setup();
124 
125  BuildElements();
126 
127  GrowLayers();
128 
129  Shrink();
130 
132  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
133  {
134  vector<blInfoSharedPtr> infos = m_nToNInfo[bit->first];
135  for (int i = 0; i < infos.size(); i++)
136  {
137  if (bit->second->bl > infos[i]->bl + 1)
138  {
139  cout << "non smooth error " << bit->second->bl << " "
140  << infos[i]->bl << endl;
141  }
142  }
143  }
144 
145  for (int i = 0; i < m_mesh->m_element[3].size(); i++)
146  {
147  ElementSharedPtr el = m_mesh->m_element[3][i];
148  if (!IsPrismValid(el))
149  {
150  cout << "validity error " << el->GetId() << endl;
151  }
152  }
153 
154  /*m_mesh->m_element[2] = m_psuedoSurface;
155  m_mesh->m_element[3].clear();
156  m_mesh->m_expDim = 2;*/
157 }
MeshSharedPtr m_mesh
mesh object containing surface mesh
Definition: BLMesh.h:121
std::map< NodeSharedPtr, blInfoSharedPtr > m_blData
data structure used to store and develop bl information
Definition: BLMesh.h:133
bool IsPrismValid(ElementSharedPtr el)
Definition: BLMesh.cpp:698
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::map< NodeSharedPtr, std::vector< blInfoSharedPtr > > m_nToNInfo
Definition: BLMesh.h:135
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
NekDouble Nektar::NekMeshUtils::BLMesh::Proximity ( NodeSharedPtr  n,
ElementSharedPtr  el 
)
private

Definition at line 362 of file BLMesh.cpp.

References Vmath::Dot(), class_topology::Node, class_topology::P, and Vmath::Vsub().

363 {
364  vector<NodeSharedPtr> ns = el->GetVertexList();
365  Array<OneD, NekDouble> B = ns[0]->GetLoc();
366  Array<OneD, NekDouble> E0(3);
367  E0[0] = ns[1]->m_x - ns[0]->m_x;
368  E0[1] = ns[1]->m_y - ns[0]->m_y;
369  E0[2] = ns[1]->m_z - ns[0]->m_z;
370  Array<OneD, NekDouble> E1(3);
371  E1[0] = ns[2]->m_x - ns[0]->m_x;
372  E1[1] = ns[2]->m_y - ns[0]->m_y;
373  E1[2] = ns[2]->m_z - ns[0]->m_z;
374  Array<OneD, NekDouble> P = n->GetLoc();
375 
376  NekDouble a = Dot(E0, E0);
377  NekDouble b = Dot(E0, E1);
378  NekDouble c = Dot(E1, E1);
379 
380  Array<OneD, NekDouble> BP(3);
381  Vmath::Vsub(3, B, 1, P, 1, BP, 1);
382 
383  NekDouble d = Dot(E0, BP);
384  NekDouble e = Dot(E1, BP);
385 
386  NekDouble det = a * c - b * b;
387  NekDouble s = b * e - c * d, t = b * d - a * e;
388 
389  if (s + t <= det)
390  {
391  if (s < 0)
392  {
393  if (t < 0)
394  {
395  t = 0;
396  s = 0;
397  }
398  else
399  {
400  s = 0;
401  t = (e >= 0 ? 0 : (-e >= c ? 1 : -e / c));
402  }
403  }
404  else if (t < 0)
405  {
406  t = 0;
407  s = (d >= 0 ? 0 : (-d >= a ? 1 : -d / a));
408  }
409  else
410  {
411  NekDouble invDet = 1.0 / det;
412  s *= invDet;
413  t *= invDet;
414  }
415  }
416  else
417  {
418  if (s < 0)
419  {
420  s = 0;
421  t = 1;
422  }
423  else if (t < 0)
424  {
425  s = 1;
426  t = 0;
427  }
428  else
429  {
430  NekDouble numer = c + e - b - d;
431  if (numer <= 0)
432  {
433  s = 0;
434  }
435  else
436  {
437  NekDouble denom = a - 2 * b + c;
438  s = (numer >= denom ? 1 : numer / denom);
439  }
440  t = 1 - s;
441  }
442  }
443 
444  NodeSharedPtr point = boost::shared_ptr<Node>(
445  new Node(0, B[0] + s * E0[0] + t * E1[0], B[1] + s * E0[1] + t * E1[1],
446  B[2] + s * E0[2] + t * E1[2]));
447 
448  return n->Distance(point);
449 }
NekDouble Dot(Array< OneD, NekDouble > a, Array< OneD, NekDouble > b)
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
double NekDouble
bg::model::point< double, 3, bg::cs::cartesian > point
Definition: BLMesh.cpp:54
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:343
void Nektar::NekMeshUtils::BLMesh::Setup ( )
private

Definition at line 947 of file BLMesh.cpp.

References ASSERTL0, Nektar::LibUtilities::eNodalPrismElec, Nektar::StdRegions::find(), Nektar::LibUtilities::NodalUtil::GetVandermonde(), Nektar::LibUtilities::NodalUtil::GetVandermondeForDeriv(), Nektar::iterator, class_topology::Node, and Nektar::LibUtilities::PointsManager().

948 {
949  NekDouble a = 2.0 * (1.0 - m_prog) / (1.0 - pow(m_prog, m_layer + 1));
950  m_layerT = Array<OneD, NekDouble>(m_layer);
951  m_layerT[0] = a * m_prog * m_bl;
952  for (int i = 1; i < m_layer; i++)
953  {
954  m_layerT[i] = m_layerT[i - 1] + a * pow(m_prog, i) * m_bl;
955  }
956 
957  if(m_mesh->m_verbose)
958  {
959  cout << "First layer height " << m_layerT[0] << endl;
960  }
961 
962  // this sets up all the boundary layer normals data holder
963  set<int> symSurfs;
965  int ct = 0;
966  int failed = 0;
967 
968  // ofstream file1;
969  // file1.open("pts.3D");
970  // file1 << "X Y Z value" << endl;
971  for (it = m_mesh->m_vertexSet.begin(); it != m_mesh->m_vertexSet.end();
972  it++, ct++)
973  {
974  vector<pair<int, CADSurfSharedPtr> > ss = (*it)->GetCADSurfs();
975  vector<unsigned int> surfs;
976  for (int i = 0; i < ss.size(); i++)
977  {
978  surfs.push_back(ss[i].first);
979  }
980  sort(surfs.begin(), surfs.end());
981  vector<unsigned int> inter, diff;
982 
983  set_intersection(m_blsurfs.begin(), m_blsurfs.end(), surfs.begin(),
984  surfs.end(), back_inserter(inter));
985  set_symmetric_difference(inter.begin(), inter.end(), surfs.begin(),
986  surfs.end(), back_inserter(diff));
987 
988  // is somewhere on a bl surface
989  if (inter.size() > 0)
990  {
991  // initialise a new bl boudnary node
992  blInfoSharedPtr bln = boost::shared_ptr<blInfo>(new blInfo);
993  bln->oNode = (*it);
994  bln->stopped = false;
995 
996  // file1 << (*it)->m_x << " " << (*it)->m_y << " " << (*it)->m_z <<
997  // " " << ss.size() << endl;
998 
999  if (diff.size() > 0)
1000  {
1001  // if the diff size is greater than 1 there is a curve that
1002  // needs remeshing
1003  ASSERTL0(diff.size() <= 1, "not setup for curve bl refinement");
1004  symSurfs.insert(diff[0]);
1005  bln->symsurf = diff[0];
1006  bln->onSym = true;
1007  }
1008  else
1009  {
1010  bln->onSym = false;
1011  }
1012 
1013  m_blData[(*it)] = bln;
1014  }
1015  }
1016  // file1.close();
1017 
1018  // need a map from vertex idx to surface elements
1019  // but do not care about triangles which are not in the bl
1020  for (int i = 0; i < m_mesh->m_element[2].size(); i++)
1021  {
1023  find(m_blsurfs.begin(), m_blsurfs.end(),
1024  m_mesh->m_element[2][i]->m_parentCAD->GetId());
1025 
1026  if (f == m_blsurfs.end())
1027  {
1028  // if this triangle is not in bl surfs continue
1029  continue;
1030  }
1031 
1032  vector<NodeSharedPtr> ns = m_mesh->m_element[2][i]->GetVertexList();
1033  for (int j = 0; j < ns.size(); j++)
1034  {
1035  m_blData[ns[j]]->els.push_back(m_mesh->m_element[2][i]);
1036  m_blData[ns[j]]->surfs.insert(
1037  m_mesh->m_element[2][i]->m_parentCAD->GetId());
1038  }
1039  }
1040 
1042  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1043  {
1044  // calculate mesh normal
1045  bit->second->N = GetNormal(bit->second->els);
1046 
1047  if (Visability(bit->second->els, bit->second->N) < 0.0)
1048  {
1049  cerr << "failed " << bit->first->m_x << " " << bit->first->m_y
1050  << " " << bit->first->m_z << " "
1051  << Visability(bit->second->els, bit->second->N) << endl;
1052  failed++;
1053  }
1054 
1055  Array<OneD, NekDouble> loc = bit->first->GetLoc();
1056  for (int k = 0; k < 3; k++)
1057  {
1058  loc[k] += bit->second->N[k] * m_layerT[0];
1059  }
1060 
1061  bit->second->pNode = boost::shared_ptr<Node>(
1062  new Node(m_mesh->m_numNodes++, loc[0], loc[1], loc[2]));
1063  bit->second->bl = 0;
1064  }
1065 
1066  m_symSurfs = vector<unsigned int>(symSurfs.begin(), symSurfs.end());
1067 
1068  // now need to enforce that all symmetry plane nodes have their normal
1069  // forced onto the symmetry surface
1070  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1071  {
1072  if (!bit->second->onSym)
1073  {
1074  continue;
1075  }
1076 
1077  Array<OneD, NekDouble> uv(2);
1078  Array<OneD, NekDouble> loc = bit->second->pNode->GetLoc();
1079  m_mesh->m_cad->GetSurf(bit->second->symsurf)->ProjectTo(loc, uv);
1080 
1081  Array<OneD, NekDouble> nl =
1082  m_mesh->m_cad->GetSurf(bit->second->symsurf)->P(uv);
1083 
1084  Array<OneD, NekDouble> N(3);
1085  N[0] = nl[0] - bit->first->m_x;
1086  N[1] = nl[1] - bit->first->m_y;
1087  N[2] = nl[2] - bit->first->m_z;
1088 
1089  NekDouble mag = sqrt(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
1090  N[0] /= mag;
1091  N[1] /= mag;
1092  N[2] /= mag;
1093 
1094  bit->second->N = N;
1095  bit->second->AlignNode(m_layerT[0]);
1096  }
1097 
1098  // now smooth all the normals by distance weighted average
1099  // keep normals on curves constant
1100  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1101  {
1102  set<int> added;
1103  added.insert(bit->first->m_id);
1104  for (int i = 0; i < bit->second->els.size(); i++)
1105  {
1106  vector<NodeSharedPtr> ns = bit->second->els[i]->GetVertexList();
1107  for (int j = 0; j < ns.size(); j++)
1108  {
1109  set<int>::iterator t = added.find(ns[j]->m_id);
1110  if (t == added.end())
1111  {
1112  m_nToNInfo[bit->first].push_back(m_blData[ns[j]]);
1113  }
1114  }
1115  }
1116  }
1117 
1118  for (int l = 0; l < 10; l++)
1119  {
1120  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1121  {
1122  if (bit->first->GetNumCADSurf() > 1)
1123  {
1124  continue;
1125  }
1126 
1127  Array<OneD, NekDouble> sumV(3, 0.0);
1128  vector<blInfoSharedPtr> data = m_nToNInfo[bit->first];
1129  NekDouble Dtotal = 0.0;
1130  for (int i = 0; i < data.size(); i++)
1131  {
1132  NekDouble d = bit->first->Distance(data[i]->oNode);
1133  Dtotal += d;
1134  sumV[0] += data[i]->N[0] / d;
1135  sumV[1] += data[i]->N[1] / d;
1136  sumV[2] += data[i]->N[2] / d;
1137  }
1138  sumV[0] *= Dtotal;
1139  sumV[1] *= Dtotal;
1140  sumV[2] *= Dtotal;
1141  NekDouble mag =
1142  sqrt(sumV[0] * sumV[0] + sumV[1] * sumV[1] + sumV[2] * sumV[2]);
1143  sumV[0] /= mag;
1144  sumV[1] /= mag;
1145  sumV[2] /= mag;
1146 
1147  Array<OneD, NekDouble> N(3);
1148 
1149  N[0] = (1.0 - 0.8) * bit->second->N[0] + 0.8 * sumV[0];
1150  N[1] = (1.0 - 0.8) * bit->second->N[1] + 0.8 * sumV[1];
1151  N[2] = (1.0 - 0.8) * bit->second->N[2] + 0.8 * sumV[2];
1152 
1153  mag = sqrt(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
1154  N[0] /= mag;
1155  N[1] /= mag;
1156  N[2] /= mag;
1157 
1158  bit->second->N = N;
1159  bit->second->AlignNode(m_layerT[0]);
1160  }
1161  }
1162 
1163  /*ofstream file;
1164  file.open("bl.lines");
1165  for(bit = m_blData.begin(); bit != m_blData.end(); bit++)
1166  {
1167  NekDouble l = 0.05;
1168  file << bit->first->m_x << ", " << bit->first->m_y << ", " <<
1169  bit->first->m_z << endl;
1170  file << bit->first->m_x + bit->second->N[0]*l << ", "
1171  << bit->first->m_y + bit->second->N[1]*l << ", "
1172  << bit->first->m_z + bit->second->N[2]*l << endl;
1173  file << endl;
1174  }
1175  file.close();*/
1176 
1177  ASSERTL0(failed == 0, "some normals failed to generate");
1178 
1179  LibUtilities::PointsKey pkey1(2, LibUtilities::eNodalPrismElec);
1180 
1181  Array<OneD, NekDouble> u1, v1, w1;
1182  LibUtilities::PointsManager()[pkey1]->GetPoints(u1, v1, w1);
1183 
1184  LibUtilities::NodalUtilPrism nodalPrism(1, u1, v1, w1);
1185 
1186  NekMatrix<NekDouble> Vandermonde = *nodalPrism.GetVandermonde();
1187  NekMatrix<NekDouble> VandermondeI = Vandermonde;
1188  VandermondeI.Invert();
1189 
1190  m_deriv[0] = *nodalPrism.GetVandermondeForDeriv(0) * VandermondeI;
1191  m_deriv[1] = *nodalPrism.GetVandermondeForDeriv(1) * VandermondeI;
1192  m_deriv[2] = *nodalPrism.GetVandermondeForDeriv(2) * VandermondeI;
1193 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Array< OneD, NekDouble > GetNormal(std::vector< ElementSharedPtr > tris)
Definition: BLMesh.cpp:844
std::vector< unsigned int > m_symSurfs
list of surfaces to be remeshed due to the boundary layer
Definition: BLMesh.h:131
MeshSharedPtr m_mesh
mesh object containing surface mesh
Definition: BLMesh.h:121
NekDouble Visability(std::vector< ElementSharedPtr > tris, Array< OneD, NekDouble > N)
Definition: BLMesh.cpp:830
Array< OneD, NekDouble > m_layerT
Definition: BLMesh.h:129
NekDouble m_bl
thickness of the boundary layer
Definition: BLMesh.h:125
std::map< NodeSharedPtr, blInfoSharedPtr > m_blData
data structure used to store and develop bl information
Definition: BLMesh.h:133
NekMatrix< NekDouble > m_deriv[3]
Definition: BLMesh.h:138
PointsManagerT & PointsManager(void)
double NekDouble
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::map< NodeSharedPtr, std::vector< blInfoSharedPtr > > m_nToNInfo
Definition: BLMesh.h:135
boost::shared_ptr< blInfo > blInfoSharedPtr
Definition: BLMesh.h:105
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:316
3D electrostatically spaced points on a Prism
Definition: PointsType.h:76
std::vector< unsigned int > m_blsurfs
List of surfaces onto which boundary layers are placed.
Definition: BLMesh.h:123
void Nektar::NekMeshUtils::BLMesh::Shrink ( )
private

Definition at line 622 of file BLMesh.cpp.

References ASSERTL0, and Nektar::iterator.

623 {
625  bool smsh = true;
626  while (smsh)
627  {
628  smsh = false;
629 
630  vector<ElementSharedPtr> inv;
631  for (int i = 0; i < m_mesh->m_element[3].size(); i++)
632  {
633  ElementSharedPtr el = m_mesh->m_element[3][i];
634  if (!IsPrismValid(el))
635  {
636  inv.push_back(el);
637  }
638  }
639 
640  smsh = (inv.size() > 0);
641 
642  for (int i = 0; i < inv.size(); i++)
643  {
644  ElementSharedPtr t = m_priToTri[inv[i]];
645  vector<blInfoSharedPtr> bls;
646  vector<NodeSharedPtr> ns = t->GetVertexList();
647  for (int j = 0; j < ns.size(); j++)
648  {
649  bls.push_back(m_blData[ns[j]]);
650  }
651  bool repeat = true;
652  while (repeat)
653  {
654  repeat = false;
655  int mx = 0;
656  for (int j = 0; j < 3; j++)
657  {
658  mx = max(mx, bls[j]->bl);
659  }
660  ASSERTL0(mx > 0, "shrinking to nothing");
661  for (int j = 0; j < 3; j++)
662  {
663  if (bls[j]->bl < mx)
664  {
665  continue;
666  }
667  bls[j]->bl--;
668  bls[j]->AlignNode(m_layerT[bls[j]->bl]);
669  }
670  if (!IsPrismValid(inv[i]))
671  {
672  repeat = true;
673  }
674  }
675  }
676 
677  bool repeat = true;
678  while (repeat)
679  {
680  repeat = false;
681  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
682  {
683  vector<blInfoSharedPtr> infos = m_nToNInfo[bit->first];
684  for (int i = 0; i < infos.size(); i++)
685  {
686  if (bit->second->bl > infos[i]->bl + 1)
687  {
688  bit->second->bl--;
689  bit->second->AlignNode(m_layerT[bit->second->bl]);
690  repeat = true;
691  }
692  }
693  }
694  }
695  }
696 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
MeshSharedPtr m_mesh
mesh object containing surface mesh
Definition: BLMesh.h:121
Array< OneD, NekDouble > m_layerT
Definition: BLMesh.h:129
std::map< NodeSharedPtr, blInfoSharedPtr > m_blData
data structure used to store and develop bl information
Definition: BLMesh.h:133
bool IsPrismValid(ElementSharedPtr el)
Definition: BLMesh.cpp:698
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::map< NodeSharedPtr, std::vector< blInfoSharedPtr > > m_nToNInfo
Definition: BLMesh.h:135
std::map< ElementSharedPtr, ElementSharedPtr > m_priToTri
Definition: BLMesh.h:136
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
bool Nektar::NekMeshUtils::BLMesh::TestIntersectionEl ( ElementSharedPtr  e1,
ElementSharedPtr  e2 
)
private

Definition at line 451 of file BLMesh.cpp.

References sign.

452 {
453  vector<NodeSharedPtr> ns1 = e1->GetVertexList();
454  vector<NodeSharedPtr> ns2 = e2->GetVertexList();
455  if (ns1[0] == ns2[0] || ns1[0] == ns2[1] || ns1[0] == ns2[2] ||
456  ns1[1] == ns2[0] || ns1[1] == ns2[1] || ns1[1] == ns2[2] ||
457  ns1[2] == ns2[0] || ns1[2] == ns2[1] || ns1[2] == ns2[2])
458  {
459  return false;
460  }
461 
462  Array<OneD, NekDouble> N1(3), N2(3);
463  NekDouble d1, d2;
464 
465  N1[0] = (ns1[1]->m_y - ns1[0]->m_y) * (ns1[2]->m_z - ns1[0]->m_z) -
466  (ns1[2]->m_y - ns1[0]->m_y) * (ns1[1]->m_z - ns1[0]->m_z);
467  N1[1] = -1.0 * ((ns1[1]->m_x - ns1[0]->m_x) * (ns1[2]->m_z - ns1[0]->m_z) -
468  (ns1[2]->m_x - ns1[0]->m_x) * (ns1[1]->m_z - ns1[0]->m_z));
469  N1[2] = (ns1[1]->m_x - ns1[0]->m_x) * (ns1[2]->m_y - ns1[0]->m_y) -
470  (ns1[2]->m_x - ns1[0]->m_x) * (ns1[1]->m_y - ns1[0]->m_y);
471 
472  N2[0] = (ns2[1]->m_y - ns2[0]->m_y) * (ns2[2]->m_z - ns2[0]->m_z) -
473  (ns2[2]->m_y - ns2[0]->m_y) * (ns2[1]->m_z - ns2[0]->m_z);
474  N2[1] = -1.0 * ((ns2[1]->m_x - ns2[0]->m_x) * (ns2[2]->m_z - ns2[0]->m_z) -
475  (ns2[2]->m_x - ns2[0]->m_x) * (ns2[1]->m_z - ns2[0]->m_z));
476  N2[2] = (ns2[1]->m_x - ns2[0]->m_x) * (ns2[2]->m_y - ns2[0]->m_y) -
477  (ns2[2]->m_x - ns2[0]->m_x) * (ns2[1]->m_y - ns2[0]->m_y);
478 
479  d1 = -1.0 *
480  (N1[0] * ns1[0]->m_x + N1[1] * ns1[0]->m_y + N1[2] * ns1[0]->m_z);
481  d2 = -1.0 *
482  (N2[0] * ns2[0]->m_x + N2[1] * ns2[0]->m_y + N2[2] * ns2[0]->m_z);
483 
484  Array<OneD, NekDouble> dv1(3), dv2(3);
485 
486  dv1[0] =
487  N2[0] * ns1[0]->m_x + N2[1] * ns1[0]->m_y + N2[2] * ns1[0]->m_z + d2;
488  dv1[1] =
489  N2[0] * ns1[1]->m_x + N2[1] * ns1[1]->m_y + N2[2] * ns1[1]->m_z + d2;
490  dv1[2] =
491  N2[0] * ns1[2]->m_x + N2[1] * ns1[2]->m_y + N2[2] * ns1[2]->m_z + d2;
492 
493  dv2[0] =
494  N1[0] * ns2[0]->m_x + N1[1] * ns2[0]->m_y + N1[2] * ns2[0]->m_z + d1;
495  dv2[1] =
496  N1[0] * ns2[1]->m_x + N1[1] * ns2[1]->m_y + N1[2] * ns2[1]->m_z + d1;
497  dv2[2] =
498  N1[0] * ns2[2]->m_x + N1[1] * ns2[2]->m_y + N1[2] * ns2[2]->m_z + d1;
499 
500  if (sign(dv1[0], dv1[1]) && sign(dv1[1], dv1[2]))
501  {
502  return false;
503  }
504  if (sign(dv2[0], dv2[1]) && sign(dv2[1], dv2[2]))
505  {
506  return false;
507  }
508 
509  Array<OneD, NekDouble> D(3);
510  D[0] = N1[1] * N2[2] - N1[2] * N2[1];
511  D[1] = -1.0 * (N1[0] * N2[2] - N1[2] * N2[0]);
512  D[2] = N1[0] * N2[1] - N1[0] * N2[1];
513 
514  int base1 = 0, base2 = 0;
515  if (!sign(dv2[0], dv2[1]) && sign(dv2[1], dv2[2]))
516  {
517  base2 = 0;
518  }
519  else if (!sign(dv2[1], dv2[2]) && sign(dv2[2], dv2[0]))
520  {
521  base2 = 1;
522  }
523  else if (!sign(dv2[2], dv2[0]) && sign(dv2[0], dv2[1]))
524  {
525  base2 = 2;
526  }
527  else
528  {
529  cout << "base not set" << endl;
530  }
531 
532  if (!sign(dv1[0], dv1[1]) && sign(dv1[1], dv1[2]))
533  {
534  base1 = 0;
535  }
536  else if (!sign(dv1[1], dv1[2]) && sign(dv1[2], dv1[0]))
537  {
538  base1 = 1;
539  }
540  else if (!sign(dv1[2], dv1[0]) && sign(dv1[0], dv1[1]))
541  {
542  base1 = 2;
543  }
544  else
545  {
546  cout << "base not set" << endl;
547  }
548 
549  Array<OneD, NekDouble> p1(3), p2(3);
550 
551  p1[0] = D[0] * ns1[0]->m_x + D[1] * ns1[0]->m_y + D[2] * ns1[0]->m_z;
552  p1[1] = D[0] * ns1[1]->m_x + D[1] * ns1[1]->m_y + D[2] * ns1[1]->m_z;
553  p1[2] = D[0] * ns1[2]->m_x + D[1] * ns1[2]->m_y + D[2] * ns1[2]->m_z;
554 
555  p2[0] = D[0] * ns2[0]->m_x + D[1] * ns2[0]->m_y + D[2] * ns2[0]->m_z;
556  p2[1] = D[0] * ns2[1]->m_x + D[1] * ns2[1]->m_y + D[2] * ns2[1]->m_z;
557  p2[2] = D[0] * ns2[2]->m_x + D[1] * ns2[2]->m_y + D[2] * ns2[2]->m_z;
558 
559  NekDouble t11, t12, t21, t22;
560  int o1 = 0, o2 = 0;
561  if (base1 == 0)
562  {
563  o1 = 1;
564  o2 = 2;
565  }
566  else if (base1 == 1)
567  {
568  o1 = 2;
569  o2 = 0;
570  }
571  else if (base1 == 2)
572  {
573  o1 = 0;
574  o2 = 1;
575  }
576 
577  t11 = p1[o1] + (p1[base1] - p1[o1]) * dv1[o1] / (dv1[o1] - dv1[base1]);
578  t12 = p1[o2] + (p1[base1] - p1[o2]) * dv1[o2] / (dv1[o2] - dv1[base1]);
579 
580  if (base2 == 0)
581  {
582  o1 = 1;
583  o2 = 2;
584  }
585  else if (base2 == 1)
586  {
587  o1 = 2;
588  o2 = 0;
589  }
590  else if (base2 == 2)
591  {
592  o1 = 0;
593  o2 = 1;
594  }
595 
596  t21 = p2[o1] + (p2[base2] - p2[o1]) * dv2[o1] / (dv2[o1] - dv2[base2]);
597  t22 = p2[o2] + (p2[base2] - p2[o2]) * dv2[o2] / (dv2[o2] - dv2[base2]);
598 
599  if (t11 > t12)
600  {
601  swap(t11, t12);
602  }
603  if (t21 > t22)
604  {
605  swap(t21, t22);
606  }
607 
608  if (t21 < t11)
609  {
610  swap(t11, t21);
611  swap(t12, t22);
612  }
613 
614  if (!sign(t21 - t11, t22 - t11) || !sign(t21 - t12, t22 - t12))
615  {
616  return true;
617  }
618 
619  return false;
620 }
bool sign(NekDouble a, NekDouble b)
Definition: BLMesh.cpp:352
double NekDouble
NekDouble Nektar::NekMeshUtils::BLMesh::Visability ( std::vector< ElementSharedPtr tris,
Array< OneD, NekDouble N 
)
private

Definition at line 830 of file BLMesh.cpp.

832 {
833  NekDouble mn = numeric_limits<double>::max();
834 
835  for (int i = 0; i < tris.size(); i++)
836  {
837  Array<OneD, NekDouble> tmp = tris[i]->Normal(true);
838  NekDouble dt = tmp[0] * N[0] + tmp[1] * N[1] + tmp[2] * N[2];
839  mn = min(mn, dt);
840  }
841  return mn;
842 }
double NekDouble

Friends And Related Function Documentation

friend class MemoryManager< BLMesh >
friend

Definition at line 52 of file BLMesh.h.

Member Data Documentation

NekDouble Nektar::NekMeshUtils::BLMesh::m_bl
private

thickness of the boundary layer

Definition at line 125 of file BLMesh.h.

std::map<NodeSharedPtr, blInfoSharedPtr> Nektar::NekMeshUtils::BLMesh::m_blData
private

data structure used to store and develop bl information

Definition at line 133 of file BLMesh.h.

std::vector<unsigned int> Nektar::NekMeshUtils::BLMesh::m_blsurfs
private

List of surfaces onto which boundary layers are placed.

Definition at line 123 of file BLMesh.h.

Referenced by GetBLSurfs().

NekMatrix<NekDouble> Nektar::NekMeshUtils::BLMesh::m_deriv[3]
private

Definition at line 138 of file BLMesh.h.

int Nektar::NekMeshUtils::BLMesh::m_id
private

Definition at line 128 of file BLMesh.h.

int Nektar::NekMeshUtils::BLMesh::m_layer
private

Definition at line 127 of file BLMesh.h.

Array<OneD, NekDouble> Nektar::NekMeshUtils::BLMesh::m_layerT
private

Definition at line 129 of file BLMesh.h.

MeshSharedPtr Nektar::NekMeshUtils::BLMesh::m_mesh
private

mesh object containing surface mesh

Definition at line 121 of file BLMesh.h.

std::map<NodeSharedPtr, std::vector<blInfoSharedPtr> > Nektar::NekMeshUtils::BLMesh::m_nToNInfo
private

Definition at line 135 of file BLMesh.h.

std::map<ElementSharedPtr, ElementSharedPtr> Nektar::NekMeshUtils::BLMesh::m_priToTri
private

Definition at line 136 of file BLMesh.h.

NekDouble Nektar::NekMeshUtils::BLMesh::m_prog
private

Definition at line 126 of file BLMesh.h.

std::vector<ElementSharedPtr> Nektar::NekMeshUtils::BLMesh::m_psuedoSurface
private

Definition at line 137 of file BLMesh.h.

Referenced by GetPseudoSurface().

std::vector<unsigned int> Nektar::NekMeshUtils::BLMesh::m_symSurfs
private

list of surfaces to be remeshed due to the boundary layer

Definition at line 131 of file BLMesh.h.

Referenced by GetSymSurfs().