From 3c446eaae24e41c04c9d24dd828a2239b094c333 Mon Sep 17 00:00:00 2001 From: Oliver Gloth Date: Mon, 24 Sep 2012 09:18:49 +0200 Subject: [PATCH] new vertex type for feature corners + geometric recogintion --- src/libengrid/egvtkobject.h | 8 +- src/libengrid/laplacesmoother.cpp | 153 ++++++++-------- src/libengrid/laplacesmoother.h | 3 + src/libengrid/surfaceoperation.cpp | 356 +++++++++++++++++++------------------ src/libengrid/surfaceoperation.h | 17 +- src/libengrid/swaptriangles.cpp | 25 +-- src/libengrid/utilities.cpp | 18 +- 7 files changed, 301 insertions(+), 279 deletions(-) diff --git a/src/libengrid/egvtkobject.h b/src/libengrid/egvtkobject.h index b0c5209..9ac6db7 100644 --- a/src/libengrid/egvtkobject.h +++ b/src/libengrid/egvtkobject.h @@ -45,10 +45,10 @@ class BezierTriangle; #include #define EG_SIMPLE_VERTEX 0 -#define EG_FIXED_VERTEX 1 -#define EG_FEATURE_EDGE_VERTEX 2 -#define EG_BOUNDARY_EDGE_VERTEX 3 -#define EG_DETECTED_EDGE_VERTEX 4 +#define EG_FEATURE_EDGE_VERTEX 1 +#define EG_BOUNDARY_EDGE_VERTEX 2 +#define EG_FEATURE_CORNER_VERTEX 3 +#define EG_FIXED_VERTEX 4 class EgVtkObject { diff --git a/src/libengrid/laplacesmoother.cpp b/src/libengrid/laplacesmoother.cpp index 98b2132..0678df1 100644 --- a/src/libengrid/laplacesmoother.cpp +++ b/src/libengrid/laplacesmoother.cpp @@ -1,4 +1,4 @@ -// +// // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // + + // + This file is part of enGrid. + @@ -26,8 +26,8 @@ #include #include "guimainwindow.h" -#include "globalnodegraphinterface.h" -#include "surfacenodemovementcheck.h" +#include "localnodegraphinterface.h" +#include "checkerboardgraphiterator.h" using namespace GeometryTools; @@ -55,16 +55,24 @@ bool LaplaceSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new) vec3_t x_old; m_Grid->GetPoint(id_node, x_old.data()); m_Grid->GetPoints()->SetPoint(id_node, x_new.data()); + + EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type"); bool move = true; - if(m_NoCheck) return move; + //if(m_NoCheck || node_type->GetValue(id_node) != EG_SIMPLE_VERTEX) { + if(m_NoCheck) { + return move; + } + + //move = m_Check.checkNode(id_node, x_new); + // compute the extrusion vector to compute the tetrahedrons for volume checking // start with an average of all adjacent cell normals and count the number of // adjacent boundary codes (for one boundary code an alternative vector can be // computed with the help of a SurfaceProjection // - vec3_t n(0,0,0); + vec3_t n1(0,0,0); QVector cell_normals(m_Part.n2cGSize(id_node)); QSet bcs; EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code"); @@ -80,25 +88,23 @@ bool LaplaceSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new) for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) { double A = fabs(GeometryTools::cellVA(m_Grid, m_Part.n2cGG(id_node, i))); if (A > 0.01*A_max) { - n += cell_normals[i]; + n1 += cell_normals[i]; ++N; } } - /* - SurfaceProjection* proj = NULL; - if (bcs.size() == 1) { - proj = GuiMainWindow::pointer()->getSurfProj(*bcs.begin()); - if (proj) { - proj->project(x_new, id_node, false, n); - n = proj->lastProjNormal(); - } - } - */ - if (N == 0) { - //EG_BUG; + n1.normalise(); + if (!checkVector(n1)) { move = false; } else { - n.normalise(); + vec3_t n2 = n1; + SurfaceProjection* proj = NULL; + if (bcs.size() == 1) { + proj = GuiMainWindow::pointer()->getSurfProj(*bcs.begin()); + if (proj) { + proj->project(x_new, id_node, false, n1); + n2 = proj->lastProjNormal(); + } + } double L_max = 0;// distance to the furthest neighbour node of id_node for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) { vec3_t xn; @@ -107,15 +113,10 @@ bool LaplaceSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new) L_max = max(L, L_max); } - vec3_t x_summit; - if(m_CorrectCurvature) { - // better for mesher with interpolation - x_summit = x_new + L_max*n; - } - else { - // better for mesher without interpolation - x_summit = x_old + L_max*n; - } + vec3_t x_summit_old1 = x_old + L_max*n1; + vec3_t x_summit_new1 = x_new + L_max*n1; + vec3_t x_summit_old2 = x_old + L_max*n2; + vec3_t x_summit_new2 = x_new + L_max*n2; for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) { vec3_t x[3]; @@ -129,33 +130,28 @@ bool LaplaceSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new) for (int j = 0; j < N_pts; ++j) { m_Grid->GetPoint(pts[j], x[j].data()); } - if (GeometryTools::tetraVol(x[0], x[1], x[2], x_summit, false) <= 0) { + if (GeometryTools::tetraVol(x[0], x[1], x[2], x_summit_old1, false) <= 0) { move = false; - break; } - } - } - /* - if (move) { - for (int i = 0; i < cell_normals.size(); ++i) { - for (int j = 0; j < cell_normals.size(); ++j) { - if (cell_normals[i]*cell_normals[j] < -100.0) {//Why -100.0? - saveGrid(m_Grid, "after_move2"); - qWarning()<<"Cannot move point "<GetPoints()->SetPoint(id_node, x_old.data()); - saveGrid(m_Grid, "before_move2"); - EG_BUG; - break; - } + if (GeometryTools::tetraVol(x[0], x[1], x[2], x_summit_new1, false) <= 0) { + move = false; + } + if (GeometryTools::tetraVol(x[0], x[1], x[2], x_summit_old2, false) <= 0) { + //move = false; + } + if (GeometryTools::tetraVol(x[0], x[1], x[2], x_summit_new2, false) <= 0) { + //move = false; + } + if (!move) { + break; } } } - */ + if (!move) { m_Grid->GetPoints()->SetPoint(id_node, x_old.data()); } - return move; + return move; } void LaplaceSmoother::featureCorrection(vtkIdType id_node, SurfaceProjection* proj, vec3_t &x_new) @@ -188,7 +184,7 @@ void LaplaceSmoother::featureCorrection(vtkIdType id_node, SurfaceProjection* pr double L1 = 0; double L2 = m_FeatureMagic*cl->GetValue(id_node); for (int i = 0; i < 30; ++i) { - x_new = x1 - 0.5*(L1 + L2)*magic_vector; + x_new = x1 + 0.5*(L1 + L2)*magic_vector; x_new = proj->project(x_new, id_node, m_CorrectCurvature, n); double displacement = fabs((x_new - x1)*n); if (displacement > 0.01*cl->GetValue(id_node)) { @@ -207,7 +203,7 @@ void LaplaceSmoother::featureCorrection(vtkIdType id_node, SurfaceProjection* pr L1 = 0.5*(L1 + L2); } } - x_new = x1 - L1*magic_vector; + x_new = x1 + L1*magic_vector; x_new = proj->project(x_new, id_node, m_CorrectCurvature, n); } else { @@ -247,7 +243,7 @@ bool LaplaceSmoother::moveNode(vtkIdType id_node, vec3_t &Dx) */ // TESTING END - for (int i_relaxation = 0; i_relaxation < 1; ++i_relaxation) { + for (int i_relaxation = 0; i_relaxation < 10; ++i_relaxation) { vec3_t x_new = x_old + Dx; if (m_UseProjection) { int i_nodes = m_Part.localNode(id_node); @@ -327,7 +323,7 @@ void LaplaceSmoother::operate() GuiMainWindow::pointer()->getSurfProj(bc)->setForegroundGrid(m_Grid); } } - UpdatePotentialSnapPoints(false, false); + updateNodeInfo(); EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code"); EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type" ); EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired"); @@ -366,45 +362,58 @@ void LaplaceSmoother::operate() for (int i_iter = 0; i_iter < m_NumberOfIterations; ++i_iter) { m_Success = true; computeNormals(); - for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) { - vtkIdType id_node = nodes[i_nodes]; - if (!m_Fixed[id_node] && !blocked[i_nodes]) { + //m_Check.setGrid(m_Grid); + LocalNodeGraphInterface graph; + graph.setMeshPartition(&m_Part); + CheckerBoardGraphIterator iter; + iter.setGraph(&graph); + + //for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) { + for (iter = 0; iter < nodes.size(); ++iter) { + if (iter.updateRequired()) { + //m_Check.setGrid(m_Grid); + } + vtkIdType id_node = nodes[*iter]; + if (!m_Fixed[id_node] && !blocked[*iter]) { if (smooth_node[id_node] && node_type->GetValue(id_node) != EG_FIXED_VERTEX) { if (node_type->GetValue(id_node) != EG_FIXED_VERTEX) { QVector snap_points = getPotentialSnapPoints(id_node); vec3_t n(0,0,0); + vec3_t x_old; + m_Grid->GetPoint(id_node, x_old.data()); + if (snap_points.size() > 0) { - vec3_t x_old; vec3_t x; - x_new[i_nodes] = vec3_t(0,0,0); - m_Grid->GetPoint(id_node, x_old.data()); + x_new[*iter] = vec3_t(0,0,0); double w_tot = 0; double L_min = 1e99; foreach (vtkIdType id_snap_node, snap_points) { m_Grid->GetPoint(id_snap_node, x.data()); double w = 1.0; w_tot += w; - x_new[i_nodes] += w*x; + x_new[*iter] += w*x; n += m_NodeNormal[id_snap_node]; double L = (x - x_old).abs(); L_min = min(L, L_min); } n.normalise(); - x_new[i_nodes] *= 1.0/w_tot; + x_new[*iter] *= 1.0/w_tot; + } else { + x_new[*iter] = x_old; + } - if (m_UseNormalCorrection) { - vec3_t dx = x_new[i_nodes] - x_old; - double scal = dx*m_NodeNormal[id_node]; - x_new[i_nodes] += scal*m_NodeNormal[id_node]; - } - vec3_t Dx = x_new[i_nodes] - x_old; - Dx *= m_UnderRelaxation; - if (moveNode(id_node, Dx)) { - x_new[i_nodes] = x_old + Dx; - } else { - x_new[i_nodes] = x_old; - m_Success = false; - } + if (m_UseNormalCorrection) { + vec3_t dx = x_new[*iter] - x_old; + double scal = dx*m_NodeNormal[id_node]; + x_new[*iter] += scal*m_NodeNormal[id_node]; + } + vec3_t Dx = x_new[*iter] - x_old; + Dx *= m_UnderRelaxation; + if (moveNode(id_node, Dx)) { + x_new[*iter] = x_old + Dx; + } else { + x_new[*iter] = x_old; + m_Success = false; } } } diff --git a/src/libengrid/laplacesmoother.h b/src/libengrid/laplacesmoother.h index 49dcd1b..f399111 100644 --- a/src/libengrid/laplacesmoother.h +++ b/src/libengrid/laplacesmoother.h @@ -25,6 +25,7 @@ #include "surfaceoperation.h" #include "surfaceprojection.h" +#include "surfacenodemovementcheck.h" class LaplaceSmoother : public SurfaceOperation { @@ -49,6 +50,8 @@ private: QSet m_AllowedCellTypes; QVector m_Fixed; + SurfaceNodeMovementCheck m_Check; + private: // methods diff --git a/src/libengrid/surfaceoperation.cpp b/src/libengrid/surfaceoperation.cpp index 70be49c..476613c 100644 --- a/src/libengrid/surfaceoperation.cpp +++ b/src/libengrid/surfaceoperation.cpp @@ -141,13 +141,66 @@ int SurfaceOperation::UpdateCurrentMeshDensity() return( 0 ); ///\todo what for??? } -int SurfaceOperation::UpdatePotentialSnapPoints(bool update_node_types, bool fix_unselected) +void SurfaceOperation::readVMD() { - setAllSurfaceCells(); + QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/table").replace("\n", " "); + int row_count = 0; + int column_count = 0; + m_VMDvector.clear(); + + if(!buffer.isEmpty()) { + QTextStream in(&buffer, QIODevice::ReadOnly); + in >> row_count >> column_count; + QVector tmp_bcs; + GuiMainWindow::pointer()->getAllBoundaryCodes(tmp_bcs); + if (column_count == tmp_bcs.size() + 3) { + m_VMDvector.fill(VertexMeshDensity(), row_count); + for (int i = 0; i < row_count; ++i) { + int row, column; + QString formula; + foreach (int bc, tmp_bcs) { + in >> row >> column >> formula; + m_VMDvector[row].BCmap[bc] = formula.toInt(); + } + in >> row >> column >> formula; + m_VMDvector[row].type = Str2VertexType(formula); + in >> row >> column >> formula; + if (formula == "{{{empty}}}") { + formula = ""; + } + m_VMDvector[i].setNodes(formula); + in >> row >> column >> formula; + m_VMDvector[i].density = formula.toDouble(); + } + } else { + EG_ERR_RETURN(QObject::tr("Mismatch of number of boundary codes!")); + } + } +} - l2g_t nodes = getPartNodes(); - l2g_t cells = getPartCells(); +void SurfaceOperation::updateNodeInfo() +{ + setAllCells(); + readVMD(); + l2g_t nodes = getPartNodes(); + computeNormals(); + foreach (vtkIdType id_node, nodes) { + EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type"); + node_type->SetValue(id_node, getNodeType(id_node, true)); + //density index from table + EG_VTKDCN(vtkIntArray, node_specified_density, m_Grid, "node_specified_density"); + + VertexMeshDensity nodeVMD = getVMD(id_node); + int idx = nodeVMD.findSmallestVMD(m_VMDvector); + node_specified_density->SetValue(id_node, idx); + } + updatePotentialSnapPoints(); +} + +void SurfaceOperation::updatePotentialSnapPoints() +{ + setAllSurfaceCells(); m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints()); if (m_UniformSnapPoints) { @@ -157,201 +210,135 @@ int SurfaceOperation::UpdatePotentialSnapPoints(bool update_node_types, bool fix m_PotentialSnapPoints[id_node].append(m_Part.n2nGG(id_node, i)); } } - return 0; + return; } - //initialize default values - EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" ); - foreach( vtkIdType id_node, nodes ) { - if ( update_node_types ) node_type->SetValue( id_node, EG_SIMPLE_VERTEX ); - m_PotentialSnapPoints[id_node].clear(); - } - - int num_edges = 0; - - // loop through edges - foreach (vtkIdType id_cell, cells) { - vtkIdType *pts, Npts; - m_Grid->GetCellPoints( id_cell, Npts, pts ); - for ( int i = 0; i < Npts; i++ ) { - - QVector EdgeCells; - int Ncells = getEdgeCells( pts[i], pts[(i+1)%Npts], EdgeCells ); - - bool already_visited_edge = false; - for(int i_cell=0;i_cellGetValue(id_node1) == EG_SIMPLE_VERTEX) - { - m_PotentialSnapPoints[id_node1].clear(); - m_PotentialSnapPoints[id_node1].push_back( id_node2 ); - if ( update_node_types ) node_type->SetValue( id_node1, edge ); - - } else if (( edge && node_type->GetValue(id_node1) == EG_BOUNDARY_EDGE_VERTEX) || - ( edge && node_type->GetValue(id_node1) == EG_FEATURE_EDGE_VERTEX) || - (!edge && node_type->GetValue(id_node1) == EG_SIMPLE_VERTEX)) + EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type"); + EG_FORALL_NODES(id_node1, m_Grid) { + m_PotentialSnapPoints[id_node1].clear(); + char type1 = node_type->GetValue(id_node1); + for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) { + vtkIdType id_node2 = m_Part.n2nGG(id_node1, i); + char type2 = node_type->GetValue(id_node2); + if ( (type1 == EG_SIMPLE_VERTEX) + || (type1 == EG_FEATURE_EDGE_VERTEX && (type2 == EG_FEATURE_EDGE_VERTEX || type2 == EG_FEATURE_CORNER_VERTEX)) + || (type1 == EG_BOUNDARY_EDGE_VERTEX && (type2 == EG_BOUNDARY_EDGE_VERTEX || type2 == EG_FIXED_VERTEX))) { - m_PotentialSnapPoints[id_node1].push_back( id_node2 ); - if ( node_type->GetValue( id_node1 ) && edge == EG_BOUNDARY_EDGE_VERTEX ) { - if ( update_node_types ) node_type->SetValue( id_node1, EG_BOUNDARY_EDGE_VERTEX );//EG_BOUNDARY_EDGE_VERTEX has priority over EG_FEATURE_EDGE_VERTEX - } - } - - if ( edge && node_type->GetValue( id_node2 ) == EG_SIMPLE_VERTEX ) { - m_PotentialSnapPoints[id_node2].clear(); - m_PotentialSnapPoints[id_node2].push_back( id_node1 ); - if ( update_node_types ) node_type->SetValue( id_node2, edge ); - } - else if (( edge && node_type->GetValue( id_node2 ) == EG_BOUNDARY_EDGE_VERTEX ) || - ( edge && node_type->GetValue( id_node2 ) == EG_FEATURE_EDGE_VERTEX ) || - ( !edge && node_type->GetValue( id_node2 ) == EG_SIMPLE_VERTEX ) ) { - m_PotentialSnapPoints[id_node2].push_back( id_node1 ); - if ( node_type->GetValue( id_node2 ) && edge == EG_BOUNDARY_EDGE_VERTEX ) { - if ( update_node_types ) node_type->SetValue( id_node2, EG_BOUNDARY_EDGE_VERTEX );//EG_BOUNDARY_EDGE_VERTEX has priority over EG_FEATURE_EDGE_VERTEX - } + m_PotentialSnapPoints[id_node1].append(id_node2); } } } - //cout<<"num_edges="<m_EdgeAngle); - //cout<<"===post-processing==="<GetValue( id_node ) == EG_FEATURE_EDGE_VERTEX || node_type->GetValue( id_node ) == EG_BOUNDARY_EDGE_VERTEX ) { //see how many edges; if two, what the angle is - - if ( !this->m_BoundarySmoothing && node_type->GetValue( id_node ) == EG_BOUNDARY_EDGE_VERTEX ) { - if ( update_node_types ) node_type->SetValue( id_node, EG_FIXED_VERTEX ); - } else if ( m_PotentialSnapPoints[id_node].size() != 2 ) { - if ( update_node_types ) node_type->SetValue( id_node, EG_FIXED_VERTEX ); - } - else { //check angle between edges - double x1[3], x2[3], x3[3], l1[3], l2[3]; - m_Grid->GetPoint( m_PotentialSnapPoints[id_node][0], x1 ); - m_Grid->GetPoint( id_node, x2 ); - m_Grid->GetPoint( m_PotentialSnapPoints[id_node][1], x3 ); - for ( int k = 0; k < 3; k++ ) { - l1[k] = x2[k] - x1[k]; - l2[k] = x3[k] - x2[k]; - } - if ( vtkMath::Normalize( l1 ) >= 0.0 && - vtkMath::Normalize( l2 ) >= 0.0 && - vtkMath::Dot( l1, l2 ) < CosEdgeAngle ) { - if ( update_node_types ) node_type->SetValue( id_node, EG_FIXED_VERTEX ); - } - }//if along edge - }//if edge vertex - } - //cout<<"m_PotentialSnapPoints.size()="<GetCellType(m_Part.n2cGG(id_node, i)) != VTK_TRIANGLE) { - return EG_FIXED_VERTEX; - } + if (m_Part.n2nGSize(id_node) == 0) { + return EG_FIXED_VERTEX; } - // experimental ... - EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type"); - char type = node_type->GetValue(id_node); - if (type != EG_SIMPLE_VERTEX) { - return type; - } - - // initialize default value //char type = EG_SIMPLE_VERTEX; + EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code"); - // loop through edges around id_node - - QVector edges; - double cos_edge_angle = cos(this->m_EdgeAngle); - - // safety switch to make sure edge angle fixing is disabled - if (cos_edge_angle < -0.99) { - cos_edge_angle = -2; - } - - foreach (int i_node2, n2n[_nodes[id_node]]) { + QSet bcs; + for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) { + vtkIdType id_cell = m_Part.n2cGG(id_node, i); - vtkIdType id_node2 = nodes[i_node2]; + // fix all vertices that are part of a volume element or a quad + if (m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) { + return EG_FIXED_VERTEX; + } - // determine edge type - char edge = getEdgeType(id_node2, id_node, fix_unselected); + int bc = cell_code->GetValue(id_cell); - // determine node type pre-processing - // count number of complex edges if the node is complex, - // otherwise, just count the number of edges - // - if (edge != EG_SIMPLE_VERTEX && type == EG_SIMPLE_VERTEX) { + // fix nodes which belong to faces with unselected boundary codes + if (!m_BoundaryCodes.contains(bc) && fix_unselected) { + return EG_FIXED_VERTEX; + } - edges.clear(); - edges.push_back(id_node2); - type = edge; + bcs.insert(bc); + } - } else if ((edge != EG_SIMPLE_VERTEX && type == EG_BOUNDARY_EDGE_VERTEX) || - (edge != EG_SIMPLE_VERTEX && type == EG_FEATURE_EDGE_VERTEX) || - (edge == EG_SIMPLE_VERTEX && type == EG_SIMPLE_VERTEX)) - { - edges.push_back(id_node2); - if (type != EG_SIMPLE_VERTEX && edge == EG_BOUNDARY_EDGE_VERTEX) { - // EG_BOUNDARY_EDGE_VERTEX has priority over EG_FEATURE_EDGE_VERTEX - type = EG_BOUNDARY_EDGE_VERTEX; - } - } + if (bcs.size() >= 3 || bcs.size() == 0) { + return EG_FIXED_VERTEX; + } + if (bcs.size() == 2) { + return EG_BOUNDARY_EDGE_VERTEX; } - // determine node type post-processing - if (type == EG_FEATURE_EDGE_VERTEX || type == EG_BOUNDARY_EDGE_VERTEX) { + SurfaceProjection* proj = GuiMainWindow::pointer()->getSurfProj(*bcs.begin(), true); + if (proj) { - // check how how many edges; if two, what is the angle? + vec3_t x, x0; + m_Grid->GetPoint(id_node, x0.data()); - if (!this->m_BoundarySmoothing && type == EG_BOUNDARY_EDGE_VERTEX) { - type = EG_FIXED_VERTEX; - } else if (edges.size() != 2) { - //type = EG_FIXED_VERTEX; + // compute average edge length + double L = 0; + for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) { + vec3_t xn; + m_Grid->GetPoint(m_Part.n2nGG(id_node, i), xn.data()); + L += (x0 - xn).abs(); + } + L /= m_Part.n2nGSize(id_node); + L *= 0.1; + bool convex = isConvexNode(id_node); + + x0 = proj->project(x0, id_node, true, m_NodeNormal[id_node]); + double radius = proj->lastProjRadius(); + if (convex) { + x = x0 + L*m_NodeNormal[id_node]; } else { + x = x0 - L*m_NodeNormal[id_node]; + } + vec3_t n = proj->lastProjNormal(); + if (GeometryTools::angle(n, m_NodeNormal[id_node]) > 0.5*m_FeatureAngle && !proj->lastProjFailed()) { + double d = L/tan(0.5*m_FeatureAngle); + static const int num_steps = 36; + double D_alpha = 2*M_PI/num_steps; + vec3_t v; + + v = GeometryTools::orthogonalVector(m_NodeNormal[id_node]); + int num_miss = 0; + for (int i = 0; i < num_steps; ++i) { + v = GeometryTools::rotate(v, m_NodeNormal[id_node], D_alpha); + vec3_t xp = proj->project(x, id_node, true, v, true); + if (proj->lastProjFailed()) { + ++num_miss; + } else { + double l = (x - xp).abs(); + if (l >= d) { + ++num_miss; + } + } + } + if (num_miss == 0) { + return EG_FEATURE_CORNER_VERTEX; + } - // check angle between edges - vec3_t x1, x2, x3, l1, l2; - m_Grid->GetPoint(edges[0], x1.data()); - m_Grid->GetPoint(id_node, x2.data()); - m_Grid->GetPoint(edges[1], x3.data()); - l1 = x2 - x1; - l2 = x3 - x2; - if (l1.abs() >= 0.0 && l2.abs() >= 0.0 && l1*l2 < cos_edge_angle) { - //type = EG_FIXED_VERTEX; + if (convex) { + x = x0 + L*n; + } else { + x = x0 - L*n; + } + v = GeometryTools::orthogonalVector(n); + int num_hit = 0; + for (int i = 0; i < num_steps; ++i) { + v = GeometryTools::rotate(v, n, D_alpha); + vec3_t xp = proj->project(x, id_node, true, v, true); + if (!proj->lastProjFailed()) { + double l = (x - xp).abs(); + if (l < d) { + ++num_hit; + } + } + } + if (num_hit > 0) { + return EG_FEATURE_EDGE_VERTEX; } - } // if along edge - } // if edge vertex + } + } - return type; + return EG_SIMPLE_VERTEX; } int SurfaceOperation::getEdgeCells(vtkIdType id_node1, vtkIdType id_node2, QVector &EdgeCells) @@ -414,7 +401,8 @@ char SurfaceOperation::getEdgeType(vtkIdType a_node1, vtkIdType a_node2, bool fi } else if (numNei == 1) { // check angle between cell1 and cell2 against FeatureAngle - if (cosAngle(m_Grid, neighbour_cells[0], neighbour_cells[1] ) <= cos_feature_angle && !feature_edges_disabled) { + double cosa = cosAngle(m_Grid, neighbour_cells[0], neighbour_cells[1]); + if (cosa <= cos_feature_angle && !feature_edges_disabled) { edge = EG_FEATURE_EDGE_VERTEX; } @@ -594,7 +582,7 @@ void SurfaceOperation::computeNormals() vec3_t n = u.cross(v); n.normalise(); if (checkVector(n)) { - normal[bcmap[bc]] += alpha*n; + normal[bcmap[bc]] -= alpha*n; } } } @@ -641,3 +629,25 @@ double SurfaceOperation::normalIrregularity(vtkIdType id_node) } return nirr; } + +bool SurfaceOperation::isConvexNode(vtkIdType id_node) +{ + int N = m_Part.n2nGSize(id_node); + if (N == 0) { + return false; + } + vec3_t x1, x2(0,0,0); + m_Grid->GetPoint(id_node, x1.data()); + for (int i = 0; i < N; ++i) { + vec3_t x; + m_Grid->GetPoint(m_Part.n2nGG(id_node, i), x.data()); + x2 += x; + } + x2 *= 1.0/N; + if ((x1 - x2)*m_NodeNormal[id_node] < 0) { + return true; + } + return false; +} + + diff --git a/src/libengrid/surfaceoperation.h b/src/libengrid/surfaceoperation.h index 01af493..a6a8274 100644 --- a/src/libengrid/surfaceoperation.h +++ b/src/libengrid/surfaceoperation.h @@ -47,9 +47,12 @@ class SurfaceOperation : public Operation private: - ///Vector used to store the "Potential Snap Points" of each point, i.e. the neighbour points belonging to the same edge (boundary or feature) in case of edge points, all neighbour points in case of simple points and the points belonging to edges in case of fixed points + /** Vector used to store the "Potential Snap Points" of each point, + * i.e. the neighbour points belonging to the same edge (boundary or feature) in case of edge points, + * all neighbour points in case of simple points and the points belonging to edges in case of fixed points */ QVector < QVector > m_PotentialSnapPoints; + void updatePotentialSnapPoints(); protected: // attributes @@ -59,13 +62,20 @@ protected: // attributes bool m_UniformSnapPoints; QVector m_NodeNormal; ///< node normal vectors + QVector m_VMDvector; + double m_StretchingFactor; protected: // methods - void computeNormals(); + void computeNormals(); + bool isConvexNode(vtkIdType id_node); + char geometricNodeType(vtkIdType id_node); double normalIrregularity(vtkIdType id_node); + void readVMD(); + void updateNodeInfo(); + public: @@ -76,14 +86,11 @@ public: int UpdateCurrentMeshDensity(); - /// Updates the m_PotentialSnapPoints structure + updates node types if desired (faster than loop through nodes with getNodeType) - int UpdatePotentialSnapPoints(bool update_node_types, bool fix_unselected = true); void setFeatureAngle(double FA) { m_FeatureAngle = FA; } void setEdgeAngle(double EA) { m_EdgeAngle = EA; } void setBoundarySmoothing(int BS) { m_BoundarySmoothing = BS; } - double currentVertexAvgDist(vtkIdType id_node); ///< Returns the average distance of id_node to its neighbours double CurrentMeshDensity( vtkIdType id_node ); ///< Returns 1/CurrentVertexAvgDist(id_node) char getNodeType(vtkIdType a_node, bool fix_unselected = true); ///< Returns the node type diff --git a/src/libengrid/swaptriangles.cpp b/src/libengrid/swaptriangles.cpp index 845c290..5964237 100644 --- a/src/libengrid/swaptriangles.cpp +++ b/src/libengrid/swaptriangles.cpp @@ -326,6 +326,7 @@ int SwapTriangles::swap() setAllSurfaceCells(); //computeAverageSurfaceError(); EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code"); + EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type"); QVector marked(m_Grid->GetNumberOfCells(), false); for (int i = 0; i < m_Part.getNumberOfCells(); ++i) { vtkIdType id_cell = m_Part.globalCell(i); @@ -333,7 +334,7 @@ int SwapTriangles::swap() if (!marked[id_cell] && !m_Swapped[id_cell]) { for (int j = 0; j < 3; ++j) { bool swap = false; - bool surface_error_swap = false; + bool feature_improvement_swap = false; stencil_t S = getStencil(id_cell, j); if(S.id_cell.size() == 2 && S.sameBC) { if (S.type_cell[1] == VTK_TRIANGLE) { @@ -360,23 +361,14 @@ int SwapTriangles::swap() force_swap = A1 < m_SmallAreaRatio*A2 || A2 < m_SmallAreaRatio*A1; } } - if (m_FeatureSwap || GeometryTools::angle(n1, n2) < m_FeatureAngle || force_swap) { - - // surface errors - double se1, se2; - bool surf_block = false; - computeSurfaceErrors(x3, cell_code->GetValue(id_cell), se1, se2); - if (se2 > se1) { // && se2 > m_SurfErrorHigherThreshold && se1 < m_SurfErrorLowerThreshold) { - surf_block = true; - } else { - if (se2 < se1) { // && se1 > m_SurfErrorHigherThreshold && se2 < m_SurfErrorLowerThreshold) { - swap = true; - surface_error_swap = true; - } + if (node_type->GetValue(S.p1) != EG_FEATURE_EDGE_VERTEX && node_type->GetValue(S.p2) != EG_FEATURE_EDGE_VERTEX) { + + if (node_type->GetValue(S.id_node[0]) == EG_FEATURE_EDGE_VERTEX && node_type->GetValue(S.id_node[1]) == EG_FEATURE_EDGE_VERTEX) { + //swap = true; } // Delaunay - if (!swap && !surf_block) { + if (!swap) { vec3_t n = n1 + n2; n.normalise(); vec3_t ex = orthogonalVector(n); @@ -411,7 +403,6 @@ int SwapTriangles::swap() } } } - //}// end of testswap } //end of if feature angle } //end of if l_marked } //end of if TestSwap @@ -419,7 +410,7 @@ int SwapTriangles::swap() } //end of S valid if (swap) { - if (testOrientation(S) || surface_error_swap) { + if (testOrientation(S)) { marked[S.id_cell[0]] = true; marked[S.id_cell[1]] = true; for (int k = 0; k < m_Part.n2cGSize(S.id_node[0]); ++k) { diff --git a/src/libengrid/utilities.cpp b/src/libengrid/utilities.cpp index 650a08c..9e65856 100644 --- a/src/libengrid/utilities.cpp +++ b/src/libengrid/utilities.cpp @@ -307,18 +307,20 @@ pair OrderedPair(vtkIdType a, vtkIdType b) { } const char* VertexType2Str(char T) { - if (T == EG_SIMPLE_VERTEX) return("EG_SIMPLE_VERTEX"); - if (T == EG_FIXED_VERTEX) return("EG_FIXED_VERTEX"); - if (T == EG_FEATURE_EDGE_VERTEX) return("EG_FEATURE_EDGE_VERTEX"); - if (T == EG_BOUNDARY_EDGE_VERTEX) return("EG_BOUNDARY_EDGE_VERTEX"); + if (T == EG_SIMPLE_VERTEX) return("EG_SIMPLE_VERTEX"); + if (T == EG_FIXED_VERTEX) return("EG_FIXED_VERTEX"); + if (T == EG_FEATURE_EDGE_VERTEX) return("EG_FEATURE_EDGE_VERTEX"); + if (T == EG_FEATURE_CORNER_VERTEX) return("EG_FEATURE_CORNER_VERTEX"); + if (T == EG_BOUNDARY_EDGE_VERTEX) return("EG_BOUNDARY_EDGE_VERTEX"); else return("Unknown vertex type"); } char Str2VertexType(QString S) { - if (S == "EG_SIMPLE_VERTEX") return(EG_SIMPLE_VERTEX); - if (S == "EG_FIXED_VERTEX") return(EG_FIXED_VERTEX); - if (S == "EG_FEATURE_EDGE_VERTEX") return(EG_FEATURE_EDGE_VERTEX); - if (S == "EG_BOUNDARY_EDGE_VERTEX") return(EG_BOUNDARY_EDGE_VERTEX); + if (S == "EG_SIMPLE_VERTEX") return(EG_SIMPLE_VERTEX); + if (S == "EG_FIXED_VERTEX") return(EG_FIXED_VERTEX); + if (S == "EG_FEATURE_EDGE_VERTEX") return(EG_FEATURE_EDGE_VERTEX); + if (S == "EG_FEATURE_CORNER_VERTEX") return(EG_FEATURE_CORNER_VERTEX); + if (S == "EG_BOUNDARY_EDGE_VERTEX") return(EG_BOUNDARY_EDGE_VERTEX); else return((char) - 1); } -- 2.11.4.GIT