diff --git a/applications/utilities/surface/surfaceFeatureExtract/Allwmake b/applications/utilities/surface/surfaceFeatureExtract/Allwmake new file mode 100755 index 0000000000..e54f15ae28 --- /dev/null +++ b/applications/utilities/surface/surfaceFeatureExtract/Allwmake @@ -0,0 +1,33 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + +set -x + +if [ ! -e "Make/files" ] || [ ! -e "Make/options" ] +then + mkdir -p Make + + if [ -n "$CGAL_ARCH_PATH" ] + then + cp -r MakeWithCGAL/* Make + + echo + echo Compiling surfaceFeatureExtract with CGAL support for curvature + echo + + wmake + else + cp -r MakeWithoutCGAL/* Make + + echo + echo Compiling surfaceFeatureExtract without CGAL support for curvature + echo + + wmake + fi +else + echo surfaceFeatureExtract already has a Make folder +fi + + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H new file mode 100644 index 0000000000..f5186c3a25 --- /dev/null +++ b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H @@ -0,0 +1,228 @@ +#ifndef CGAL_PSURF_RINGS_H_ +#define CGAL_PSURF_RINGS_H_ + +// This file adapted from +// CGAL-3.7/examples/Jet_fitting_3//PolyhedralSurf_rings.h +// Licensed under CGAL-3.7/LICENSE.FREE_USE + +// Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007 +// Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie +// Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France), +// Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute +// Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University +// (Israel). All rights reserved. + +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: + +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. + + +#include + +using namespace std; + +template +< + class TPoly, + class VertexPropertyMap +> + +class T_PolyhedralSurf_rings +{ + +protected: + + // Polyhedron + typedef typename TPoly::Vertex Vertex; + typedef typename TPoly::Halfedge Halfedge; + typedef typename TPoly::Facet Facet; + typedef typename TPoly::Halfedge_around_vertex_circulator + Halfedge_around_vertex_circulator; + typedef typename TPoly::Vertex_iterator Vertex_iterator; + + // vertex indices are initialised to -1 + static void reset_ring_indices + ( + std::vector< Vertex* >& vces, + VertexPropertyMap& vpm + ); + + // i >= 1; from a start vertex on the current i-1 ring, push non-visited + // neighbors of start in the nextRing and set indices to i. Also add these + // vertices in all. + static void push_neighbours_of + ( + Vertex* start, + int ith, + std::vector< Vertex* >& nextRing, + std::vector< Vertex* >& all, + VertexPropertyMap& vpm + ); + + // i >= 1, from a currentRing i-1, collect all neighbors, set indices + // to i and store them in nextRing and all. + static void collect_ith_ring + ( + int ith, + std::vector< Vertex* >& currentRing, + std::vector< Vertex* >& nextRing, + std::vector< Vertex* >& all, + VertexPropertyMap& vpm + ); + + +public: + + // collect i>=1 rings : all neighbours up to the ith ring, + static void collect_i_rings + ( + Vertex* v, + int ring_i, + std::vector< Vertex* >& all, + VertexPropertyMap& vpm + ); + + //collect enough rings (at least 1), to get at least min_nb of neighbors + static void collect_enough_rings + ( + Vertex* v, + unsigned int min_nb, + std::vector< Vertex* >& all, + VertexPropertyMap& vpm + ); +}; + + +////IMPLEMENTATION//////////////////////////////////////////////////// + +template < class TPoly , class VertexPropertyMap> +void T_PolyhedralSurf_rings :: +push_neighbours_of(Vertex* start, int ith, + std::vector< Vertex* >& nextRing, + std::vector< Vertex* >& all, + VertexPropertyMap& vpm) +{ + Vertex *v; + Halfedge_around_vertex_circulator + hedgeb = start->vertex_begin(), hedgee = hedgeb; + + CGAL_For_all(hedgeb, hedgee) + { + v = &*(hedgeb->opposite()->vertex()); + if (get(vpm, v) != -1) continue; //if visited: next + + put(vpm, v, ith); + nextRing.push_back(v); + all.push_back(v); + } +} + + +template +void T_PolyhedralSurf_rings :: +collect_ith_ring(int ith, std::vector< Vertex* >& currentRing, + std::vector< Vertex* >& nextRing, + std::vector< Vertex* >& all, + VertexPropertyMap& vpm) +{ + typename std::vector< Vertex* >::iterator + itb = currentRing.begin(), ite = currentRing.end(); + + CGAL_For_all(itb, ite) + { + push_neighbours_of(*itb, ith, nextRing, all, vpm); + } +} + + +template +void T_PolyhedralSurf_rings :: +reset_ring_indices(std::vector< Vertex* >& vces, + VertexPropertyMap& vpm) +{ + typename std::vector< Vertex* >::iterator + itb = vces.begin(), ite = vces.end(); + + CGAL_For_all(itb, ite) + { + put(vpm, *itb, -1); + } +} + + +template +void T_PolyhedralSurf_rings :: +collect_i_rings(Vertex* v, + int ring_i, + std::vector< Vertex* >& all, + VertexPropertyMap& vpm) +{ + std::vector current_ring, next_ring; + std::vector *p_current_ring, *p_next_ring; + assert(ring_i >= 1); + + //initialize + p_current_ring = ¤t_ring; + p_next_ring = &next_ring; + put(vpm, v, 0); + current_ring.push_back(v); + all.push_back(v); + + for (int i=1; i<=ring_i; i++) + { + collect_ith_ring(i, *p_current_ring, *p_next_ring, all, vpm); + //next round must be launched from p_nextRing... + p_current_ring->clear(); + std::swap(p_current_ring, p_next_ring); + } + + //clean up + reset_ring_indices(all, vpm); +} + + +template +void T_PolyhedralSurf_rings :: +collect_enough_rings(Vertex* v, + unsigned int min_nb, + std::vector< Vertex* >& all, + VertexPropertyMap& vpm) +{ + std::vector current_ring, next_ring; + std::vector *p_current_ring, *p_next_ring; + + //initialize + p_current_ring = ¤t_ring; + p_next_ring = &next_ring; + put(vpm, v, 0); + current_ring.push_back(v); + all.push_back(v); + + int i = 1; + + while ( (all.size() < min_nb) && (p_current_ring->size() != 0) ) + { + collect_ith_ring(i, *p_current_ring, *p_next_ring, all, vpm); + //next round must be launched from p_nextRing... + p_current_ring->clear(); + std::swap(p_current_ring, p_next_ring); + i++; + } + + //clean up + reset_ring_indices(all, vpm); +} + + +#endif diff --git a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C new file mode 100644 index 0000000000..63e2b904bf --- /dev/null +++ b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C @@ -0,0 +1,89 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "buildCGALPolyhedron.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::buildCGALPolyhedron::buildCGALPolyhedron +( + const Foam::triSurface& surf +) +: + CGAL::Modifier_base(), + surf_(surf) +{} + + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::buildCGALPolyhedron::~buildCGALPolyhedron() +{} + + +// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * // + +void Foam::buildCGALPolyhedron::operator() +( + HalfedgeDS& hds +) +{ + typedef HalfedgeDS::Traits Traits; + typedef Traits::Point_3 Point; + + // Postcondition: `hds' is a valid polyhedral surface. + CGAL::Polyhedron_incremental_builder_3 B(hds, false); + + B.begin_surface + ( + surf_.points().size(), // n points + surf_.size(), // n facets + 2*surf_.edges().size() // n halfedges + ); + + forAll(surf_.points(), pI) + { + const Foam::point& p = surf_.points()[pI]; + + B.add_vertex(Point(p.x(), p.y(), p.z())); + } + + forAll(surf_, fI) + { + B.begin_facet(); + + B.add_vertex_to_facet(surf_[fI][0]); + B.add_vertex_to_facet(surf_[fI][1]); + B.add_vertex_to_facet(surf_[fI][2]); + + B.end_facet(); + } + + B.end_surface(); +} + + +// ************************************************************************* // diff --git a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H new file mode 100644 index 0000000000..3001a247d3 --- /dev/null +++ b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H @@ -0,0 +1,106 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::buildCGALPolyhedron + +Description + Convert a triSurface into a CGAL Polyhedron + +SourceFiles + buildCGALPolyhedron.C + +\*---------------------------------------------------------------------------*/ + +#ifndef buildCGALPolyhedron_H +#define buildCGALPolyhedron_H + +#include "triSurface.H" +#include +#include +#include + +typedef CGAL::Simple_cartesian Kernel; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef Polyhedron::HalfedgeDS HalfedgeDS; +typedef Polyhedron::Vertex Vertex; +typedef Polyhedron::Vertex_iterator Vertex_iterator; +typedef Kernel::Point_3 Point_3; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class buildCGALPolyhedron Declaration +\*---------------------------------------------------------------------------*/ + +class buildCGALPolyhedron +: + public CGAL::Modifier_base +{ + // Private data + + //- Reference to triSurface to convert + const Foam::triSurface& surf_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + buildCGALPolyhedron(const buildCGALPolyhedron&); + + //- Disallow default bitwise assignment + void operator=(const buildCGALPolyhedron&); + + +public: + + // Constructors + + //- Construct with reference to triSurface + buildCGALPolyhedron(const triSurface& surf); + + + //- Destructor + ~buildCGALPolyhedron(); + + + // Member Operators + + //- operator() of this `modifier' called by delegate function of + // Polyhedron + void operator()(HalfedgeDS& hds); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files b/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files new file mode 100644 index 0000000000..0c0f6f7966 --- /dev/null +++ b/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files @@ -0,0 +1,4 @@ +surfaceFeatureExtract.C +CGALPolyhedron/buildCGALPolyhedron.C + +EXE = $(FOAM_APPBIN)/surfaceFeatureExtract diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options b/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options new file mode 100644 index 0000000000..60833572b7 --- /dev/null +++ b/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options @@ -0,0 +1,29 @@ +EXE_FROUNDING_MATH = -frounding-math +EXE_NDEBUG = -DNDEBUG +USE_F2C = -DCGAL_USE_F2C +include $(GENERAL_RULES)/CGAL + +EXE_INC = \ + -DENABLE_CURVATURE \ + ${EXE_FROUNDING_MATH} \ + ${EXE_NDEBUG} \ + ${USE_F2C} \ + ${CGAL_INC} \ + -ICGALPolyhedron \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/edgeMesh/lnInclude \ + -I$(LIB_SRC)/triSurface/lnInclude \ + -I$(LIB_SRC)/surfMesh/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + $(CGAL_LIBS) \ + -L$(CGAL_ARCH_PATH)/lib \ + -llapack \ + -lblas \ + -lCGAL \ + -lmeshTools \ + -ledgeMesh \ + -ltriSurface \ + -lsampling diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files b/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files new file mode 100644 index 0000000000..a57125cd6a --- /dev/null +++ b/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files @@ -0,0 +1,3 @@ +surfaceFeatureExtract.C + +EXE = $(FOAM_APPBIN)/surfaceFeatureExtract diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options b/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options new file mode 100644 index 0000000000..b733d5fde9 --- /dev/null +++ b/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options @@ -0,0 +1,13 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/edgeMesh/lnInclude \ + -I$(LIB_SRC)/triSurface/lnInclude \ + -I$(LIB_SRC)/surfMesh/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + -lmeshTools \ + -ledgeMesh \ + -ltriSurface \ + -lsampling diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C index 5f6ee05fb1..c2baec519e 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C +++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -46,10 +46,263 @@ Description #include "unitConversion.H" #include "plane.H" +#ifdef ENABLE_CURVATURE +#include "buildCGALPolyhedron.H" +#include "CGALPolyhedronRings.H" +#include +#include +#include +#endif + using namespace Foam; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#ifdef ENABLE_CURVATURE +scalarField calcCurvature(const triSurface& surf) +{ + scalarField k(surf.points().size(), 0); + + Polyhedron P; + + buildCGALPolyhedron convert(surf); + P.delegate(convert); + + // Info<< "Created CGAL Polyhedron with " << label(P.size_of_vertices()) + // << " vertices and " << label(P.size_of_facets()) + // << " facets. " << endl; + + // The rest of this function adapted from + // CGAL-3.7/examples/Jet_fitting_3/Mesh_estimation.cpp + // Licensed under CGAL-3.7/LICENSE.FREE_USE + + // Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007 + // Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie + // Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France), + // Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute + // Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University + // (Israel). All rights reserved. + + // Permission is hereby granted, free of charge, to any person obtaining a + // copy of this software and associated documentation files (the + // "Software"), to deal in the Software without restriction, including + // without limitation the rights to use, copy, modify, merge, publish, + // distribute, sublicense, and/or sell copies of the Software, and to permit + // persons to whom the Software is furnished to do so, subject to the + // following conditions: + + // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + // IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + // CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT + // OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR + // THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + //Vertex property map, with std::map + typedef std::map Vertex2int_map_type; + typedef boost::associative_property_map< Vertex2int_map_type > + Vertex_PM_type; + typedef T_PolyhedralSurf_rings Poly_rings; + + typedef CGAL::Monge_via_jet_fitting Monge_via_jet_fitting; + typedef Monge_via_jet_fitting::Monge_form Monge_form; + + std::vector in_points; //container for data points + + // default parameter values and global variables + unsigned int d_fitting = 2; + unsigned int d_monge = 2; + unsigned int min_nb_points = (d_fitting + 1)*(d_fitting + 2)/2; + + //initialize the tag of all vertices to -1 + Vertex_iterator vitb = P.vertices_begin(); + Vertex_iterator vite = P.vertices_end(); + + Vertex2int_map_type vertex2props; + Vertex_PM_type vpm(vertex2props); + + CGAL_For_all(vitb, vite) + { + put(vpm, &(*vitb), -1); + } + + vite = P.vertices_end(); + + label vertI = 0; + + for (vitb = P.vertices_begin(); vitb != vite; vitb++) + { + //initialize + Vertex* v = &(*vitb); + + //gather points around the vertex using rings + // From: gather_fitting_points(v, in_points, vpm); + { + std::vector gathered; + in_points.clear(); + + Poly_rings::collect_enough_rings(v, min_nb_points, gathered, vpm); + + //store the gathered points + std::vector::iterator itb = gathered.begin(); + std::vector::iterator ite = gathered.end(); + + CGAL_For_all(itb, ite) + { + in_points.push_back((*itb)->point()); + } + } + + //skip if the nb of points is to small + if ( in_points.size() < min_nb_points ) + { + std::cerr + << "not enough pts for fitting this vertex" + << in_points.size() + << std::endl; + + continue; + } + + // perform the fitting + Monge_via_jet_fitting monge_fit; + + Monge_form monge_form = monge_fit + ( + in_points.begin(), + in_points.end(), + d_fitting, + d_monge + ); + +// std::cout<< monge_form;; +// std::cout<< "condition number : " +// << monge_fit.condition_number() << nl << std::endl; + + // Use the maximum curvature to give smaller cell sizes later. + k[vertI++] + = max + ( + mag(monge_form.principal_curvatures(0)), + mag(monge_form.principal_curvatures(1)) + ); + } + + return k; +} +#endif + + +bool edgesConnected(const edge& e1, const edge& e2) +{ + if + ( + e1.start() == e2.start() || e1.start() == e2.end() + || e1.end() == e2.start() || e1.end() == e2.end() + ) + { + return true; + } + + return false; +} + + +scalar calcProximityOfFeaturePoints +( + const List& hitList, + const scalar defaultCellSize +) +{ + scalar minDist = defaultCellSize; + + for + ( + label hI1 = 0; + hI1 < hitList.size() - 1; + ++hI1 + ) + { + const pointIndexHit& pHit1 = hitList[hI1]; + + if (pHit1.hit()) + { + for + ( + label hI2 = hI1 + 1; + hI2 < hitList.size(); + ++hI2 + ) + { + const pointIndexHit& pHit2 = hitList[hI2]; + + if (pHit2.hit()) + { + scalar curDist = mag(pHit1.hitPoint() - pHit2.hitPoint()); + + minDist = min(curDist, minDist); + } + } + } + } + + return minDist; +} + + +scalar calcProximityOfFeatureEdges +( + const extendedFeatureEdgeMesh& efem, + const List& hitList, + const scalar defaultCellSize +) +{ + scalar minDist = defaultCellSize; + + for + ( + label hI1 = 0; + hI1 < hitList.size() - 1; + ++hI1 + ) + { + const pointIndexHit& pHit1 = hitList[hI1]; + + if (pHit1.hit()) + { + const edge& e1 = efem.edges()[pHit1.index()]; + + for + ( + label hI2 = hI1 + 1; + hI2 < hitList.size(); + ++hI2 + ) + { + const pointIndexHit& pHit2 = hitList[hI2]; + + if (pHit2.hit()) + { + const edge& e2 = efem.edges()[pHit2.index()]; + + // Don't refine if the edges are connected to each other + if (!edgesConnected(e1, e2)) + { + scalar curDist + = mag(pHit1.hitPoint() - pHit2.hitPoint()); + + minDist = min(curDist, minDist); + } + } + } + } + } + + return minDist; +} + + void dumpBox(const treeBoundBox& bb, const fileName& fName) { OFstream str(fName); @@ -299,10 +552,9 @@ int main(int argc, char *argv[]) "writeVTK", "write extendedFeatureEdgeMesh vtk files" ); - argList::addOption + argList::addBoolOption ( "closeness", - "scalar", "span to look for surface closeness" ); argList::addOption @@ -331,7 +583,7 @@ int main(int argc, char *argv[]) # ifdef ENABLE_CURVATURE argList::addBoolOption ( - "calcCurvature", + "curvature", "calculate curvature and closeness fields" ); # endif @@ -527,6 +779,7 @@ int main(int argc, char *argv[]) << " will be included as feature edges."<< endl; } + surfaceFeatures newSet(surf); newSet.setFromStatus(edgeStat); @@ -564,7 +817,6 @@ int main(int argc, char *argv[]) feMesh.write(); - // Write a featureEdgeMesh for backwards compatibility { featureEdgeMesh bfeMesh @@ -595,7 +847,7 @@ int main(int argc, char *argv[]) ( sFeatFileName + ".closeness", runTime.constant(), - "extendedFeatureEdgeMesh", + "triSurface", runTime, IOobject::NO_READ, IOobject::NO_WRITE @@ -603,13 +855,6 @@ int main(int argc, char *argv[]) surf ); - if (!curvature) - { - Info<< "End\n" << endl; - - return 0; - } - // Find close features // // Dummy trim operation to mark features @@ -677,70 +922,56 @@ int main(int argc, char *argv[]) // ) // ); - Info<< "Examine curvature, feature proximity and internal and " - << "external closeness." << endl; - - // Internal and external closeness - - // Prepare start and end points for intersection tests - - const vectorField& normals = searchSurf.faceNormals(); - - scalar span = searchSurf.bounds().mag(); - - args.optionReadIfPresent("closeness", span); - - scalar externalAngleTolerance = 10; - scalar externalToleranceCosAngle = Foam::cos - ( - degToRad(180 - externalAngleTolerance) - ); - - scalar internalAngleTolerance = 45; - scalar internalToleranceCosAngle = Foam::cos - ( - degToRad(180 - internalAngleTolerance) - ); - - Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle << nl - << "internalToleranceCosAngle: " << internalToleranceCosAngle - << endl; - - // Info<< "span " << span << endl; - - pointField start = searchSurf.faceCentres() - span*normals; - pointField end = searchSurf.faceCentres() + span*normals; - const pointField& faceCentres = searchSurf.faceCentres(); - - List > allHitInfo; - - // Find all intersections (in order) - searchSurf.findLineAll(start, end, allHitInfo); - - scalarField internalCloseness(start.size(), GREAT); - scalarField externalCloseness(start.size(), GREAT); - - forAll(allHitInfo, fI) + if (args.optionFound("closeness")) { - const List& hitInfo = allHitInfo[fI]; + Info<< nl << "Extracting internal and external closeness of surface." + << endl; - if (hitInfo.size() < 1) - { - drawHitProblem(fI, surf, start, faceCentres, end, hitInfo); + // Internal and external closeness - // FatalErrorIn(args.executable()) - // << "findLineAll did not hit its own face." - // << exit(FatalError); - } - else if (hitInfo.size() == 1) + // Prepare start and end points for intersection tests + + const vectorField& normals = searchSurf.faceNormals(); + + scalar span = searchSurf.bounds().mag(); + + //args.optionReadIfPresent("closeness", span); + + scalar externalAngleTolerance = 10; + scalar externalToleranceCosAngle = Foam::cos + ( + degToRad(180 - externalAngleTolerance) + ); + + scalar internalAngleTolerance = 45; + scalar internalToleranceCosAngle = Foam::cos + ( + degToRad(180 - internalAngleTolerance) + ); + + Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle << nl + << "internalToleranceCosAngle: " << internalToleranceCosAngle + << endl; + + // Info<< "span " << span << endl; + + pointField start = searchSurf.faceCentres() - span*normals; + pointField end = searchSurf.faceCentres() + span*normals; + const pointField& faceCentres = searchSurf.faceCentres(); + + List > allHitInfo; + + // Find all intersections (in order) + searchSurf.findLineAll(start, end, allHitInfo); + + scalarField internalCloseness(start.size(), GREAT); + scalarField externalCloseness(start.size(), GREAT); + + forAll(allHitInfo, fI) { - if (!hitInfo[0].hit()) - { - // FatalErrorIn(args.executable()) - // << "findLineAll did not hit any face." - // << exit(FatalError); - } - else if (hitInfo[0].index() != fI) + const List& hitInfo = allHitInfo[fI]; + + if (hitInfo.size() < 1) { drawHitProblem(fI, surf, start, faceCentres, end, hitInfo); @@ -748,202 +979,299 @@ int main(int argc, char *argv[]) // << "findLineAll did not hit its own face." // << exit(FatalError); } - } - else - { - label ownHitI = -1; - - forAll(hitInfo, hI) + else if (hitInfo.size() == 1) { - // Find the hit on the triangle that launched the ray - - if (hitInfo[hI].index() == fI) + if (!hitInfo[0].hit()) { - ownHitI = hI; - - break; + // FatalErrorIn(args.executable()) + // << "findLineAll did not hit any face." + // << exit(FatalError); } - } - - if (ownHitI < 0) - { - drawHitProblem(fI, surf, start, faceCentres, end, hitInfo); - - // FatalErrorIn(args.executable()) - // << "findLineAll did not hit its own face." - // << exit(FatalError); - } - else if (ownHitI == 0) - { - // There are no internal hits, the first hit is the closest - // external hit - - if - ( - (normals[fI] & normals[hitInfo[ownHitI + 1].index()]) - < externalToleranceCosAngle - ) + else if (hitInfo[0].index() != fI) { - externalCloseness[fI] = mag - ( - faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint() - ); - } - } - else if (ownHitI == hitInfo.size() - 1) - { - // There are no external hits, the last but one hit is the - // closest internal hit + drawHitProblem(fI, surf, start, faceCentres, end, hitInfo); - if - ( - (normals[fI] & normals[hitInfo[ownHitI - 1].index()]) - < internalToleranceCosAngle - ) - { - internalCloseness[fI] = mag - ( - faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint() - ); + // FatalErrorIn(args.executable()) + // << "findLineAll did not hit its own face." + // << exit(FatalError); } } else { - if - ( - (normals[fI] & normals[hitInfo[ownHitI + 1].index()]) - < externalToleranceCosAngle - ) + label ownHitI = -1; + + forAll(hitInfo, hI) { - externalCloseness[fI] = mag - ( - faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint() - ); + // Find the hit on the triangle that launched the ray + + if (hitInfo[hI].index() == fI) + { + ownHitI = hI; + + break; + } } - if - ( - (normals[fI] & normals[hitInfo[ownHitI - 1].index()]) - < internalToleranceCosAngle - ) + if (ownHitI < 0) { - internalCloseness[fI] = mag + drawHitProblem(fI, surf, start, faceCentres, end, hitInfo); + + // FatalErrorIn(args.executable()) + // << "findLineAll did not hit its own face." + // << exit(FatalError); + } + else if (ownHitI == 0) + { + // There are no internal hits, the first hit is the closest + // external hit + + if ( - faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint() - ); + (normals[fI] & normals[hitInfo[ownHitI + 1].index()]) + < externalToleranceCosAngle + ) + { + externalCloseness[fI] = mag + ( + faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint() + ); + } + } + else if (ownHitI == hitInfo.size() - 1) + { + // There are no external hits, the last but one hit is the + // closest internal hit + + if + ( + (normals[fI] & normals[hitInfo[ownHitI - 1].index()]) + < internalToleranceCosAngle + ) + { + internalCloseness[fI] = mag + ( + faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint() + ); + } + } + else + { + if + ( + (normals[fI] & normals[hitInfo[ownHitI + 1].index()]) + < externalToleranceCosAngle + ) + { + externalCloseness[fI] = mag + ( + faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint() + ); + } + + if + ( + (normals[fI] & normals[hitInfo[ownHitI - 1].index()]) + < internalToleranceCosAngle + ) + { + internalCloseness[fI] = mag + ( + faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint() + ); + } } } } + + triSurfaceScalarField internalClosenessField + ( + IOobject + ( + sFeatFileName + ".internalCloseness", + runTime.constant(), + "triSurface", + runTime, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + surf, + dimLength, + internalCloseness + ); + + internalClosenessField.write(); + + triSurfaceScalarField externalClosenessField + ( + IOobject + ( + sFeatFileName + ".externalCloseness", + runTime.constant(), + "triSurface", + runTime, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + surf, + dimLength, + externalCloseness + ); + + externalClosenessField.write(); + + if (writeVTK) + { + vtkSurfaceWriter().write + ( + runTime.constant()/"triSurface", // outputDir + sFeatFileName, // surfaceName + surf.points(), + faces, + "internalCloseness", // fieldName + internalCloseness, + false, // isNodeValues + true // verbose + ); + + vtkSurfaceWriter().write + ( + runTime.constant()/"triSurface", // outputDir + sFeatFileName, // surfaceName + surf.points(), + faces, + "externalCloseness", // fieldName + externalCloseness, + false, // isNodeValues + true // verbose + ); + } } - triSurfaceScalarField internalClosenessField - ( - IOobject - ( - sFeatFileName + ".internalCloseness", - runTime.constant(), - "extendedFeatureEdgeMesh", - runTime, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - surf, - dimLength, - internalCloseness - ); - - internalClosenessField.write(); - - triSurfaceScalarField externalClosenessField - ( - IOobject - ( - sFeatFileName + ".externalCloseness", - runTime.constant(), - "extendedFeatureEdgeMesh", - runTime, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - surf, - dimLength, - externalCloseness - ); - - externalClosenessField.write(); - #ifdef ENABLE_CURVATURE - scalarField k = calcCurvature(surf); - - // Modify the curvature values on feature edges and points to be zero. - - forAll(newSet.featureEdges(), fEI) + if (args.optionFound("curvature")) { - const edge& e = surf.edges()[newSet.featureEdges()[fEI]]; + Info<< nl << "Extracting curvature of surface at the points." << endl; - k[surf.meshPoints()[e.start()]] = 0.0; - k[surf.meshPoints()[e.end()]] = 0.0; - } + scalarField k = calcCurvature(surf); - triSurfacePointScalarField kField - ( - IOobject + // Modify the curvature values on feature edges and points to be zero. + + // forAll(newSet.featureEdges(), fEI) + // { + // const edge& e = surf.edges()[newSet.featureEdges()[fEI]]; + // + // k[surf.meshPoints()[e.start()]] = 0.0; + // k[surf.meshPoints()[e.end()]] = 0.0; + // } + + triSurfacePointScalarField kField ( - sFeatFileName + ".curvature", - runTime.constant(), - "extendedFeatureEdgeMesh", - runTime, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - surf, - dimLength, - k - ); + IOobject + ( + sFeatFileName + ".curvature", + runTime.constant(), + "triSurface", + runTime, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + surf, + dimLength, + k + ); - kField.write(); + kField.write(); + + if (writeVTK) + { + vtkSurfaceWriter().write + ( + runTime.constant()/"triSurface", // outputDir + sFeatFileName, // surfaceName + surf.points(), + faces, + "curvature", // fieldName + k, + true, // isNodeValues + true // verbose + ); + } + } #endif - if (writeVTK) + + if (args.optionFound("featureProximity")) { - vtkSurfaceWriter().write + Info<< nl << "Extracting proximity of close feature points and edges " + << "to the surface" << endl; + + const scalar searchDistance = + args.optionRead("featureProximity"); + + const scalar radiusSqr = sqr(searchDistance); + + scalarField featureProximity(surf.size(), searchDistance); + + forAll(surf, fI) + { + const triPointRef& tri = surf[fI].tri(surf.points()); + const point& triCentre = tri.circumCentre(); + + List hitList; + + feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList); + + featureProximity[fI] = + calcProximityOfFeatureEdges + ( + feMesh, + hitList, + featureProximity[fI] + ); + + feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList); + + featureProximity[fI] = + calcProximityOfFeaturePoints + ( + hitList, + featureProximity[fI] + ); + } + + triSurfaceScalarField featureProximityField ( - runTime.constant()/"triSurface", // outputDir - sFeatFileName, // surfaceName - surf.points(), - faces, - "internalCloseness", // fieldName - internalCloseness, - false, // isNodeValues - true // verbose + IOobject + ( + sFeatFileName + ".featureProximity", + runTime.constant(), + "triSurface", + runTime, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + surf, + dimLength, + featureProximity ); - vtkSurfaceWriter().write - ( - runTime.constant()/"triSurface", // outputDir - sFeatFileName, // surfaceName - surf.points(), - faces, - "externalCloseness", // fieldName - externalCloseness, - false, // isNodeValues - true // verbose - ); + featureProximityField.write(); -# ifdef ENABLE_CURVATURE - vtkSurfaceWriter().write - ( - runTime.constant()/"triSurface", // outputDir - sFeatFileName, // surfaceName - surf.points(), - faces, - "curvature", // fieldName - k, - true, // isNodeValues - true // verbose - ); -# endif + if (writeVTK) + { + vtkSurfaceWriter().write + ( + runTime.constant()/"triSurface", // outputDir + sFeatFileName, // surfaceName + surf.points(), + faces, + "featureProximity", // fieldName + featureProximity, + false, // isNodeValues + true // verbose + ); + } } Info<< "End\n" << endl; diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C index 81f853aa41..cd0ac4fc94 100644 --- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C +++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C @@ -819,6 +819,39 @@ void Foam::extendedFeatureEdgeMesh::nearestFeatureEdgeByType } +void Foam::extendedFeatureEdgeMesh::allNearestFeaturePoints +( + const point& sample, + scalar searchRadiusSqr, + List& info +) const +{ + DynamicList dynPointHit; + + // Pick up all the feature points that intersect the search sphere + labelList elems = pointTree().findSphere + ( + sample, + searchRadiusSqr + ); + + forAll(elems, elemI) + { + label index = elems[elemI]; + label ptI = pointTree().shapes().pointLabels()[index]; + const point& pt = points()[ptI]; + + pointIndexHit nearHit; + + nearHit = pointIndexHit(true, pt, index); + + dynPointHit.append(nearHit); + } + + info.transfer(dynPointHit); +} + + void Foam::extendedFeatureEdgeMesh::allNearestFeatureEdges ( const point& sample, diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H index 52de62552c..28601a26dd 100644 --- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H +++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H @@ -283,6 +283,14 @@ public: List& info ) const; + //- Find all the feature points within searchDistSqr of sample + void allNearestFeaturePoints + ( + const point& sample, + scalar searchRadiusSqr, + List& info + ) const; + //- Find all the feature edges within searchDistSqr of sample void allNearestFeatureEdges ( @@ -291,6 +299,7 @@ public: List& info ) const; + // Access //- Return the index of the start of the convex feature points diff --git a/src/meshTools/indexedOctree/treeDataPoint.C b/src/meshTools/indexedOctree/treeDataPoint.C index 7e381a7577..09bf90de38 100644 --- a/src/meshTools/indexedOctree/treeDataPoint.C +++ b/src/meshTools/indexedOctree/treeDataPoint.C @@ -93,6 +93,25 @@ bool Foam::treeDataPoint::overlaps } +// Check if any point on shape is inside sphere. +bool Foam::treeDataPoint::overlaps +( + const label index, + const point& centre, + const scalar radiusSqr +) const +{ + label pointI = (useSubset_ ? pointLabels_[index] : index); + + if (magSqr(points_[pointI] - centre) <= radiusSqr) + { + return true; + } + + return false; +} + + // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex, // nearestPoint. void Foam::treeDataPoint::findNearest diff --git a/src/meshTools/indexedOctree/treeDataPoint.H b/src/meshTools/indexedOctree/treeDataPoint.H index abd0a4024a..f867a3274a 100644 --- a/src/meshTools/indexedOctree/treeDataPoint.H +++ b/src/meshTools/indexedOctree/treeDataPoint.H @@ -128,6 +128,14 @@ public: const treeBoundBox& sampleBb ) const; + //- Does shape at index overlap the sphere + bool overlaps + ( + const label index, + const point& centre, + const scalar radiusSqr + ) const; + //- Calculates nearest (to sample) point in shape. // Returns actual point and distance (squared) void findNearest