From e32512deb102a8fda319c651a50de31b17eacfb5 Mon Sep 17 00:00:00 2001 From: Oliver Gloth Date: Mon, 6 Aug 2012 20:10:42 +0200 Subject: [PATCH] fixed various problems with BRL projector and edge swapping --- src/libengrid/brlcadinterface.cpp | 16 ++++++++++------ src/libengrid/brlcadinterface.h | 3 ++- src/libengrid/brlcadprojection.cpp | 14 ++++++++++---- src/libengrid/brlcadprojection.h | 1 + src/libengrid/insertpoints.cpp | 1 - src/libengrid/meshpartition.cpp | 10 +++++++--- src/libengrid/removepoints.cpp | 3 --- src/libengrid/swaptriangles.cpp | 9 +++++++-- 8 files changed, 37 insertions(+), 20 deletions(-) diff --git a/src/libengrid/brlcadinterface.cpp b/src/libengrid/brlcadinterface.cpp index 89367d6..66b4f9f 100644 --- a/src/libengrid/brlcadinterface.cpp +++ b/src/libengrid/brlcadinterface.cpp @@ -75,7 +75,7 @@ int BrlCadInterface::hit(application *ap, struct partition *PartHeadp, seg *segs RT_HIT_NORMAL(inormal, hitp, stp, &(ap->a_ray), pp->pt_inflip); RT_CURVATURE(&cur, hitp, pp->pt_inflip, stp); - curv = max(cur.crv_c1, cur.crv_c2); + curv = max(fabs(cur.crv_c1), fabs(cur.crv_c2)); m_InRadius = 1.0/max(1e-10, curv); m_XIn[0] = pt[0]; @@ -90,7 +90,7 @@ int BrlCadInterface::hit(application *ap, struct partition *PartHeadp, seg *segs VJOIN1(pt, ap->a_ray.r_pt, hitp->hit_dist, ap->a_ray.r_dir); RT_HIT_NORMAL( onormal, hitp, stp, &(ap->a_ray), pp->pt_outflip ); RT_CURVATURE(&cur, hitp, pp->pt_inflip, stp); - curv = max(cur.crv_c1, cur.crv_c2); + curv = max(fabs(cur.crv_c1), fabs(cur.crv_c2)); m_OutRadius = 1.0/max(1e-10, curv); m_XOut[0] = pt[0]; @@ -171,12 +171,16 @@ BrlCadInterface::HitType BrlCadInterface::shootRay(vec3_t x, vec3_t v, vec3_t &x return hit_type; } -bool BrlCadInterface::isInside(vec3_t x) +BrlCadInterface::PositionType BrlCadInterface::position(vec3_t x, vec3_t n) { vec3_t x_hit, n_hit; double r_hit; - if (shootRay(x, vec3_t(1,0,0), x_hit, n_hit, r_hit) == HitOut) { - return true; + HitType hit_type = shootRay(x, vec3_t(1,0,0), x_hit, n_hit, r_hit); + if (hit_type == HitOut) { + return Inside; } - return false; + if (hit_type == HitIn) { + return Outside; + } + return Surface; } diff --git a/src/libengrid/brlcadinterface.h b/src/libengrid/brlcadinterface.h index 9053fb4..15e660c 100644 --- a/src/libengrid/brlcadinterface.h +++ b/src/libengrid/brlcadinterface.h @@ -37,6 +37,7 @@ class BrlCadInterface public: // data types enum HitType { Miss, HitIn, HitOut }; + enum PositionType { Inside, Outside, Surface }; private: // attributes @@ -67,7 +68,7 @@ protected: // methods HitType shootRay(vec3_t x, vec3_t v, vec3_t &x_hit, vec3_t &n_hit, double &r); void setupBrlCad(QString file_name, QString object_name); - bool isInside(vec3_t x); + PositionType position(vec3_t x, vec3_t n); public: diff --git a/src/libengrid/brlcadprojection.cpp b/src/libengrid/brlcadprojection.cpp index e76d7a4..207f3f4 100644 --- a/src/libengrid/brlcadprojection.cpp +++ b/src/libengrid/brlcadprojection.cpp @@ -26,6 +26,7 @@ BrlCadProjection::BrlCadProjection(QString file_name, QString object_name) { setupBrlCad(file_name, object_name); + m_ForceRay = false; } BrlCadProjection::~BrlCadProjection() @@ -36,9 +37,6 @@ BrlCadProjection::~BrlCadProjection() vec3_t BrlCadProjection::project(vec3_t x, vtkIdType id_node, bool, vec3_t v) { vec3_t n = v; - if (id_node == 0) { - cout << "break"<< endl; - } if (n.abs() < 1e-3) { if (id_node == -1) { EG_BUG; @@ -49,6 +47,8 @@ vec3_t BrlCadProjection::project(vec3_t x, vtkIdType id_node, bool, vec3_t v) EG_BUG; } if (!checkVector(n)) { + cout << "vector defect (id_node=" << id_node << endl; + //return x; EG_BUG; } @@ -57,7 +57,11 @@ vec3_t BrlCadProjection::project(vec3_t x, vtkIdType id_node, bool, vec3_t v) vec3_t x_hit, n_hit; double r_hit; - if (isInside(x)) { + PositionType pos_type = position(x, n); + if (pos_type == Surface && !m_ForceRay) { + return x; + } + if (pos_type == Outside) { n *= -1; } @@ -74,6 +78,8 @@ double BrlCadProjection::getRadius(vtkIdType id_node) { vec3_t x; m_FGrid->GetPoint(id_node, x.data()); + m_ForceRay = true; project(x, id_node); + m_ForceRay = false; return m_LastRadius; } diff --git a/src/libengrid/brlcadprojection.h b/src/libengrid/brlcadprojection.h index a9d6963..b60a7c6 100644 --- a/src/libengrid/brlcadprojection.h +++ b/src/libengrid/brlcadprojection.h @@ -34,6 +34,7 @@ class BrlCadProjection : public SurfaceProjection, public BrlCadInterface vec3_t m_LastNormal; double m_LastRadius; + bool m_ForceRay; public: diff --git a/src/libengrid/insertpoints.cpp b/src/libengrid/insertpoints.cpp index e02ad08..857db38 100644 --- a/src/libengrid/insertpoints.cpp +++ b/src/libengrid/insertpoints.cpp @@ -36,7 +36,6 @@ InsertPoints::InsertPoints() : SurfaceOperation() void InsertPoints::operate() { - return; int N1 = m_Grid->GetNumberOfPoints(); insertPoints(); int N2 = m_Grid->GetNumberOfPoints(); diff --git a/src/libengrid/meshpartition.cpp b/src/libengrid/meshpartition.cpp index 5e6185e..b8a8066 100644 --- a/src/libengrid/meshpartition.cpp +++ b/src/libengrid/meshpartition.cpp @@ -340,7 +340,7 @@ bool MeshPartition::hasBC(vtkIdType id_node, int bc) bool found = false; for (int j = 0; j < n2bcGSize(id_node); ++j) { if (n2bcG(id_node, j) == bc) { - found == true; + found = true; break; } } @@ -355,7 +355,6 @@ vtkIdType MeshPartition::getVolumeCell(vtkIdType id_face) vec3_t MeshPartition::globalNormal(vtkIdType id_node) { vec3_t normal(0,0,0); - int murks = id_node; for (int i = 0; i < n2cGSize(id_node); ++i) { vtkIdType id_cell = n2cGG(id_node, i); if (isSurface(id_cell, m_Grid)) { @@ -382,10 +381,15 @@ vec3_t MeshPartition::globalNormal(vtkIdType id_node) double alpha = GeometryTools::angle(u, v); vec3_t n = u.cross(v); n.normalise(); - normal += alpha*n; + if (checkVector(n)) { + normal -= alpha*n; + } } } normal.normalise(); + if (!checkVector(normal)) { + EG_BUG; + } return normal; } diff --git a/src/libengrid/removepoints.cpp b/src/libengrid/removepoints.cpp index 2dea365..61d101f 100644 --- a/src/libengrid/removepoints.cpp +++ b/src/libengrid/removepoints.cpp @@ -91,9 +91,6 @@ void RemovePoints::fixNodes(const QVector &fixnodes) void RemovePoints::operate() { - return; - - if (m_Fixed.size() != m_Grid->GetNumberOfPoints()) { m_Fixed.fill(false, m_Grid->GetNumberOfPoints()); } diff --git a/src/libengrid/swaptriangles.cpp b/src/libengrid/swaptriangles.cpp index 38ef206..c2c92ce 100644 --- a/src/libengrid/swaptriangles.cpp +++ b/src/libengrid/swaptriangles.cpp @@ -108,6 +108,7 @@ bool SwapTriangles::isEdge(vtkIdType id_node1, vtkIdType id_node2) void SwapTriangles::computeSurfaceErrors(const QVector &x, int bc, double &err1, double &err2) { + static int counter = 0; using namespace GeometryTools; err1 = 0; err2 = 0; @@ -138,6 +139,7 @@ void SwapTriangles::computeSurfaceErrors(const QVector &x, int bc, doubl err1 = (xe1 - xe1_proj).abs()/L; err2 = (xe2 - xe2_proj).abs()/L; + ++counter; } int SwapTriangles::swap() @@ -147,19 +149,22 @@ int SwapTriangles::swap() EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code"); QVector marked(m_Grid->GetNumberOfCells(), false); for (int i = 0; i < m_Part.getNumberOfCells(); ++i) { - //m_Timer << " cell " << i+1 << "/" << m_Part.getNumberOfCells() << Timer::endl; vtkIdType id_cell = m_Part.globalCell(i); if (!m_BoundaryCodes.contains(cell_code->GetValue(id_cell)) && m_Grid->GetCellType(id_cell) == VTK_TRIANGLE) { //if it is a selected triangle if (!marked[id_cell] && !m_Swapped[id_cell]) { for (int j = 0; j < 3; ++j) { bool swap = false; stencil_t S = getStencil(id_cell, j); + if (S.id_cell.size() == 2) { + if ((S.id_cell[1] == 2171 && S.id_cell[0] == 2110) || (S.id_cell[0] == 2171 && S.id_cell[1] == 2110)) { + cout << "break" << endl; + } + } if(S.id_cell.size() == 2 && S.sameBC) { if (S.type_cell[1] == VTK_TRIANGLE) { if(!isEdge(S.id_node[0], S.id_node[1]) ) { if (!marked[S.id_cell[1]] && !m_Swapped[S.id_cell[1]]) { QVector x3(4); - vec3_t x3_0(0,0,0); vec2_t x[4]; m_Grid->GetPoint(S.id_node[0], x3[0].data()); -- 2.11.4.GIT