From 83780665e68a43b840ecb6ab04ae2a3e83ae67d8 Mon Sep 17 00:00:00 2001 From: Oliver Gloth Date: Thu, 15 Jul 2010 09:02:26 +0200 Subject: [PATCH] using FaceFinder now -- no neighbour search anymore --- src/surfaceprojection.cpp | 113 +++++++++++++++++++++++++--------------------- src/surfaceprojection.h | 42 ++++++++++------- 2 files changed, 86 insertions(+), 69 deletions(-) diff --git a/src/surfaceprojection.cpp b/src/surfaceprojection.cpp index afad86c..423f913 100644 --- a/src/surfaceprojection.cpp +++ b/src/surfaceprojection.cpp @@ -34,7 +34,6 @@ SurfaceProjection::SurfaceProjection(int bc) : SurfaceAlgorithm() this->setGrid(m_BGrid); m_BC = bc; getSet("surface meshing", "correct curvature (experimental)", false, m_correctCurvature); - getSet("surface meshing", "number ofneighbour searches", 2, m_NumNeighSearches); m_CritDistance = 0.1; } @@ -53,71 +52,46 @@ void SurfaceProjection::searchNewTriangle(vec3_t xp, vtkIdType &id_tri, vec3_t & x_proj = vec3_t(1e99, 1e99, 1e99); r_proj = vec3_t(0, 0, 0); vec3_t xi, ri; + double d_min = 1e99; + bool x_proj_set = false; on_triangle = false; - double d_min = 1e99; - bool tri_found = false; - if (id_tri >= 0) { - QSet tris; - tris.insert(id_tri); - for (int i_search = 0; i_search < m_NumNeighSearches; ++i_search) { - //for (int i_search = 0; i_search < 100; ++i_search) { - QSet candidates; - foreach (int i_tri, tris) { - for (int k = 0; k < m_BPart.c2cGSize(i_tri); ++k) { - int new_cand = m_BPart.c2cGG(i_tri,k); - if (!tris.contains(new_cand) && new_cand >= 0) { - candidates.insert(m_BPart.c2cGG(i_tri,k)); - } - } - } - foreach (int i_cand, candidates) { - if (i_cand < 0) EG_BUG; - if (i_cand > m_BGrid->GetNumberOfCells()) EG_BUG; - int side; - double d; - if(m_Triangles[i_cand].projectOnTriangle(xp, xi, ri, d, side, m_RestrictToTriangle)) { - x_proj = xi; - r_proj = ri; - on_triangle = true; - id_tri = i_cand; - ++Nhalf; - return; - } else if (d < m_CritDistance*m_Triangles[i_cand].m_smallest_length) { - if (d < d_min) { - x_proj = xi; - r_proj = ri; - id_tri = i_cand; - tri_found = true; - } - } - } - tris += candidates; - if (tri_found) { - return; - } + QVector close_faces; + m_FaceFinder.getCloseFaces(xp, close_faces); + foreach (vtkIdType id_triangle, close_faces) { + Triangle T = m_Triangles[id_triangle]; + double d; + int side; + bool intersects = T.projectOnTriangle(xp, xi, ri, d, side, m_RestrictToTriangle); + if (d >= 1e99) { + EG_BUG; + } + if (d < d_min) { + x_proj = xi; + x_proj_set = true; + r_proj = ri; + d_min = d; + id_tri = id_triangle; + on_triangle = intersects; } } - if (tri_found) { - ++Nhalf; - } else { - // full search :-( - ++Nfull; - bool x_proj_set = false; - d_min = 1e99; - for (int i_triangles = 0; i_triangles < m_Triangles.size(); ++i_triangles) { - Triangle T = m_Triangles[i_triangles]; + if (!x_proj_set) { // should never happen + for (vtkIdType id_triangle = 0; id_triangle < m_BGrid->GetNumberOfCells(); ++id_triangle) { + Triangle T = m_Triangles[id_triangle]; double d; int side; bool intersects = T.projectOnTriangle(xp, xi, ri, d, side, m_RestrictToTriangle); if (d >= 1e99) { EG_BUG; } + if (d < 0) { + EG_BUG; + } if (d < d_min) { x_proj = xi; x_proj_set = true; r_proj = ri; d_min = d; - id_tri = i_triangles; + id_tri = id_triangle; on_triangle = intersects; } } @@ -189,6 +163,7 @@ void SurfaceProjection::updateBackgroundGridInfo() } m_BPart.setGrid(m_BGrid); m_BPart.setAllCells(); + computeSurfaceCurvature(); } @@ -288,3 +263,37 @@ vec3_t SurfaceProjection::correctCurvature(int, vec3_t g_M) { return g_M; } + +void SurfaceProjection::computeSurfaceCurvature() +{ + m_Radius.fill(1e99, m_BGrid->GetNumberOfCells()); + for (vtkIdType id_cell = 0; id_cell < m_BGrid->GetNumberOfCells(); ++id_cell) { + vtkIdType N_pts, *pts; + m_BGrid->GetCellPoints(id_cell, N_pts, pts); + vec3_t x[N_pts+1]; + vec3_t n[N_pts+1]; + for (int i = 0; i < N_pts; ++i) { + m_BGrid->GetPoint(pts[i], x[i].data()); + n[i] = m_NodeNormals[pts[i]]; + } + x[N_pts] = x[0]; + n[N_pts] = n[0]; + for (int i = 0; i < N_pts; ++i) { + double alpha = max(1e-3,acos(n[i]*n[i+1])); + double a = (x[i]-x[i+1]).abs(); + m_Radius[id_cell] = min(m_Radius[id_cell], a/alpha); + } + } +} + +double SurfaceProjection::getRadius(vtkIdType id_node) +{ + vec3_t x; + m_FGrid->GetPoint(id_node, x.data()); + projectRestricted(x, id_node); + vtkIdType id_tri = getProjTriangle(id_node); + if (id_tri == -1) { + EG_BUG; + } + return m_Radius[id_tri]; +} diff --git a/src/surfaceprojection.h b/src/surfaceprojection.h index c190c84..e15eace 100644 --- a/src/surfaceprojection.h +++ b/src/surfaceprojection.h @@ -32,6 +32,7 @@ class SurfaceProjection; #include "surfaceoperation.h" #include "surfacealgorithm.h" #include "triangle.h" +#include "facefinder.h" class SurfaceProjection : public SurfaceAlgorithm { @@ -53,21 +54,22 @@ private: // methods protected: // attributes - vtkUnstructuredGrid* m_BGrid; ///< the background grid defining the geometry - MeshPartition m_BPart; - vtkUnstructuredGrid* m_FGrid; ///< the foreground grid to project - QVector m_EdgeLength; - QVector m_Cells; - QVector m_Nodes; - QVector m_NodeNormals; ///< The surface normal at each node of m_BGrid - QVector m_Triangles; ///< All triangles of m_BGrid. One for each triangle cell of m_BGrid. - QVector > m_N2N; - bool m_correctCurvature; ///< Should correctCurvature() be used? - int m_NumNeighSearches; - int m_BC; - bool m_RestrictToTriangle; - double m_CritDistance; + vtkUnstructuredGrid* m_BGrid; ///< the background grid defining the geometry + MeshPartition m_BPart; + vtkUnstructuredGrid* m_FGrid; ///< the foreground grid to project + QVector m_EdgeLength; + QVector m_Cells; + QVector m_Nodes; + QVector m_NodeNormals; ///< The surface normal at each node of m_BGrid + QVector m_Triangles; ///< All triangles of m_BGrid. One for each triangle cell of m_BGrid. + QVector m_Radius; ///< Surface radius for mesh resolution. + QVector > m_N2N; + bool m_correctCurvature; ///< Should correctCurvature() be used? + int m_BC; + bool m_RestrictToTriangle; + double m_CritDistance; QMap m_Pindex; + FaceFinder m_FaceFinder; protected: // static attributes @@ -81,6 +83,7 @@ protected: // methods void searchNewTriangle(vec3_t xp, vtkIdType &id_tri, vec3_t &x_proj, vec3_t &r_proj, bool &on_triangle); vtkIdType getProjTriangle(vtkIdType id_node); void setProjTriangle(vtkIdType id_node, vtkIdType proj_triangle); + void computeSurfaceCurvature(); public: // methods @@ -97,9 +100,10 @@ public: // methods virtual vec3_t projectFree(vec3_t x, vtkIdType id_node = -1); vtkUnstructuredGrid* getBGrid() { return m_BGrid; } + double getRadius(vtkIdType id_node); void setCorrectCurvature(bool b) { m_correctCurvature = b; } - bool getCorrectCurvature() { return m_correctCurvature; } + bool getCorrectCurvature() { return m_correctCurvature; } public: // static methods @@ -116,12 +120,16 @@ void SurfaceProjection::setBackgroundGrid(vtkUnstructuredGrid* grid, const C& ce setForegroundGrid(grid); for (int i_cells = 0; i_cells < cells.size(); ++i_cells) { vtkIdType id_cell = cells[i_cells]; + } + for (vtkIdType id_cell = 0; id_cell < m_BGrid->GetNumberOfCells(); ++id_cell) { vtkIdType N_pts, *pts; - grid->GetCellPoints(id_cell, N_pts, pts); + m_BGrid->GetCellPoints(id_cell, N_pts, pts); for (int i = 0; i < N_pts; ++i) { - setProjTriangle(pts[i], i_cells); + setProjTriangle(pts[i], id_cell); } } + m_FaceFinder.setMaxNumFaces(100); + m_FaceFinder.setGrid(m_BGrid); } template -- 2.11.4.GIT