diff --git a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H new file mode 100644 index 0000000000..22c87008f9 --- /dev/null +++ b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H @@ -0,0 +1,203 @@ +#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 * >¤tRing, + 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 * >¤tRing, + 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..fdffb7cac2 --- /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) 2011-2011 OpenCFD Ltd. + \\/ 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 typename HalfedgeDS::Traits Traits; + typedef typename 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..4a5082de7a --- /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) 2011-2011 OpenCFD Ltd. + \\/ 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/Make/files b/applications/utilities/surface/surfaceFeatureExtract/Make/files index a57125cd6a..0c0f6f7966 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/Make/files +++ b/applications/utilities/surface/surfaceFeatureExtract/Make/files @@ -1,3 +1,4 @@ surfaceFeatureExtract.C +CGALPolyhedron/buildCGALPolyhedron.C EXE = $(FOAM_APPBIN)/surfaceFeatureExtract diff --git a/applications/utilities/surface/surfaceFeatureExtract/Make/options b/applications/utilities/surface/surfaceFeatureExtract/Make/options index f87dafaa5d..ec38ffbbda 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/Make/options +++ b/applications/utilities/surface/surfaceFeatureExtract/Make/options @@ -1,9 +1,25 @@ +EXE_FROUNDING_MATH = -frounding-math +USE_F2C = -DCGAL_USE_F2C + +include $(GENERAL_RULES)/CGAL + EXE_INC = \ + ${EXE_FROUNDING_MATH} \ + ${USE_F2C} \ + ${CGAL_INC} \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/edgeMesh/lnInclude \ - -I$(LIB_SRC)/triSurface/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 + -ltriSurface \ + -lsampling diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C index 91bbd616cb..763fe38dd3 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) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -29,7 +29,6 @@ Description \*---------------------------------------------------------------------------*/ - #include "triangle.H" #include "triSurface.H" #include "argList.H" @@ -39,6 +38,15 @@ Description #include "treeBoundBox.H" #include "meshTools.H" #include "OFstream.H" +#include "triSurfaceMesh.H" +#include "vtkSurfaceWriter.H" +#include "triSurfaceFields.H" + +#include "buildCGALPolyhedron.H" +#include "CGALPolyhedronRings.H" +#include +#include +#include using namespace Foam; @@ -89,6 +97,154 @@ void deleteBox } +scalarField curvature(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::endl; + + // std::cout<< "k1 " << monge_form.principal_curvatures(0) << std::endl; + // std::cout<< "k2 " << monge_form.principal_curvatures(1) << std::endl; + + // std::vector::iterator itbp = in_points.begin(); + // std::vector::iterator itep = in_points.end(); + + // std::cout << "in_points list : " << std::endl; + + // for (; itbp != itep; itbp++) + // { + // std::cout << *itbp << std::endl; + // } + + // std::cout << "--- vertex " << vertI + // << " : " << v->point() << std::endl + // << "number of points used : " << in_points.size() + // << std::endl; + + k[vertI++] = Foam::sqrt + ( + sqr(monge_form.principal_curvatures(0)) + + sqr(monge_form.principal_curvatures(1)) + ); + } + + return k; +} + + // Main program: int main(int argc, char *argv[]) @@ -137,6 +293,22 @@ int main(int argc, char *argv[]) "((x0 y0 z0)(x1 y1 z1))", "remove edges within specified bounding box" ); + argList::addBoolOption + ( + "writeObj", + "write featureEdgeMesh obj files" + ); + argList::addOption + ( + "closeness", + "scalar", + "span to look for surface closeness" + ); + argList::addBoolOption + ( + "writeVTK", + "write surface property VTK files" + ); # include "setRootCase.H" # include "createTime.H" @@ -163,7 +335,6 @@ int main(int argc, char *argv[]) Info<< endl; - // Either construct features from surface&featureangle or read set. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -319,9 +490,338 @@ int main(int argc, char *argv[]) feMesh.write(); - feMesh.writeObj(sFeatFileName); + if (args.optionFound("writeObj")) + { + feMesh.writeObj(runTime.constant()/"featureEdgeMesh"/sFeatFileName); + }; - Info<< "End\n" << endl; + // Examine curvature, feature proximity and internal and external closeness. + + // Internal and external closeness + + triSurfaceMesh searchSurf + ( + IOobject + ( + sFeatFileName + ".closeness", + runTime.constant(), + "featureEdgeMesh", + runTime, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + surf + ); + + // Prepare start and end points for intersection tests + + const vectorField& normals = searchSurf.faceNormals(); + + faceList faces(surf.size()); + + forAll(surf, fI) + { + faces[fI] = surf[fI].triFaceFace(); + } + + scalar span = searchSurf.bounds().mag(); + + args.optionReadIfPresent("closeness", span); + + Info<< "span " << span << endl; + + pointField start = searchSurf.faceCentres() - span*normals; + pointField end = searchSurf.faceCentres() + span*normals; + const pointField& faceCentres = searchSurf.faceCentres(); + + // { + // label testI = 192; + + // List > testAllHitInfo; + + // searchSurf.findLineAll + // ( + // pointField(1, start[testI]), + // pointField(1, end[testI]), + // testAllHitInfo + // ); + + // Info<< testAllHitInfo << endl; + + // // return 0; + // } + + 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) + { + const List& hitInfo = allHitInfo[fI]; + + if (hitInfo.size() < 1) + { + Info<< nl << "# fI " << fI + << nl << "# start " << start[fI] + << nl << "# f centre " << faceCentres[fI] + << nl << "# end " << end[fI] + << endl; + + meshTools::writeOBJ(Info, start[fI]); + meshTools::writeOBJ(Info, faceCentres[fI]); + meshTools::writeOBJ(Info, end[fI]); + + Info<< "l 1 2 3" << endl; + + meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]); + meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]); + meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]); + + Info<< "f 4 5 6" << endl; + + FatalErrorIn(args.executable()) + << "findLineAll did not hit its own face." + << exit(FatalError); + } + else if (hitInfo.size() == 1) + { + if (!hitInfo[0].hit()) + { + FatalErrorIn(args.executable()) + << "findLineAll did not hit any face." + << exit(FatalError); + } + else if (hitInfo[0].index() != fI) + { + Info<< nl << "# findLineAll did not hit its own face." + << nl << "# fI " << fI + << nl << "# start " << start[fI] + << nl << "# f centre " << faceCentres[fI] + << nl << "# end " << end[fI] + << nl << "# hitInfo " << hitInfo + << endl; + + meshTools::writeOBJ(Info, start[fI]); + meshTools::writeOBJ(Info, faceCentres[fI]); + meshTools::writeOBJ(Info, end[fI]); + + Info<< "l 1 2 3" << endl; + + meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]); + meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]); + meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]); + + Info<< "f 4 5 6" << endl; + + forAll(hitInfo, hI) + { + label hFI = hitInfo[hI].index(); + + meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]); + meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]); + meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]); + + Info<< "f " + << 3*hI + 7 << " " + << 3*hI + 8 << " " + << 3*hI + 9 + << endl; + } + // FatalErrorIn(args.executable()) + // << "findLineAll did not hit its own face." + // << exit(FatalError); + } + } + else + { + label ownHitI = -1; + + forAll(hitInfo, hI) + { + // Find the hit on the triangle that launched the ray + + if (hitInfo[hI].index() == fI) + { + ownHitI = hI; + + break; + } + } + + if (ownHitI < 0) + { + Info<< nl << "# findLineAll did not hit its own face." + << nl << "# fI " << fI + << nl << "# start " << start[fI] + << nl << "# f centre " << faceCentres[fI] + << nl << "# end " << end[fI] + << nl << "# hitInfo " << hitInfo + << endl; + + meshTools::writeOBJ(Info, start[fI]); + meshTools::writeOBJ(Info, faceCentres[fI]); + meshTools::writeOBJ(Info, end[fI]); + + Info<< "l 1 2 3" << endl; + + meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]); + meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]); + meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]); + + Info<< "f 4 5 6" << endl; + + forAll(hitInfo, hI) + { + label hFI = hitInfo[hI].index(); + + meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]); + meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]); + meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]); + + Info<< "f " + << 3*hI + 7 << " " + << 3*hI + 8 << " " + << 3*hI + 9 + << endl; + } + + // 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 + + 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 + + internalCloseness[fI] = mag + ( + faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint() + ); + } + else + { + externalCloseness[fI] = mag + ( + faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint() + ); + + internalCloseness[fI] = mag + ( + faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint() + ); + } + } + } + + triSurfaceScalarField internalClosenessField + ( + IOobject + ( + sFeatFileName + ".internalCloseness", + runTime.constant(), + "featureEdgeMesh", + runTime, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + surf, + dimLength, + internalCloseness + ); + + internalClosenessField.write(); + + triSurfaceScalarField externalClosenessField + ( + IOobject + ( + sFeatFileName + ".externalCloseness", + runTime.constant(), + "featureEdgeMesh", + runTime, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + surf, + dimLength, + externalCloseness + ); + + externalClosenessField.write(); + + scalarField k = curvature(surf); + + triSurfacePointScalarField kField + ( + IOobject + ( + sFeatFileName + ".curvature", + runTime.constant(), + "featureEdgeMesh", + runTime, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + surf, + dimLength, + k + ); + + kField.write(); + + if (args.optionFound("writeObj")) + { + 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 + ); + + vtkSurfaceWriter().write + ( + runTime.constant()/"triSurface", // outputDir + sFeatFileName, // surfaceName + surf.points(), + faces, + "curvature", // fieldName + k, + true, // isNodeValues + true // verbose + ); + } return 0; }