diff --git a/applications/utilities/surface/surfaceAutoPatch/Make/files b/applications/utilities/surface/surfaceAutoPatch/Make/files
deleted file mode 100644
index 38a0ebbbf2a913b022dd9b4302f29b4eeee0984f..0000000000000000000000000000000000000000
--- a/applications/utilities/surface/surfaceAutoPatch/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-surfaceAutoPatch.C
-
-EXE = $(FOAM_APPBIN)/surfaceAutoPatch
diff --git a/applications/utilities/surface/surfaceAutoPatch/Make/options b/applications/utilities/surface/surfaceAutoPatch/Make/options
deleted file mode 100644
index 9f08e8d2a80f6f1fb269c1815170453ae49cf783..0000000000000000000000000000000000000000
--- a/applications/utilities/surface/surfaceAutoPatch/Make/options
+++ /dev/null
@@ -1,7 +0,0 @@
-EXE_INC = \
- -I$(LIB_SRC)/meshTools/lnInclude \
- -I$(LIB_SRC)/triSurface/lnInclude
-
-EXE_LIBS = \
- -lmeshTools \
- -ltriSurface
diff --git a/applications/utilities/surface/surfaceAutoPatch/surfaceAutoPatch.C b/applications/utilities/surface/surfaceAutoPatch/surfaceAutoPatch.C
deleted file mode 100644
index f53ec829635b309221fda439443aa83e1bd820eb..0000000000000000000000000000000000000000
--- a/applications/utilities/surface/surfaceAutoPatch/surfaceAutoPatch.C
+++ /dev/null
@@ -1,125 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 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 <http://www.gnu.org/licenses/>.
-
-Application
- surfaceAutoPatch
-
-Description
- Patches surface according to feature angle. Like autoPatch.
-
-\*---------------------------------------------------------------------------*/
-
-#include "triangle.H"
-#include "triSurface.H"
-#include "argList.H"
-#include "surfaceFeatures.H"
-#include "treeBoundBox.H"
-#include "meshTools.H"
-#include "OFstream.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-int main(int argc, char *argv[])
-{
- argList::noParallel();
- argList::validArgs.append("input surfaceFile");
- argList::validArgs.append("output surfaceFile");
- argList::validArgs.append("includedAngle [0..180]");
- argList args(argc, argv);
-
- const fileName inFileName = args[1];
- const fileName outFileName = args[2];
- const scalar includedAngle = args.argRead<scalar>(3);
-
- Info<< "Surface : " << inFileName << nl << endl;
-
-
- // Read
- // ~~~~
-
- Info<< "Reading : " << inFileName << endl;
- triSurface surf(inFileName);
-
- Info<< "Read surface:" << endl;
- surf.writeStats(Info);
- Info<< endl;
-
-
-
- // Construct features from surface&featureangle
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- Info<< "Constructing feature set from included angle " << includedAngle
- << endl;
-
- surfaceFeatures set(surf, includedAngle);
-
- Info<< nl
- << "Feature set:" << nl
- << " feature points : " << set.featurePoints().size() << nl
- << " feature edges : " << set.featureEdges().size() << nl
- << " of which" << nl
- << " region edges : " << set.nRegionEdges() << nl
- << " external edges : " << set.nExternalEdges() << nl
- << " internal edges : " << set.nInternalEdges() << nl
- << endl;
-
- // Get per-edge status.
- boolList borderEdge(surf.nEdges(), false);
- forAll(set.featureEdges(), i)
- {
- borderEdge[set.featureEdges()[i]] = true;
- }
-
- labelList faceRegion(surf.size());
- label nRegions = surf.markZones(borderEdge, faceRegion);
-
- // Reregion triangles.
- forAll(surf, i)
- {
- surf[i].region() = faceRegion[i];
- }
-
- // Create some patches
- surf.patches().setSize(nRegions);
-
- forAll(surf.patches(), patchI)
- {
- surf.patches()[patchI].name() = "patch" + Foam::name(patchI);
- surf.patches()[patchI].geometricType() = "empty";
- }
-
-
- Info<< "Writing : " << outFileName << endl;
- surf.write(outFileName, true);
-
- Info<< "End\n" << endl;
-
- return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/Allwmake b/applications/utilities/surface/surfaceBooleanFeatures/Allwmake
new file mode 100755
index 0000000000000000000000000000000000000000..50d640f2b259e90e7b898ee8e01761780709b6e6
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/Allwmake
@@ -0,0 +1,15 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+set -x
+
+if [ -z "$CGAL_ARCH_PATH" ]
+then
+ export COMPILE_FLAGS="-DNO_CGAL"
+else
+ wmake PolyhedronReader
+ export COMPILE_FLAGS='-IPolyhedronReader'
+ export LINK_FLAGS='${CGAL_LIBS} -lPolyhedronReader'
+fi
+wmake
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/CGAL3DKernel.H b/applications/utilities/surface/surfaceBooleanFeatures/CGAL3DKernel.H
new file mode 100644
index 0000000000000000000000000000000000000000..e753060ce0f161ee68bdc4e2c5fd29eb53ada8ea
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/CGAL3DKernel.H
@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+Typedefs
+ CGAL3DKernel
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CGAL3DKernel_H
+#define CGAL3DKernel_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef CGAL_INEXACT
+
+ // Fast kernel using a double as the storage type
+ #include "CGAL/Exact_predicates_inexact_constructions_kernel.h"
+ typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+#else
+
+ // Very robust but expensive kernel
+ #include "CGAL/Exact_predicates_exact_constructions_kernel.h"
+ typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+
+#endif
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/CGALIndexedPolyhedron.H b/applications/utilities/surface/surfaceBooleanFeatures/CGALIndexedPolyhedron.H
new file mode 100644
index 0000000000000000000000000000000000000000..0eba1bee2f3511207ae228de722e73f686e9c78b
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/CGALIndexedPolyhedron.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+Typedefs
+ IndexedPolyhedron
+
+Description
+ CGAL data structures used for triSurface handling
+
+ Define CGAL_INEXACT to use Exact_predicates_inexact_constructions kernel
+ otherwise the more robust but much less efficient
+ Exact_predicates_exact_constructions will be used.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IndexedPolyhedron_H
+#define IndexedPolyhedron_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Include uint.H before CGAL headers to define __STDC_LIMIT_MACROS
+#include "uint.H"
+
+#include "CGAL3DKernel.H"
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Nef_polyhedron_3.h>
+
+#include "label.H"
+
+
+typedef CGAL::Point_3<K> Point;
+typedef CGAL::Segment_3<K> Segment;
+typedef CGAL::Direction_3<K> Direction;
+typedef CGAL::Plane_3<K> Plane;
+typedef CGAL::Triangle_3<K> Triangle;
+
+// Define new class with color and define the polyhedron types
+template<class Refs>
+struct IndexedFace
+:
+ public CGAL::HalfedgeDS_face_base<Refs>
+{
+ Foam::label index;
+ Foam::label region;
+};
+struct My_items
+:
+ public CGAL::Polyhedron_items_3
+{
+ template<class Refs, class Traits>
+ struct Face_wrapper
+ {
+ typedef IndexedFace<Refs> Face;
+ };
+};
+
+
+//typedef CGAL::Polyhedron_3<K> Polyhedron;
+typedef CGAL::Polyhedron_3<K, My_items> Polyhedron;
+
+typedef Polyhedron::HalfedgeDS HalfedgeDS;
+typedef Polyhedron::Edge_iterator Edge_iterator;
+typedef Polyhedron::Vertex Vertex;
+typedef Polyhedron::Vertex_iterator Vertex_iterator;
+typedef Polyhedron::Halfedge_handle Halfedge_handle;
+typedef Polyhedron::Edge_iterator Edge_iterator;
+typedef Polyhedron::Facet_iterator Facet_iterator;
+typedef Polyhedron::Halfedge_around_facet_const_circulator HFCC;
+typedef Polyhedron::Vertex_const_iterator VCI;
+
+typedef CGAL::Nef_polyhedron_3<K> Nef_polyhedron;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/Make/options b/applications/utilities/surface/surfaceBooleanFeatures/Make/options
index 727dbd1c2b0081d68cb522eaf36925d415c5495a..ed05fef039a9070d0287ea2f6c2167acbc15427f 100644
--- a/applications/utilities/surface/surfaceBooleanFeatures/Make/options
+++ b/applications/utilities/surface/surfaceBooleanFeatures/Make/options
@@ -1,9 +1,29 @@
+EXE_NDEBUG = -DNDEBUG
+/* EXE_NDEBUG = -g -O0 -DFULLDEBUG */
+
+
+c++CGALWARN = -Wno-old-style-cast
+
+/*-- Define NO_CGAL to avoid using CGAL altogether */
+/*-- Define CGAL_INEXACT to use inexact CGAL constructions */
+
+include $(GENERAL_RULES)/CGAL
+
EXE_INC = \
- -I$(LIB_SRC)/triSurface/lnInclude \
+ ${ROUNDING_MATH} \
+ ${EXE_NDEBUG} \
+ ${CGAL_INC} \
+ ${c++CGALWARN} \
+ $(COMPILE_FLAGS) \
+ -IPolyhedronReader \
+ -I$(FOAM_SRC)/surfMesh/lnInclude \
+ -I$(FOAM_SRC)/triSurface/lnInclude \
-I$(LIB_SRC)/edgeMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
+ -lsurfMesh \
-ltriSurface \
-ledgeMesh \
- -lmeshTools
+ -lmeshTools \
+ $(LINK_FLAGS)
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/files b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..d0e01d2747a4c42de80f72ea5f2d0b8ab62cb4fc
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/files
@@ -0,0 +1,3 @@
+PolyhedronReader.C
+
+LIB = $(FOAM_LIBBIN)/libPolyhedronReader
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/options b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..6603adda0ff4df4762f57ce04094904a0f32664b
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/Make/options
@@ -0,0 +1,25 @@
+EXE_NDEBUG = -DNDEBUG
+/* EXE_NDEBUG = -g -O0 -DFULLDEBUG */
+
+
+c++CGALWARN = -Wno-old-style-cast
+
+/*-- Define CGAL_INEXACT to use inexact CGAL constructions */
+
+include $(GENERAL_RULES)/CGAL
+
+EXE_INC = \
+ ${ROUNDING_MATH} \
+ ${EXE_NDEBUG} \
+ ${CGAL_INC} \
+ ${c++CGALWARN} \
+ -I.. \
+ -I$(FOAM_SRC)/surfMesh/lnInclude \
+ -I$(FOAM_SRC)/triSurface/lnInclude \
+ -I$(LIB_SRC)/edgeMesh/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I/usr/include/Qt
+
+LIB_LIBS = \
+ -L${CGAL_ARCH_PATH}/lib \
+ -ltriSurface
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.C b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.C
new file mode 100644
index 0000000000000000000000000000000000000000..462b9a6f05ab278d1cc76188cf930b6bb6a94bda
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.C
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PolyhedronReader.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::PolyhedronReader::PolyhedronReader(const triSurface& s, Polyhedron& p)
+{
+ Build_triangle<HalfedgeDS> triangle(s);
+ p.delegate(triangle);
+ // Populate index and region
+ Foam::label nTris = 0;
+ for
+ (
+ Facet_iterator fi = p.facets_begin();
+ fi != p.facets_end();
+ ++fi
+ )
+ {
+ fi->index = nTris++;
+ fi->region = s[fi->index].region();
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.H b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.H
new file mode 100644
index 0000000000000000000000000000000000000000..ef3e88a0ac3ecce27e3e27306807c80754d411df
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReader.H
@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::PolyhedronReader
+
+Description
+
+SourceFiles
+ PolyhedronReader.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PolyhedronReader_H
+#define PolyhedronReader_H
+
+#include "CGALIndexedPolyhedron.H"
+#include "triSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class PolyhedronReader Declaration
+\*---------------------------------------------------------------------------*/
+
+class PolyhedronReader
+{
+ // Private Classes
+
+ template<class HDS>
+ class Build_triangle
+ :
+ public CGAL::Modifier_base<HDS>
+ {
+ const triSurface& s_;
+
+ public:
+
+ Build_triangle(const triSurface& s);
+
+ void operator()(HDS& hds);
+ };
+
+
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ PolyhedronReader(const PolyhedronReader&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const PolyhedronReader&);
+
+
+public:
+
+ // Constructors
+
+ //- Populate polyhedron from surface
+ PolyhedronReader(const triSurface& s, Polyhedron& p);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+# include "PolyhedronReaderTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReaderTemplates.C b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReaderTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..53ea3b1b87b5349f578cb7f622121ad1404e0ac6
--- /dev/null
+++ b/applications/utilities/surface/surfaceBooleanFeatures/PolyhedronReader/PolyhedronReaderTemplates.C
@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PolyhedronReader.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class HDS>
+Foam::PolyhedronReader::Build_triangle<HDS>::Build_triangle
+(
+ const triSurface& s
+)
+:
+ s_(s)
+{}
+
+
+// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
+
+template<class HDS>
+void Foam::PolyhedronReader::Build_triangle<HDS>::operator()(HDS& hds)
+{
+ // Postcondition: hds is a valid polyhedral surface.
+ CGAL::Polyhedron_incremental_builder_3<HDS> B(hds, true);
+
+ B.begin_surface(s_.nPoints(), s_.size());
+
+ typedef typename HDS::Vertex Vertex;
+ typedef typename Vertex::Point Point;
+
+ const Foam::pointField& pts = s_.points();
+ forAll(pts, i)
+ {
+ const Foam::point& pt = pts[i];
+ B.add_vertex(Point(pt.x(), pt.y(), pt.z()));
+ }
+ forAll(s_, i)
+ {
+ const Foam::labelledTri& t = s_[i];
+ B.begin_facet();
+ B.add_vertex_to_facet(t[0]);
+ B.add_vertex_to_facet(t[1]);
+ B.add_vertex_to_facet(t[2]);
+ B.end_facet();
+ }
+ B.end_surface();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
index 6c4e67dd7d14bf423c146ec0dc4a60909e15041c..e948cd1fdbc59b7a7dfc48eac2fe7d86bc1d94e4 100644
--- a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
+++ b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -74,11 +74,34 @@ Description
#include "featureEdgeMesh.H"
#include "extendedFeatureEdgeMesh.H"
#include "triSurfaceSearch.H"
+#include "triSurfaceMesh.H"
#include "OFstream.H"
+#include "OBJstream.H"
#include "booleanSurface.H"
#include "edgeIntersections.H"
#include "meshTools.H"
-#include "labelPair.H"
+#include "DynamicField.H"
+
+
+#ifndef NO_CGAL
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include "CGALIndexedPolyhedron.H"
+#include "PolyhedronReader.H"
+typedef CGAL::AABB_face_graph_triangle_primitive
+<
+ Polyhedron, CGAL::Default, CGAL::Tag_false
+> Primitive;
+typedef CGAL::AABB_traits<K, Primitive> Traits;
+typedef CGAL::AABB_tree<Traits> Tree;
+
+typedef boost::optional<Tree::Intersection_and_primitive_id<Segment>::Type >
+Segment_intersection;
+
+#endif // NO_CGAL
+
using namespace Foam;
@@ -105,7 +128,7 @@ bool intersectSurfaces
scalarField surfPointTol
(
- 1e-3*edgeIntersections::minEdgeLength(surf)
+ max(1e-3*edgeIntersections::minEdgeLength(surf), SMALL)
);
// Determine raw intersections
@@ -178,7 +201,7 @@ bool intersectSurfaces
scalarField surf1PointTol
(
- 1e-3*edgeIntersections::minEdgeLength(surf1)
+ max(1e-3*edgeIntersections::minEdgeLength(surf1), SMALL)
);
// Determine raw intersections
@@ -221,7 +244,7 @@ bool intersectSurfaces
scalarField surf2PointTol
(
- 1e-3*edgeIntersections::minEdgeLength(surf2)
+ max(1e-3*edgeIntersections::minEdgeLength(surf2), SMALL)
);
// Determine raw intersections
@@ -311,8 +334,8 @@ void calcEdgeCuts
triSurface& surf1,
triSurface& surf2,
const bool perturb,
- edgeIntersections& edge1Cuts,
- edgeIntersections& edge2Cuts
+ edgeIntersections& edgeCuts1,
+ edgeIntersections& edgeCuts2
)
{
if (perturb)
@@ -320,9 +343,9 @@ void calcEdgeCuts
intersectSurfaces
(
surf1,
- edge1Cuts,
+ edgeCuts1,
surf2,
- edge2Cuts
+ edgeCuts2
);
}
else
@@ -332,11 +355,11 @@ void calcEdgeCuts
Info<< "Determining intersections of surf1 edges with surf2 faces"
<< endl;
- edge1Cuts = edgeIntersections
+ edgeCuts1 = edgeIntersections
(
surf1,
querySurf2,
- 1e-3*edgeIntersections::minEdgeLength(surf1)
+ max(1e-3*edgeIntersections::minEdgeLength(surf1), SMALL)
);
triSurfaceSearch querySurf1(surf1);
@@ -344,308 +367,1026 @@ void calcEdgeCuts
Info<< "Determining intersections of surf2 edges with surf1 faces"
<< endl;
- edge2Cuts = edgeIntersections
+ edgeCuts2 = edgeIntersections
(
surf2,
querySurf1,
- 1e-3*edgeIntersections::minEdgeLength(surf2)
+ max(1e-3*edgeIntersections::minEdgeLength(surf2), SMALL)
);
}
}
-void calcFeaturePoints(const pointField& points, const edgeList& edges)
-{
- edgeMesh eMesh(points, edges);
+// CGAL variants
- const labelListList& pointEdges = eMesh.pointEdges();
+#ifndef NO_CGAL
- // Get total number of feature points
- label nFeaturePoints = 0;
- forAll(pointEdges, pI)
- {
- const labelList& pEdges = pointEdges[pI];
+void visitPointRegion
+(
+ const triSurface& s,
+ const label zoneI,
+ const label pointI,
+ const label startEdgeI,
+ const label startFaceI,
+ labelList& pFacesZone
+)
+{
+ const labelList& eFaces = s.edgeFaces()[startEdgeI];
- if (pEdges.size() == 1)
+ if (eFaces.size() == 2)
+ {
+ label nextFaceI;
+ if (eFaces[0] == startFaceI)
{
- nFeaturePoints++;
+ nextFaceI = eFaces[1];
}
- }
-
-
- // Calculate addressing from feature point to cut point and cut edge
- labelList featurePointToCutPoint(nFeaturePoints);
- labelList featurePointToCutEdge(nFeaturePoints);
-
- label nFeatPts = 0;
- forAll(pointEdges, pI)
- {
- const labelList& pEdges = pointEdges[pI];
-
- if (pEdges.size() == 1)
+ else if (eFaces[1] == startFaceI)
{
- featurePointToCutPoint[nFeatPts] = pI;
- featurePointToCutEdge[nFeatPts] = pEdges[0];
- nFeatPts++;
+ nextFaceI = eFaces[0];
+ }
+ else
+ {
+ FatalErrorIn("visitPointRegion(..)")
+ << "problem" << exit(FatalError);
+ nextFaceI = -1;
}
- }
-}
-int main(int argc, char *argv[])
-{
- argList::noParallel();
- argList::validArgs.append("action");
- argList::validArgs.append("surface file");
- argList::validArgs.append("surface file");
- argList::addBoolOption
- (
- "surf1Baffle",
- "Mark surface 1 as a baffle"
- );
+ label index = findIndex(s.pointFaces()[pointI], nextFaceI);
- argList::addBoolOption
- (
- "surf2Baffle",
- "Mark surface 2 as a baffle"
- );
+ if (pFacesZone[index] == -1)
+ {
+ // Mark face as been visited.
+ pFacesZone[index] = zoneI;
- argList::addBoolOption
- (
- "perturb",
- "Perturb surface points to escape degenerate intersections"
- );
+ // Step to next edge on face which is still using pointI
+ const labelList& fEdges = s.faceEdges()[nextFaceI];
- argList::addBoolOption
- (
- "invertedSpace",
- "do the surfaces have inverted space orientation, "
- "i.e. a point at infinity is considered inside. "
- "This is only sensible for union and intersection."
- );
+ label nextEdgeI = -1;
- #include "setRootCase.H"
- #include "createTime.H"
+ forAll(fEdges, i)
+ {
+ label edgeI = fEdges[i];
+ const edge& e = s.edges()[edgeI];
- word action(args.args()[1]);
+ if (edgeI != startEdgeI && (e[0] == pointI || e[1] == pointI))
+ {
+ nextEdgeI = edgeI;
- HashTable<booleanSurface::booleanOpType> validActions;
- validActions.insert("intersection", booleanSurface::INTERSECTION);
- validActions.insert("union", booleanSurface::UNION);
- validActions.insert("difference", booleanSurface::DIFFERENCE);
+ break;
+ }
+ }
- if (!validActions.found(action))
- {
- FatalErrorIn(args.executable())
- << "Unsupported action " << action << endl
- << "Supported actions:" << validActions.toc() << abort(FatalError);
- }
+ if (nextEdgeI == -1)
+ {
+ FatalErrorIn("visitPointRegion()")
+ << "Problem: cannot find edge out of " << fEdges
+ << "on face " << nextFaceI << " that uses point " << pointI
+ << " and is not edge " << startEdgeI << abort(FatalError);
+ }
- fileName surf1Name(args.args()[2]);
- Info<< "Reading surface " << surf1Name << endl;
- triSurface surf1(surf1Name);
- Info<< surf1Name << " statistics:" << endl;
- surf1.writeStats(Info);
- Info<< endl;
+ visitPointRegion
+ (
+ s,
+ zoneI,
+ pointI,
+ nextEdgeI,
+ nextFaceI,
+ pFacesZone
+ );
+ }
+ }
+}
- fileName surf2Name(args[3]);
- Info<< "Reading surface " << surf2Name << endl;
- triSurface surf2(surf2Name);
- Info<< surf2Name << " statistics:" << endl;
- surf2.writeStats(Info);
- Info<< endl;
+label dupNonManifoldPoints(triSurface& s, labelList& pointMap)
+{
+ const labelListList& pf = s.pointFaces();
+ const labelListList& fe = s.faceEdges();
+ const edgeList& edges = s.edges();
- const bool surf1Baffle = args.optionFound("surf1Baffle");
- const bool surf2Baffle = args.optionFound("surf2Baffle");
- edgeIntersections edge1Cuts;
- edgeIntersections edge2Cuts;
+ DynamicField<point> newPoints(s.points());
+ // From dupSurf back to s.pointa
+ DynamicList<label> newPointMap(identity(newPoints.size()));
+ List<labelledTri> newFaces(s);
+ label nNonManifold = 0;
- bool invertedSpace = args.optionFound("invertedSpace");
- if (invertedSpace && validActions[action] == booleanSurface::DIFFERENCE)
+ forAll(pf, pointI)
{
- FatalErrorIn(args.executable())
- << "Inverted space only makes sense for union or intersection."
- << exit(FatalError);
- }
-
- // Calculate the points where the edges are cut by the other surface
- calcEdgeCuts
- (
- surf1,
- surf2,
- args.optionFound("perturb"),
- edge1Cuts,
- edge2Cuts
- );
+ const labelList& pFaces = pf[pointI];
- // Determine intersection edges from the edge cuts
- surfaceIntersection inter
- (
- surf1,
- edge1Cuts,
- surf2,
- edge2Cuts
- );
+ // Visited faces (as indices into pFaces)
+ labelList pFacesZone(pFaces.size(), -1);
- fileName sFeatFileName
- (
- surf1Name.lessExt().name()
- + "_"
- + surf2Name.lessExt().name()
- + "_"
- + action
- );
+ label nZones = 0;
+ label index = 0;
+ for (; index < pFacesZone.size(); index++)
+ {
+ if (pFacesZone[index] == -1)
+ {
+ label zoneI = nZones++;
+ pFacesZone[index] = zoneI;
- label nFeatEds = inter.cutEdges().size();
+ label faceI = pFaces[index];
+ const labelList& fEdges = fe[faceI];
- DynamicList<vector> normals(2*nFeatEds);
- vectorField edgeDirections(nFeatEds, vector::zero);
- DynamicList<extendedFeatureEdgeMesh::sideVolumeType> normalVolumeTypes
- (
- 2*nFeatEds
- );
- List<DynamicList<label> > edgeNormals(nFeatEds);
- List<DynamicList<label> > normalDirections(nFeatEds);
+ // Find starting edge
+ forAll(fEdges, fEdgeI)
+ {
+ label edgeI = fEdges[fEdgeI];
+ const edge& e = edges[edgeI];
+
+ if (e[0] == pointI || e[1] == pointI)
+ {
+ visitPointRegion
+ (
+ s,
+ zoneI,
+ pointI,
+ edgeI,
+ faceI,
+ pFacesZone
+ );
+ }
+ }
+ }
+ }
- forAllConstIter(labelPairLookup, inter.facePairToEdge(), iter)
- {
- const label& cutEdgeI = iter();
- const labelPair& facePair = iter.key();
- const edge& fE = inter.cutEdges()[cutEdgeI];
+ // Subset
+ if (nZones > 1)
+ {
+ for (label zoneI = 1; zoneI < nZones; zoneI++)
+ {
+ label newPointI = newPoints.size();
+ newPointMap.append(s.meshPoints()[pointI]);
+ newPoints.append(s.points()[s.meshPoints()[pointI]]);
- const vector& norm1 = surf1.faceNormals()[facePair.first()];
- const vector& norm2 = surf2.faceNormals()[facePair.second()];
+ forAll(pFacesZone, index)
+ {
+ if (pFacesZone[index] == zoneI)
+ {
+ label faceI = pFaces[index];
+ const labelledTri& localF = s.localFaces()[faceI];
+ forAll(localF, fp)
+ {
+ if (localF[fp] == pointI)
+ {
+ newFaces[faceI][fp] = newPointI;
+ }
+ }
+ }
+ }
+ }
+ nNonManifold++;
+ }
+ }
- DynamicList<label>& eNormals = edgeNormals[cutEdgeI];
- DynamicList<label>& nDirections = normalDirections[cutEdgeI];
- edgeDirections[cutEdgeI] = fE.vec(inter.cutPoints());
+ Info<< "Duplicating " << nNonManifold << " points out of " << s.nPoints()
+ << endl;
+ if (nNonManifold > 0)
+ {
+ triSurface dupSurf(newFaces, s.patches(), newPoints, true);
- normals.append(norm1);
- eNormals.append(normals.size() - 1);
+ // Create map from dupSurf localPoints to s.localPoints
+ const Map<label> mpm = s.meshPointMap();
- if (surf1Baffle)
- {
- normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
+ const labelList& dupMp = dupSurf.meshPoints();
- nDirections.append(1);
- }
- else
+ labelList dupPointMap(dupMp.size());
+ forAll(dupMp, pointI)
{
- normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
- nDirections.append
- (
- calcNormalDirection
- (
- norm1,
- norm2,
- edgeDirections[cutEdgeI],
- surf1[facePair.first()].centre(surf1.points()),
- inter.cutPoints()[fE.start()]
- )
- );
+ label dupMeshPointI = dupMp[pointI];
+ label meshPointI = newPointMap[dupMeshPointI];
+ dupPointMap[pointI] = mpm[meshPointI];
}
- normals.append(norm2);
- eNormals.append(normals.size() - 1);
- if (surf2Baffle)
+ forAll(dupPointMap, pointI)
{
- normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
+ const point& dupPt = dupSurf.points()[dupMp[pointI]];
+ label sLocalPointI = dupPointMap[pointI];
+ label sMeshPointI = s.meshPoints()[sLocalPointI];
+ const point& sPt = s.points()[sMeshPointI];
- nDirections.append(1);
+ if (mag(dupPt-sPt) > SMALL)
+ {
+ FatalErrorIn("dupNonManifoldPoints(..)")
+ << "dupPt:" << dupPt
+ << " sPt:" << sPt
+ << exit(FatalError);
+ }
}
- else
- {
- normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
- nDirections.append
- (
- calcNormalDirection
- (
- norm2,
- norm1,
- edgeDirections[cutEdgeI],
- surf2[facePair.second()].centre(surf2.points()),
- inter.cutPoints()[fE.start()]
- )
- );
- }
+ //s.transfer(dupSurf);
+ s = dupSurf;
+ pointMap = UIndirectList<label>(pointMap, dupPointMap)();
+ }
- if (surf1Baffle)
- {
- normals.append(norm2);
+ return nNonManifold;
+}
- if (surf2Baffle)
- {
- normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
- nDirections.append(1);
- }
- else
+// Find intersections of surf1 by edges of surf2. Return number of degenerate
+// cuts (cuts through faces/edges/points)
+labelPair edgeIntersectionsCGAL
+(
+ const Tree& tree,
+ const triSurface& surf,
+ const pointField& points,
+ edgeIntersections& edgeCuts
+)
+{
+ const edgeList& edges = surf.edges();
+ const labelList& meshPoints = surf.meshPoints();
+
+ //Info<< "Intersecting CGAL surface ..." << endl;
+ List<List<pointIndexHit> > intersections(edges.size());
+ labelListList classifications(edges.size());
+
+ label nPoints = 0;
+ label nSegments = 0;
+
+ std::vector<Segment_intersection> segments;
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+
+ const point& a = points[meshPoints[e[0]]];
+ const point& b = points[meshPoints[e[1]]];
+
+ K::Segment_3 segment_query
+ (
+ Point(a.x(), a.y(), a.z()),
+ Point(b.x(), b.y(), b.z())
+ );
+
+ segments.clear();
+ tree.all_intersections(segment_query, std::back_inserter(segments));
+
+ for
+ (
+ std::vector<Segment_intersection>::const_iterator iter =
+ segments.begin(),
+ end = segments.end();
+ iter != end;
+ ++iter
+ )
+ {
+ // Get intersection object
+ if (const Point* ptPtr = boost::get<Point>(&((*iter)->first)))
{
- normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
+ point pt
+ (
+ CGAL::to_double(ptPtr->x()),
+ CGAL::to_double(ptPtr->y()),
+ CGAL::to_double(ptPtr->z())
+ );
- nDirections.append
+ Polyhedron::Face_handle f = (*iter)->second;
+
+ intersections[edgeI].append
(
- calcNormalDirection
+ pointIndexHit
(
- norm2,
- norm1,
- edgeDirections[cutEdgeI],
- surf2[facePair.second()].centre(surf2.points()),
- inter.cutPoints()[fE.start()]
+ true,
+ pt,
+ f->index
)
);
+ // Intersection on edge interior
+ classifications[edgeI].append(-1);
+ nPoints++;
}
+ else if
+ (
+ const Segment* sPtr = boost::get<Segment>(&((*iter)->first))
+ )
+ {
+ //std::cout
+ // << "intersection object is a segment:" << sPtr->source()
+ // << " " << sPtr->target() << std::endl;
- eNormals.append(normals.size() - 1);
- }
-
- if (surf2Baffle)
- {
- normals.append(norm1);
+ Polyhedron::Face_handle f = (*iter)->second;
+ //std::cout<< "triangle:" << f->index
+ // << " region:" << f->region << std::endl;
- if (surf1Baffle)
- {
- normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
+ const point source
+ (
+ CGAL::to_double(sPtr->source().x()),
+ CGAL::to_double(sPtr->source().y()),
+ CGAL::to_double(sPtr->source().z())
+ );
- nDirections.append(1);
- }
- else
- {
- normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
+ const point target
+ (
+ CGAL::to_double(sPtr->target().x()),
+ CGAL::to_double(sPtr->target().y()),
+ CGAL::to_double(sPtr->target().z())
+ );
- nDirections.append
+ // Make up some intersection point
+ intersections[edgeI].append
(
- calcNormalDirection
+ pointIndexHit
(
- norm1,
- norm2,
- edgeDirections[cutEdgeI],
- surf1[facePair.first()].centre(surf1.points()),
- inter.cutPoints()[fE.start()]
+ true,
+ 0.5*(source+target),
+ f->index
)
);
+ // Intersection aligned with face. Tbd: enums
+ classifications[edgeI].append(2);
+ nSegments++;
}
-
- eNormals.append(normals.size() - 1);
}
}
+ edgeCuts = edgeIntersections(intersections, classifications);
- label internalStart = -1;
- label nIntOrExt = 0;
- label nFlat = 0;
- label nOpen = 0;
- label nMultiple = 0;
+ return labelPair(nPoints, nSegments);
+}
+
+
+// Intersect edges of surf1 until no more degenerate intersections. Return
+// number of degenerates
+labelPair edgeIntersectionsAndShuffleCGAL
+(
+ Random& rndGen,
+ const triSurface& surf2,
+ const scalarField& surf1PointTol,
+ triSurface& surf1,
+ edgeIntersections& edgeCuts1
+)
+{
+ //Info<< "Constructing CGAL surface ..." << endl;
+ Polyhedron p;
+ PolyhedronReader(surf2, p);
+
+
+ //Info<< "Constructing CGAL tree ..." << endl;
+ const Tree tree(p.facets_begin(), p.facets_end(), p);
+
+
+ labelPair cuts(0, 0);
+
+ // Update surface1 points until no longer intersecting
+ pointField surf1Points(surf1.points());
+
+ for (label iter = 0; iter < 5; iter++)
+ {
+ // See which edges of 1 intersect 2
+ Info<< "Determining intersections of " << surf1.nEdges()
+ << " edges with surface of " << label(tree.size()) << " triangles"
+ << endl;
+ cuts = edgeIntersectionsCGAL
+ (
+ tree,
+ surf1,
+ surf1Points,
+ edgeCuts1
+ );
+ Info<< "Determined intersections:" << nl
+ << " edges : " << surf1.nEdges() << nl
+ << " non-degenerate cuts : " << cuts.first() << nl
+ << " degenerate cuts : " << cuts.second() << nl
+ << endl;
+
+ if (cuts.second() == 0)
+ {
+ break;
+ }
+
+ Info<< "Shuffling conflicting points" << endl;
+
+ const labelListList& edgeStat = edgeCuts1.classification();
+ const edgeList& edges = surf1.edges();
+ const labelList& mp = surf1.meshPoints();
+ const point p05(0.5, 0.5, 0.5);
+
+ forAll(edgeStat, edgeI)
+ {
+ const labelList& stat = edgeStat[edgeI];
+ forAll(stat, i)
+ {
+ if (stat[i] == 2)
+ {
+ const edge& e = edges[edgeI];
+ forAll(e, eI)
+ {
+ vector d = rndGen.vector01()-p05;
+ surf1Points[mp[e[eI]]] += surf1PointTol[e[eI]]*d;
+ }
+ }
+ }
+ }
+ }
+ surf1.movePoints(surf1Points);
+
+ return cuts;
+}
+
+
+// Return map from subSurf.edges() back to surf.edges()
+labelList matchEdges
+(
+ const triSurface& surf,
+ const triSurface& subSurf,
+ const labelList& pointMap
+)
+{
+ if (pointMap.size() != subSurf.nPoints())
+ {
+ FatalErrorIn("findEdges(..)") << "problem" << exit(FatalError);
+ }
+
+ labelList edgeMap(subSurf.nEdges(), -1);
+
+ const edgeList& edges = surf.edges();
+ const labelListList& pointEdges = surf.pointEdges();
+
+ const edgeList& subEdges = subSurf.edges();
+
+
+ forAll(subEdges, subEdgeI)
+ {
+ const edge& subE = subEdges[subEdgeI];
+
+ // Match points on edge to those on surf
+ label start = pointMap[subE[0]];
+ label end = pointMap[subE[1]];
+ const labelList& pEdges = pointEdges[start];
+ forAll(pEdges, pEdgeI)
+ {
+ label edgeI = pEdges[pEdgeI];
+ const edge& e = edges[edgeI];
+
+ if (e.otherVertex(start) == end)
+ {
+ if (edgeMap[subEdgeI] == -1)
+ {
+ edgeMap[subEdgeI] = edgeI;
+ }
+ else if (edgeMap[subEdgeI] != edgeI)
+ {
+ WarningIn("findEdges(..)") << "sub edge "
+ << subE << " points:"
+ << subE.line(subSurf.localPoints())
+ << " matches to " << edgeI
+ << " points:" << e.line(surf.localPoints())
+ << " but also " << edgeMap[subEdgeI]
+ << " points:"
+ << edges[edgeMap[subEdgeI]].line(surf.localPoints())
+ << exit(FatalError);
+ }
+ break;
+ }
+ }
+
+ if (edgeMap[subEdgeI] == -1)
+ {
+ FatalErrorIn("findEdges(..)") << "did not find edge matching "
+ << subE << " at:" << subSurf.localPoints()[subE[0]]
+ << subSurf.localPoints()[subE[1]]
+ << exit(FatalError);
+ }
+ }
+
+ return edgeMap;
+}
+
+
+void calcEdgeCutsCGAL
+(
+ triSurface& surf1,
+ triSurface& surf2,
+ const bool perturb,
+ edgeIntersections& edgeCuts1,
+ edgeIntersections& edgeCuts2
+)
+{
+ if (!perturb)
+ {
+ // See which edges of 1 intersect 2
+ {
+ Info<< "Constructing CGAL surface ..." << endl;
+ Polyhedron p;
+ PolyhedronReader(surf2, p);
+
+ Info<< "Constructing CGAL tree ..." << endl;
+ const Tree tree(p.facets_begin(), p.facets_end(), p);
+
+ edgeIntersectionsCGAL
+ (
+ tree,
+ surf1,
+ surf1.points(),
+ edgeCuts1
+ );
+ }
+ // See which edges of 2 intersect 1
+ {
+ Info<< "Constructing CGAL surface ..." << endl;
+ Polyhedron p;
+ PolyhedronReader(surf1, p);
+
+ Info<< "Constructing CGAL tree ..." << endl;
+ const Tree tree(p.facets_begin(), p.facets_end(), p);
+
+ edgeIntersectionsCGAL
+ (
+ tree,
+ surf2,
+ surf2.points(),
+ edgeCuts2
+ );
+ }
+ }
+ else
+ {
+ const scalarField surf1PointTol
+ (
+ max(1e-8*edgeIntersections::minEdgeLength(surf1), SMALL)
+ );
+ const scalarField surf2PointTol
+ (
+ max(1e-8*edgeIntersections::minEdgeLength(surf2), SMALL)
+ );
+
+
+ Random rndGen(0);
+
+ labelPair cuts1;
+ labelPair cuts2;
+
+ for (label iter = 0; iter < 10; iter++)
+ {
+ // Find intersections of surf1 edges with surf2 triangles
+ cuts1 = edgeIntersectionsAndShuffleCGAL
+ (
+ rndGen,
+ surf2,
+ surf1PointTol,
+ surf1,
+ edgeCuts1
+ );
+
+ // Find intersections of surf2 edges with surf1 triangles
+ cuts2 = edgeIntersectionsAndShuffleCGAL
+ (
+ rndGen,
+ surf1,
+ surf2PointTol,
+ surf2,
+ edgeCuts2
+ );
+
+ if (cuts1.second() + cuts2.second() == 0)
+ {
+ break;
+ }
+ }
+ }
+}
+
+
+void calcEdgeCutsBitsCGAL
+(
+ triSurface& surf1,
+ triSurface& surf2,
+ const bool perturb,
+ edgeIntersections& edgeCuts1,
+ edgeIntersections& edgeCuts2
+)
+{
+ label nZones1 = 1;
+ labelList zone1;
+ {
+ labelHashSet orientationEdge(surf1.size()/1000);
+ PatchTools::checkOrientation(surf1, false, &orientationEdge);
+ nZones1 = PatchTools::markZones(surf1, orientationEdge, zone1);
+
+ Info<< "Surface triangles " << surf1.size()
+ << " number of zones (orientation or non-manifold):"
+ << nZones1 << endl;
+ }
+
+ label nZones2 = 1;
+ labelList zone2;
+ {
+ labelHashSet orientationEdge(surf2.size()/1000);
+ PatchTools::checkOrientation(surf2, false, &orientationEdge);
+ nZones2 = PatchTools::markZones(surf2, orientationEdge, zone2);
+
+ Info<< "Surface triangles " << surf2.size()
+ << " number of zones (orientation or non-manifold):"
+ << nZones2 << endl;
+ }
+
+
+ if (nZones1 == 1 && nZones2 == 1)
+ {
+ calcEdgeCutsCGAL
+ (
+ surf1,
+ surf2,
+ perturb,
+ edgeCuts1,
+ edgeCuts2
+ );
+ }
+ else
+ {
+ edgeCuts1 = edgeIntersections
+ (
+ List<List<pointIndexHit> >(surf1.nEdges()),
+ labelListList(surf1.nEdges())
+ );
+ edgeCuts2 = edgeIntersections
+ (
+ List<List<pointIndexHit> >(surf2.nEdges()),
+ labelListList(surf2.nEdges())
+ );
+
+
+ for (label zone1I = 0; zone1I < nZones1; zone1I++)
+ {
+ // Generate sub surface for zone1I
+ boolList includeMap1(surf1.size(), false);
+
+ forAll(zone1, faceI)
+ {
+ if (zone1[faceI] == zone1I)
+ {
+ includeMap1[faceI] = true;
+ }
+ }
+
+ // Subset. Map from local points on subset to local points on
+ // original
+ labelList pointMap1;
+ labelList faceMap1;
+ triSurface subSurf1
+ (
+ surf1.subsetMesh
+ (
+ includeMap1,
+ pointMap1,
+ faceMap1
+ )
+ );
+
+ // Split non-manifold points; update pointMap
+ dupNonManifoldPoints(subSurf1, pointMap1);
+
+ const boundBox subBb1
+ (
+ subSurf1.points(),
+ subSurf1.meshPoints(),
+ false
+ );
+
+ const labelList edgeMap1
+ (
+ matchEdges
+ (
+ surf1,
+ subSurf1,
+ pointMap1
+ )
+ );
+
+
+ for (label zone2I = 0; zone2I < nZones2; zone2I++)
+ {
+ // Generate sub surface for zone2I
+ boolList includeMap2(surf2.size(), false);
+
+ forAll(zone2, faceI)
+ {
+ if (zone2[faceI] == zone2I)
+ {
+ includeMap2[faceI] = true;
+ }
+ }
+
+ labelList pointMap2;
+ labelList faceMap2;
+ triSurface subSurf2
+ (
+ surf2.subsetMesh
+ (
+ includeMap2,
+ pointMap2,
+ faceMap2
+ )
+ );
+
+
+ const boundBox subBb2
+ (
+ subSurf2.points(),
+ subSurf2.meshPoints(),
+ false
+ );
+
+ // Short-circuit expensive calculations
+ if (!subBb2.overlaps(subBb1))
+ {
+ continue;
+ }
+
+
+ // Split non-manifold points; update pointMap
+ dupNonManifoldPoints(subSurf2, pointMap2);
+
+ const labelList edgeMap2
+ (
+ matchEdges
+ (
+ surf2,
+ subSurf2,
+ pointMap2
+ )
+ );
+
+
+ // Do cuts
+ edgeIntersections subEdgeCuts1;
+ edgeIntersections subEdgeCuts2;
+ calcEdgeCutsCGAL
+ (
+ subSurf1,
+ subSurf2,
+ perturb,
+ subEdgeCuts1,
+ subEdgeCuts2
+ );
+
+ // Move original surface
+ {
+ pointField points2(surf2.points());
+ forAll(pointMap2, i)
+ {
+ label subMeshPointI = subSurf2.meshPoints()[i];
+ label meshPointI = surf2.meshPoints()[pointMap2[i]];
+ points2[meshPointI] = subSurf2.points()[subMeshPointI];
+ }
+ surf2.movePoints(points2);
+ }
+
+ // Insert into main structure
+ edgeCuts1.merge(subEdgeCuts1, edgeMap1, faceMap2);
+ edgeCuts2.merge(subEdgeCuts2, edgeMap2, faceMap1);
+ }
+
+
+ // Move original surface
+ {
+ pointField points1(surf1.points());
+ forAll(pointMap1, i)
+ {
+ label subMeshPointI = subSurf1.meshPoints()[i];
+ label meshPointI = surf1.meshPoints()[pointMap1[i]];
+ points1[meshPointI] = subSurf1.points()[subMeshPointI];
+ }
+ surf1.movePoints(points1);
+ }
+ }
+ }
+}
+
+
+#endif // NO_CGAL
+
+
+//void calcFeaturePoints(const pointField& points, const edgeList& edges)
+//{
+// edgeMesh eMesh(points, edges);
+//
+// const labelListList& pointEdges = eMesh.pointEdges();
+//
+//
+// // Get total number of feature points
+// label nFeaturePoints = 0;
+// forAll(pointEdges, pI)
+// {
+// const labelList& pEdges = pointEdges[pI];
+//
+// if (pEdges.size() == 1)
+// {
+// nFeaturePoints++;
+// }
+// }
+//
+//
+// // Calculate addressing from feature point to cut point and cut edge
+// labelList featurePointToCutPoint(nFeaturePoints);
+// labelList featurePointToCutEdge(nFeaturePoints);
+//
+// label nFeatPts = 0;
+// forAll(pointEdges, pI)
+// {
+// const labelList& pEdges = pointEdges[pI];
+//
+// if (pEdges.size() == 1)
+// {
+// featurePointToCutPoint[nFeatPts] = pI;
+// featurePointToCutEdge[nFeatPts] = pEdges[0];
+// nFeatPts++;
+// }
+// }
+//
+//
+//
+// label concaveStart = 0;
+// label mixedStart = 0;
+// label nonFeatureStart = nFeaturePoints;
+//
+//
+// labelListList featurePointNormals(nFeaturePoints);
+// labelListList featurePointEdges(nFeaturePoints);
+// labelList regionEdges;
+//}
+
+
+autoPtr<extendedFeatureEdgeMesh> createEdgeMesh
+(
+ const IOobject& io,
+ const booleanSurface::booleanOpType action,
+ const bool surf1Baffle,
+ const bool surf2Baffle,
+ const bool invertedSpace,
+ const triSurface& surf1,
+ const edgeIntersections& edgeCuts1,
+ const triSurface& surf2,
+ const edgeIntersections& edgeCuts2
+)
+{
+ // Determine intersection edges from the edge cuts
+ surfaceIntersection inter
+ (
+ surf1,
+ edgeCuts1,
+ surf2,
+ edgeCuts2
+ );
+
+ label nFeatEds = inter.cutEdges().size();
+
+ DynamicList<vector> normals(2*nFeatEds);
+ vectorField edgeDirections(nFeatEds, vector::zero);
+ DynamicList<extendedFeatureEdgeMesh::sideVolumeType> normalVolumeTypes
+ (
+ 2*nFeatEds
+ );
+ List<DynamicList<label> > edgeNormals(nFeatEds);
+ List<DynamicList<label> > normalDirections(nFeatEds);
+
+
+ const triSurface& s1 = surf1;
+ const triSurface& s2 = surf2;
+
+ forAllConstIter(labelPairLookup, inter.facePairToEdge(), iter)
+ {
+ const label& cutEdgeI = iter();
+ const labelPair& facePair = iter.key();
+
+ const edge& fE = inter.cutEdges()[cutEdgeI];
+
+ const vector& norm1 = s1.faceNormals()[facePair.first()];
+ const vector& norm2 = s2.faceNormals()[facePair.second()];
+
+ DynamicList<label>& eNormals = edgeNormals[cutEdgeI];
+ DynamicList<label>& nDirections = normalDirections[cutEdgeI];
+
+ edgeDirections[cutEdgeI] = fE.vec(inter.cutPoints());
+
+ normals.append(norm1);
+ eNormals.append(normals.size() - 1);
+
+ if (surf1Baffle)
+ {
+ normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
+
+ nDirections.append(1);
+ }
+ else
+ {
+ normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
+ nDirections.append
+ (
+ calcNormalDirection
+ (
+ norm1,
+ norm2,
+ edgeDirections[cutEdgeI],
+ s1[facePair.first()].centre(s1.points()),
+ inter.cutPoints()[fE.start()]
+ )
+ );
+ }
+
+ normals.append(norm2);
+ eNormals.append(normals.size() - 1);
+
+ if (surf2Baffle)
+ {
+ normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
+
+ nDirections.append(1);
+ }
+ else
+ {
+ normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
+
+ nDirections.append
+ (
+ calcNormalDirection
+ (
+ norm2,
+ norm1,
+ edgeDirections[cutEdgeI],
+ s2[facePair.second()].centre(s2.points()),
+ inter.cutPoints()[fE.start()]
+ )
+ );
+ }
+
+
+ if (surf1Baffle)
+ {
+ normals.append(norm2);
+
+ if (surf2Baffle)
+ {
+ normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
+
+ nDirections.append(1);
+ }
+ else
+ {
+ normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
+
+ nDirections.append
+ (
+ calcNormalDirection
+ (
+ norm2,
+ norm1,
+ edgeDirections[cutEdgeI],
+ s2[facePair.second()].centre(s2.points()),
+ inter.cutPoints()[fE.start()]
+ )
+ );
+ }
+
+ eNormals.append(normals.size() - 1);
+ }
+
+ if (surf2Baffle)
+ {
+ normals.append(norm1);
+
+ if (surf1Baffle)
+ {
+ normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
+
+ nDirections.append(1);
+ }
+ else
+ {
+ normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
+
+ nDirections.append
+ (
+ calcNormalDirection
+ (
+ norm1,
+ norm2,
+ edgeDirections[cutEdgeI],
+ s1[facePair.first()].centre(s1.points()),
+ inter.cutPoints()[fE.start()]
+ )
+ );
+ }
+
+ eNormals.append(normals.size() - 1);
+ }
+ }
+
+
+ label internalStart = -1;
+ label nIntOrExt = 0;
+ label nFlat = 0;
+ label nOpen = 0;
+ label nMultiple = 0;
forAll(edgeNormals, eI)
{
@@ -675,7 +1416,7 @@ int main(int argc, char *argv[])
}
}
- if (validActions[action] == booleanSurface::UNION)
+ if (action == booleanSurface::UNION)
{
if (!invertedSpace)
{
@@ -688,7 +1429,7 @@ int main(int argc, char *argv[])
internalStart = nIntOrExt;
}
}
- else if (validActions[action] == booleanSurface::INTERSECTION)
+ else if (action == booleanSurface::INTERSECTION)
{
if (!invertedSpace)
{
@@ -701,14 +1442,14 @@ int main(int argc, char *argv[])
internalStart = 0;
}
}
- else if (validActions[action] == booleanSurface::DIFFERENCE)
+ else if (action == booleanSurface::DIFFERENCE)
{
// All edges are external
internalStart = nIntOrExt;
}
else
{
- FatalErrorIn(args.executable())
+ FatalErrorIn("createEdgeMesh(..)")
<< "Unsupported booleanSurface:booleanOpType and space "
<< action << " " << invertedSpace
<< abort(FatalError);
@@ -734,42 +1475,257 @@ int main(int argc, char *argv[])
normalDirectionsTmp[i] = normalDirections[i];
}
- calcFeaturePoints(inter.cutPoints(), inter.cutEdges());
+ //calcFeaturePoints(inter.cutPoints(), inter.cutEdges());
+
+ return autoPtr<extendedFeatureEdgeMesh>
+ (
+ new extendedFeatureEdgeMesh
+ (
+ io,
+ inter.cutPoints(),
+ inter.cutEdges(),
+
+ 0, // concaveStart,
+ 0, // mixedStart,
+ 0, // nonFeatureStart,
+
+ internalStart, // internalStart,
+ nIntOrExt, // flatStart,
+ nIntOrExt + nFlat, // openStart,
+ nIntOrExt + nFlat + nOpen, // multipleStart,
+
+ normalsTmp,
+ normalVolumeTypesTmp,
+ edgeDirections,
+ normalDirectionsTmp,
+ edgeNormalsTmp,
+
+ labelListList(0), // featurePointNormals,
+ labelListList(0), // featurePointEdges,
+ labelList(0) // regionEdges
+ )
+ );
+}
+
+int main(int argc, char *argv[])
+{
+ argList::noParallel();
+ argList::validArgs.append("action");
+ argList::validArgs.append("surface file");
+ argList::validArgs.append("surface file");
- extendedFeatureEdgeMesh feMesh
+ argList::addBoolOption
+ (
+ "surf1Baffle",
+ "Mark surface 1 as a baffle"
+ );
+
+ argList::addBoolOption
+ (
+ "surf2Baffle",
+ "Mark surface 2 as a baffle"
+ );
+
+ argList::addBoolOption
+ (
+ "perturb",
+ "Perturb surface points to escape degenerate intersections"
+ );
+
+ argList::addBoolOption
+ (
+ "invertedSpace",
+ "do the surfaces have inverted space orientation, "
+ "i.e. a point at infinity is considered inside. "
+ "This is only sensible for union and intersection."
+ );
+
+ argList::addOption
+ (
+ "trim",
+ "((surface1 volumeType) .. (surfaceN volumeType))",
+ "Trim resulting intersection with additional surfaces;"
+ " volumeType is 'inside' (keep (parts of) edges that are inside)"
+ ", 'outside' (keep (parts of) edges that are outside) or"
+ " 'mixed' (keep all)"
+ );
+
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+
+ word action(args.args()[1]);
+
+ HashTable<booleanSurface::booleanOpType> validActions;
+ validActions.insert("intersection", booleanSurface::INTERSECTION);
+ validActions.insert("union", booleanSurface::UNION);
+ validActions.insert("difference", booleanSurface::DIFFERENCE);
+
+ if (!validActions.found(action))
+ {
+ FatalErrorIn(args.executable())
+ << "Unsupported action " << action << endl
+ << "Supported actions:" << validActions.toc() << abort(FatalError);
+ }
+
+
+ List<Pair<word> > surfaceAndSide;
+ if (args.optionReadIfPresent("trim", surfaceAndSide))
+ {
+ Info<< "Trimming edges with " << surfaceAndSide << endl;
+ }
+
+
+ const word surf1Name(args.args()[2]);
+ Info<< "Reading surface " << surf1Name << endl;
+ triSurfaceMesh surf1
+ (
+ IOobject
+ (
+ surf1Name,
+ runTime.constant(),
+ triSurfaceMesh::meshSubDir,
+ runTime
+ )
+ );
+
+ Info<< surf1Name << " statistics:" << endl;
+ surf1.writeStats(Info);
+ Info<< endl;
+
+ const word surf2Name(args[3]);
+ Info<< "Reading surface " << surf2Name << endl;
+ triSurfaceMesh surf2
(
IOobject
(
- sFeatFileName + ".extendedFeatureEdgeMesh",
+ surf2Name,
runTime.constant(),
- "extendedFeatureEdgeMesh",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- inter.cutPoints(),
- inter.cutEdges(),
-
- 0, // concaveStart,
- 0, // mixedStart,
- 0, // nonFeatureStart,
-
- internalStart, // internalStart,
- nIntOrExt, // flatStart,
- nIntOrExt + nFlat, // openStart,
- nIntOrExt + nFlat + nOpen, // multipleStart,
-
- normalsTmp,
- normalVolumeTypesTmp,
- edgeDirections,
- normalDirectionsTmp,
- edgeNormalsTmp,
-
- labelListList(0), // featurePointNormals,
- labelListList(0), // featurePointEdges,
- labelList(0) // regionEdges
+ triSurfaceMesh::meshSubDir,
+ runTime
+ )
);
+ Info<< surf2Name << " statistics:" << endl;
+ surf2.writeStats(Info);
+ Info<< endl;
+
+ const bool surf1Baffle = args.optionFound("surf1Baffle");
+ const bool surf2Baffle = args.optionFound("surf2Baffle");
+
+ edgeIntersections edgeCuts1;
+ edgeIntersections edgeCuts2;
+
+ bool invertedSpace = args.optionFound("invertedSpace");
+
+ if (invertedSpace && validActions[action] == booleanSurface::DIFFERENCE)
+ {
+ FatalErrorIn(args.executable())
+ << "Inverted space only makes sense for union or intersection."
+ << exit(FatalError);
+ }
+
+
+#ifdef NO_CGAL
+ // Calculate the points where the edges are cut by the other surface
+ calcEdgeCuts
+ (
+ surf1,
+ surf2,
+ args.optionFound("perturb"),
+ edgeCuts1,
+ edgeCuts2
+ );
+#else
+ //calcEdgeCutsCGAL
+ calcEdgeCutsBitsCGAL
+ (
+ surf1,
+ surf2,
+ args.optionFound("perturb"),
+ edgeCuts1,
+ edgeCuts2
+ );
+#endif // NO_CGAL
+
+
+ const fileName sFeatFileName
+ (
+ fileName(surf1Name).lessExt().name()
+ + "_"
+ + fileName(surf2Name).lessExt().name()
+ + "_"
+ + action
+ );
+
+ autoPtr<extendedFeatureEdgeMesh> feMeshPtr
+ (
+ createEdgeMesh
+ (
+ IOobject
+ (
+ sFeatFileName + ".extendedFeatureEdgeMesh",
+ runTime.constant(),
+ "extendedFeatureEdgeMesh",
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ booleanSurface::booleanOpTypeNames[action],
+ surf1Baffle,
+ surf2Baffle,
+ invertedSpace,
+ surf1,
+ edgeCuts1,
+ surf2,
+ edgeCuts2
+ )
+ );
+
+
+ // Trim intersections
+ forAll(surfaceAndSide, trimI)
+ {
+ const word& trimName = surfaceAndSide[trimI].first();
+ const volumeType trimType
+ (
+ volumeType::names[surfaceAndSide[trimI].second()]
+ );
+
+ Info<< "Reading trim surface " << trimName << endl;
+ const triSurfaceMesh trimSurf
+ (
+ IOobject
+ (
+ trimName,
+ runTime.constant(),
+ triSurfaceMesh::meshSubDir,
+ runTime,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ )
+ );
+
+ Info<< trimName << " statistics:" << endl;
+ trimSurf.writeStats(Info);
+ Info<< endl;
+
+ labelList pointMap;
+ labelList edgeMap;
+ feMeshPtr().trim
+ (
+ trimSurf,
+ trimType,
+ pointMap,
+ edgeMap
+ );
+ }
+
+
+ const extendedFeatureEdgeMesh& feMesh = feMeshPtr();
+
+ feMesh.writeStats(Info);
+
feMesh.write();
feMesh.writeObj(feMesh.path()/sFeatFileName);
@@ -782,7 +1738,7 @@ int main(int argc, char *argv[])
(
sFeatFileName + ".eMesh", // name
runTime.constant(), // instance
- "triSurface",
+ triSurfaceMesh::meshSubDir,
runTime, // registry
IOobject::NO_READ,
IOobject::NO_WRITE,
diff --git a/applications/utilities/surface/surfaceCheck/surfaceCheck.C b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
index 6ec067aa84763457164d4e3f90599d5136e76e25..c292ebe6b5ca8e577d96bdb32d8363cf7e0184d2 100644
--- a/applications/utilities/surface/surfaceCheck/surfaceCheck.C
+++ b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -170,6 +170,101 @@ labelList countBins
}
+
+void writeZoning
+(
+ const triSurface& surf,
+ const labelList& faceZone,
+ const word& fieldName,
+ const fileName& surfFilePath,
+ const fileName& surfFileNameBase
+)
+{
+ Info<< "Writing zoning to "
+ << fileName
+ (
+ surfFilePath
+ / fieldName
+ + '_'
+ + surfFileNameBase
+ + '.'
+ + vtkSurfaceWriter::typeName
+ )
+ << "..." << endl << endl;
+
+ // Convert data
+ scalarField scalarFaceZone(faceZone.size());
+ forAll(faceZone, i)
+ {
+ scalarFaceZone[i] = faceZone[i];
+ }
+ faceList faces(surf.size());
+ forAll(surf, i)
+ {
+ faces[i] = surf[i].triFaceFace();
+ }
+
+ vtkSurfaceWriter().write
+ (
+ surfFilePath,
+ surfFileNameBase,
+ surf.points(),
+ faces,
+ fieldName,
+ scalarFaceZone,
+ false // face based data
+ );
+}
+
+
+void writeParts
+(
+ const triSurface& surf,
+ const label nFaceZones,
+ const labelList& faceZone,
+ const fileName& surfFilePath,
+ const fileName& surfFileNameBase
+)
+{
+ for (label zone = 0; zone < nFaceZones; zone++)
+ {
+ boolList includeMap(surf.size(), false);
+
+ forAll(faceZone, faceI)
+ {
+ if (faceZone[faceI] == zone)
+ {
+ includeMap[faceI] = true;
+ }
+ }
+
+ labelList pointMap;
+ labelList faceMap;
+
+ triSurface subSurf
+ (
+ surf.subsetMesh
+ (
+ includeMap,
+ pointMap,
+ faceMap
+ )
+ );
+
+ fileName subName
+ (
+ surfFilePath
+ /surfFileNameBase + "_" + name(zone) + ".obj"
+ );
+
+ Info<< "writing part " << zone << " size " << subSurf.size()
+ << " to " << subName << endl;
+
+ subSurf.write(subName);
+ }
+}
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
@@ -218,6 +313,20 @@ int main(int argc, char *argv[])
surf.writeStats(Info);
Info<< endl;
+
+ // Determine path and extension
+ fileName surfFileNameBase(surfFileName.name());
+ const word fileType = surfFileNameBase.ext();
+ // Strip extension
+ surfFileNameBase = surfFileNameBase.lessExt();
+ // If extension was .gz strip original extension
+ if (fileType == "gz")
+ {
+ surfFileNameBase = surfFileNameBase.lessExt();
+ }
+ const fileName surfFilePath(surfFileName.path());
+
+
// write bounding box corners
if (args.optionFound("blockMesh"))
{
@@ -515,12 +624,12 @@ int main(int argc, char *argv[])
DynamicList<label> problemFaces(surf.size()/100 + 1);
- const labelListList& eFaces = surf.edgeFaces();
+ const labelListList& edgeFaces = surf.edgeFaces();
label nSingleEdges = 0;
- forAll(eFaces, edgeI)
+ forAll(edgeFaces, edgeI)
{
- const labelList& myFaces = eFaces[edgeI];
+ const labelList& myFaces = edgeFaces[edgeI];
if (myFaces.size() == 1)
{
@@ -531,9 +640,9 @@ int main(int argc, char *argv[])
}
label nMultEdges = 0;
- forAll(eFaces, edgeI)
+ forAll(edgeFaces, edgeI)
{
- const labelList& myFaces = eFaces[edgeI];
+ const labelList& myFaces = edgeFaces[edgeI];
if (myFaces.size() > 2)
{
@@ -549,7 +658,8 @@ int main(int argc, char *argv[])
if ((nSingleEdges != 0) || (nMultEdges != 0))
{
- Info<< "Surface is not closed since not all edges connected to "
+ Info<< "Surface is not closed since not all edges ("
+ << edgeFaces.size() << ") connected to "
<< "two faces:" << endl
<< " connected to one face : " << nSingleEdges << endl
<< " connected to >2 faces : " << nMultEdges << endl;
@@ -578,10 +688,9 @@ int main(int argc, char *argv[])
boolList borderEdge(surf.nEdges(), false);
if (splitNonManifold)
{
- const labelListList& eFaces = surf.edgeFaces();
- forAll(eFaces, edgeI)
+ forAll(edgeFaces, edgeI)
{
- if (eFaces[edgeI].size() > 2)
+ if (edgeFaces[edgeI].size() > 2)
{
borderEdge[edgeI] = true;
}
@@ -597,85 +706,15 @@ int main(int argc, char *argv[])
{
Info<< "Splitting surface into parts ..." << endl << endl;
- fileName surfFileNameBase(surfFileName.name());
- const word fileType = surfFileNameBase.ext();
- // Strip extension
- surfFileNameBase = surfFileNameBase.lessExt();
- // If extension was .gz strip original extension
- if (fileType == "gz")
- {
- surfFileNameBase = surfFileNameBase.lessExt();
- }
-
-
- {
- Info<< "Writing zoning to "
- << fileName
- (
- "zone_"
- + surfFileNameBase
- + '.'
- + vtkSurfaceWriter::typeName
- )
- << "..." << endl << endl;
-
- // Convert data
- scalarField scalarFaceZone(faceZone.size());
- forAll(faceZone, i)
- {
- scalarFaceZone[i] = faceZone[i];
- }
- faceList faces(surf.size());
- forAll(surf, i)
- {
- faces[i] = surf[i].triFaceFace();
- }
-
- vtkSurfaceWriter().write
- (
- surfFileName.path(),
- surfFileNameBase,
- surf.points(),
- faces,
- "zone",
- scalarFaceZone,
- false // face based data
- );
- }
-
-
- for (label zone = 0; zone < numZones; zone++)
- {
- boolList includeMap(surf.size(), false);
-
- forAll(faceZone, faceI)
- {
- if (faceZone[faceI] == zone)
- {
- includeMap[faceI] = true;
- }
- }
-
- labelList pointMap;
- labelList faceMap;
-
- triSurface subSurf
- (
- surf.subsetMesh
- (
- includeMap,
- pointMap,
- faceMap
- )
- );
-
- fileName subName(surfFileNameBase + "_" + name(zone) + ".obj");
-
- Info<< "writing part " << zone << " size " << subSurf.size()
- << " to " << subName << endl;
-
- subSurf.write(subName);
- }
+ writeZoning(surf, faceZone, "zone", surfFilePath, surfFileNameBase);
+ writeParts
+ (
+ surf,
+ numZones,
+ faceZone,
+ surfFilePath,
+ surfFileNameBase
+ );
}
}
@@ -700,6 +739,15 @@ int main(int argc, char *argv[])
if (numNormalZones > 1)
{
Info<< "More than one normal orientation." << endl;
+ writeZoning(surf, normalZone, "normal", surfFilePath, surfFileNameBase);
+ writeParts
+ (
+ surf,
+ numNormalZones,
+ normalZone,
+ surfFilePath,
+ surfFileNameBase + "_normal"
+ );
}
Info<< endl;
diff --git a/applications/utilities/surface/surfaceFeatureConvert/surfaceFeatureConvert.C b/applications/utilities/surface/surfaceFeatureConvert/surfaceFeatureConvert.C
index fadf9a7ec7816880a3e4dcb6636f82bb01e4a619..e62bee9a68f5e86486d75d39aa016ee8fb6a7a81 100644
--- a/applications/utilities/surface/surfaceFeatureConvert/surfaceFeatureConvert.C
+++ b/applications/utilities/surface/surfaceFeatureConvert/surfaceFeatureConvert.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -99,7 +99,6 @@ int main(int argc, char *argv[])
mesh.write(exportName);
mesh.writeStats(Info);
- Info<< endl;
Info<< "\nEnd\n" << endl;
diff --git a/applications/utilities/surface/surfaceInertia/surfaceInertia.C b/applications/utilities/surface/surfaceInertia/surfaceInertia.C
index 6b30f610b134dfa1e69a455a7f73fda1f6879c69..14ff3cc5c8b46999f065bc5c79d23c8f1de485b9 100644
--- a/applications/utilities/surface/surfaceInertia/surfaceInertia.C
+++ b/applications/utilities/surface/surfaceInertia/surfaceInertia.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -406,7 +406,7 @@ int main(int argc, char *argv[])
str << "l " << 1 << ' ' << i + 1 << endl;
}
- Info<< nl << "End" << nl << endl;
+ Info<< "\nEnd\n" << endl;
return 0;
}
diff --git a/applications/utilities/surface/surfaceInflate/Make/files b/applications/utilities/surface/surfaceInflate/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..8677a01c390fc67be2ff0b9a905d1d9e80ff6512
--- /dev/null
+++ b/applications/utilities/surface/surfaceInflate/Make/files
@@ -0,0 +1,3 @@
+surfaceInflate.C
+
+EXE = $(FOAM_APPBIN)/surfaceInflate
diff --git a/applications/utilities/surface/surfaceInflate/Make/options b/applications/utilities/surface/surfaceInflate/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..6c6cd122873520ca7edceea5bf8b5be49d77f4ef
--- /dev/null
+++ b/applications/utilities/surface/surfaceInflate/Make/options
@@ -0,0 +1,10 @@
+EXE_INC = \
+ -I$(FOAM_SRC)/triSurface/lnInclude \
+ -I$(FOAM_SRC)/surfMesh/lnInclude \
+ -I$(FOAM_SRC)/meshTools/lnInclude
+
+
+EXE_LIBS = \
+ -ltriSurface \
+ -lsurfMesh \
+ -lmeshTools
diff --git a/applications/utilities/surface/surfaceInflate/surfaceInflate.C b/applications/utilities/surface/surfaceInflate/surfaceInflate.C
new file mode 100644
index 0000000000000000000000000000000000000000..dfc415432178c878c1e389e1b421fec1e82f2783
--- /dev/null
+++ b/applications/utilities/surface/surfaceInflate/surfaceInflate.C
@@ -0,0 +1,879 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+Application
+ surfaceInflate
+
+Description
+ Inflates surface. WIP. Checks for overlaps and locally lowers
+ inflation distance
+
+Usage
+ - surfaceInflate [OPTION]
+
+ \param -checkSelfIntersection \n
+ Includes checks for self-intersection.
+
+ \param -nSmooth
+ Specify number of smoothing iterations
+
+ \param -featureAngle
+ Specify a feature angle
+
+
+ E.g. inflate surface by 2cm with 1.5 safety factor:
+ surfaceInflate DTC-scaled.obj 0.02 1.5 -featureAngle 45 -nSmooth 2
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "triSurfaceFields.H"
+#include "triSurfaceMesh.H"
+#include "triSurfaceGeoMesh.H"
+#include "PackedBoolList.H"
+#include "OBJstream.H"
+#include "surfaceFeatures.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scalar calcVertexNormalWeight
+(
+ const triFace& f,
+ const label pI,
+ const vector& fN,
+ const pointField& points
+)
+{
+ label index = findIndex(f, pI);
+
+ if (index == -1)
+ {
+ FatalErrorIn("calcVertexNormals()")
+ << "Point not in face" << abort(FatalError);
+ }
+
+ const vector e1 = points[f[index]] - points[f[f.fcIndex(index)]];
+ const vector e2 = points[f[index]] - points[f[f.rcIndex(index)]];
+
+ return mag(fN)/(magSqr(e1)*magSqr(e2) + VSMALL);
+}
+
+
+tmp<vectorField> calcVertexNormals(const triSurface& surf)
+{
+ // Weighted average of normals of faces attached to the vertex
+ // Weight = fA / (mag(e0)^2 * mag(e1)^2);
+ tmp<vectorField> tpointNormals
+ (
+ new pointField(surf.nPoints(), vector::zero)
+ );
+ vectorField& pointNormals = tpointNormals();
+
+ const pointField& points = surf.points();
+ const labelListList& pointFaces = surf.pointFaces();
+ const labelList& meshPoints = surf.meshPoints();
+
+ forAll(pointFaces, pI)
+ {
+ const labelList& pFaces = pointFaces[pI];
+
+ forAll(pFaces, fI)
+ {
+ const label faceI = pFaces[fI];
+ const triFace& f = surf[faceI];
+
+ vector fN = f.normal(points);
+
+ scalar weight = calcVertexNormalWeight
+ (
+ f,
+ meshPoints[pI],
+ fN,
+ points
+ );
+
+ pointNormals[pI] += weight*fN;
+ }
+
+ pointNormals[pI] /= mag(pointNormals[pI]) + VSMALL;
+ }
+
+ return tpointNormals;
+}
+
+
+tmp<vectorField> calcPointNormals
+(
+ const triSurface& s,
+ const PackedBoolList& isFeaturePoint,
+ const List<surfaceFeatures::edgeStatus>& edgeStat
+)
+{
+ //const pointField pointNormals(s.pointNormals());
+ tmp<vectorField> tpointNormals(calcVertexNormals(s));
+ vectorField& pointNormals = tpointNormals();
+
+
+ // feature edges: create edge normals from edgeFaces only.
+ {
+ const labelListList& edgeFaces = s.edgeFaces();
+
+ labelList nNormals(s.nPoints(), 0);
+ forAll(edgeStat, edgeI)
+ {
+ if (edgeStat[edgeI] != surfaceFeatures::NONE)
+ {
+ const edge& e = s.edges()[edgeI];
+ forAll(e, i)
+ {
+ if (!isFeaturePoint[e[i]])
+ {
+ pointNormals[e[i]] = vector::zero;
+ }
+ }
+ }
+ }
+
+ forAll(edgeStat, edgeI)
+ {
+ if (edgeStat[edgeI] != surfaceFeatures::NONE)
+ {
+ const labelList& eFaces = edgeFaces[edgeI];
+
+ // Get average edge normal
+ vector n = vector::zero;
+ forAll(eFaces, i)
+ {
+ n += s.faceNormals()[eFaces[i]];
+ }
+ n /= eFaces.size();
+
+
+ // Sum to points
+ const edge& e = s.edges()[edgeI];
+ forAll(e, i)
+ {
+ if (!isFeaturePoint[e[i]])
+ {
+ pointNormals[e[i]] += n;
+ nNormals[e[i]]++;
+ }
+ }
+ }
+ }
+
+ forAll(nNormals, pointI)
+ {
+ if (nNormals[pointI] > 0)
+ {
+ pointNormals[pointI] /= mag(pointNormals[pointI]);
+ }
+ }
+ }
+
+
+ forAll(pointNormals, pointI)
+ {
+ if (mag(mag(pointNormals[pointI])-1) > SMALL)
+ {
+ FatalErrorIn("calcPointNormals()") << "unitialised"
+ << exit(FatalError);
+ }
+ }
+
+ return tpointNormals;
+}
+
+
+void detectSelfIntersections
+(
+ const triSurfaceMesh& s,
+ PackedBoolList& isEdgeIntersecting
+)
+{
+ const edgeList& edges = s.edges();
+ const indexedOctree<treeDataTriSurface>& tree = s.tree();
+ const labelList& meshPoints = s.meshPoints();
+ const pointField& points = s.points();
+
+ isEdgeIntersecting.setSize(edges.size());
+ isEdgeIntersecting = false;
+
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+
+ pointIndexHit hitInfo
+ (
+ tree.findLine
+ (
+ points[meshPoints[e[0]]],
+ points[meshPoints[e[1]]],
+ treeDataTriSurface::findSelfIntersectOp(tree, edgeI)
+ )
+ );
+
+ if (hitInfo.hit())
+ {
+ isEdgeIntersecting[edgeI] = true;
+ }
+ }
+}
+
+
+label detectIntersectionPoints
+(
+ const scalar tolerance,
+ const triSurfaceMesh& s,
+
+ const scalar extendFactor,
+ const pointField& initialPoints,
+ const vectorField& displacement,
+
+ const bool checkSelfIntersect,
+ const PackedBoolList& initialIsEdgeIntersecting,
+
+ PackedBoolList& isPointOnHitEdge,
+ scalarField& scale
+)
+{
+ const pointField initialLocalPoints(initialPoints, s.meshPoints());
+ const List<labelledTri>& localFaces = s.localFaces();
+
+
+ label nHits = 0;
+ isPointOnHitEdge.setSize(s.nPoints());
+ isPointOnHitEdge = false;
+
+
+ // 1. Extrusion offset vectors intersecting new surface location
+ {
+ scalar tol = max(tolerance, 10*s.tolerance());
+
+ // Collect all the edge vectors. Slightly shorten the edges to prevent
+ // finding lots of intersections. The fast triangle intersection routine
+ // has problems with rays passing through a point of the
+ // triangle.
+ pointField start(initialLocalPoints+tol*displacement);
+ pointField end(initialLocalPoints+extendFactor*displacement);
+
+ List<pointIndexHit> hits;
+ s.findLineAny(start, end, hits);
+
+ forAll(hits, pointI)
+ {
+ if
+ (
+ hits[pointI].hit()
+ && findIndex(localFaces[hits[pointI].index()], pointI) == -1
+ )
+ {
+ scale[pointI] = max(0.0, scale[pointI]-0.2);
+
+ isPointOnHitEdge[pointI] = true;
+ nHits++;
+ }
+ }
+ }
+
+
+ // 2. (new) surface self intersections
+ if (checkSelfIntersect)
+ {
+ PackedBoolList isEdgeIntersecting;
+ detectSelfIntersections(s, isEdgeIntersecting);
+
+ const edgeList& edges = s.edges();
+ const pointField& points = s.points();
+
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+
+ if (isEdgeIntersecting[edgeI] && !initialIsEdgeIntersecting[edgeI])
+ {
+ if (isPointOnHitEdge.set(e[0]))
+ {
+ label start = s.meshPoints()[e[0]];
+ const point& pt = points[start];
+
+ Pout<< "Additional self intersection at "
+ << pt
+ << endl;
+
+ scale[e[0]] = max(0.0, scale[e[0]]-0.2);
+ nHits++;
+ }
+ if (isPointOnHitEdge.set(e[1]))
+ {
+ label end = s.meshPoints()[e[1]];
+ const point& pt = points[end];
+
+ Pout<< "Additional self intersection at "
+ << pt
+ << endl;
+
+ scale[e[1]] = max(0.0, scale[e[1]]-0.2);
+ nHits++;
+ }
+ }
+ }
+ }
+
+ return nHits;
+}
+
+
+tmp<scalarField> avg
+(
+ const triSurface& s,
+ const scalarField& fld,
+ const scalarField& edgeWeights
+)
+{
+ tmp<scalarField> tres(new scalarField(s.nPoints(), 0.0));
+ scalarField& res = tres();
+
+ scalarField sumWeight(s.nPoints(), 0.0);
+
+ const edgeList& edges = s.edges();
+
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+ const scalar w = edgeWeights[edgeI];
+
+ res[e[0]] += w*fld[e[1]];
+ sumWeight[e[0]] += w;
+
+ res[e[1]] += w*fld[e[0]];
+ sumWeight[e[1]] += w;
+ }
+
+ // Average
+ // ~~~~~~~
+
+ forAll(res, pointI)
+ {
+ if (mag(sumWeight[pointI]) < VSMALL)
+ {
+ // Unconnected point. Take over original value
+ res[pointI] = fld[pointI];
+ }
+ else
+ {
+ res[pointI] /= sumWeight[pointI];
+ }
+ }
+
+ return tres;
+}
+
+
+void minSmooth
+(
+ const triSurface& s,
+ const PackedBoolList& isAffectedPoint,
+ const scalarField& fld,
+ scalarField& newFld
+)
+{
+
+ const edgeList& edges = s.edges();
+ const pointField& points = s.points();
+ const labelList& mp = s.meshPoints();
+
+ scalarField edgeWeights(edges.size());
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+ scalar w = mag(points[mp[e[0]]]-points[mp[e[1]]]);
+
+ edgeWeights[edgeI] = 1.0/(max(w, SMALL));
+ }
+
+ tmp<scalarField> tavgFld = avg(s, fld, edgeWeights);
+
+ const scalarField& avgFld = tavgFld();
+
+ forAll(fld, pointI)
+ {
+ if (isAffectedPoint.get(pointI) == 1)
+ {
+ newFld[pointI] = min
+ (
+ fld[pointI],
+ 0.5*fld[pointI] + 0.5*avgFld[pointI]
+ //avgFld[pointI]
+ );
+ }
+ }
+}
+
+
+void lloydsSmoothing
+(
+ const label nSmooth,
+ triSurface& s,
+ const PackedBoolList& isFeaturePoint,
+ const List<surfaceFeatures::edgeStatus>& edgeStat,
+ const PackedBoolList& isAffectedPoint
+)
+{
+ const labelList& meshPoints = s.meshPoints();
+ const edgeList& edges = s.edges();
+
+
+ PackedBoolList isSmoothPoint(isAffectedPoint);
+ // Extend isSmoothPoint
+ {
+ PackedBoolList newIsSmoothPoint(isSmoothPoint);
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+ if (isSmoothPoint[e[0]])
+ {
+ newIsSmoothPoint[e[1]] = true;
+ }
+ if (isSmoothPoint[e[1]])
+ {
+ newIsSmoothPoint[e[0]] = true;
+ }
+ }
+ Info<< "From points-to-smooth " << isSmoothPoint.count()
+ << " to " << newIsSmoothPoint.count()
+ << endl;
+ isSmoothPoint.transfer(newIsSmoothPoint);
+ }
+
+ // Do some smoothing (Lloyds algorithm) around problematic points
+ for (label i = 0; i < nSmooth; i++)
+ {
+ const labelListList& pointFaces = s.pointFaces();
+ const pointField& faceCentres = s.faceCentres();
+
+ pointField newPoints(s.points());
+ forAll(isSmoothPoint, pointI)
+ {
+ if (isSmoothPoint[pointI] && !isFeaturePoint[pointI])
+ {
+ const labelList& pFaces = pointFaces[pointI];
+
+ point avg = point::zero;
+ forAll(pFaces, pFaceI)
+ {
+ avg += faceCentres[pFaces[pFaceI]];
+ }
+ avg /= pFaces.size();
+ newPoints[meshPoints[pointI]] = avg;
+ }
+ }
+
+
+ // Move points on feature edges only according to feature edges.
+
+ const pointField& points = s.points();
+
+ vectorField pointSum(s.nPoints(), vector::zero);
+ labelList nPointSum(s.nPoints(), 0);
+
+ forAll(edges, edgeI)
+ {
+ if (edgeStat[edgeI] != surfaceFeatures::NONE)
+ {
+ const edge& e = edges[edgeI];
+ const point& start = points[meshPoints[e[0]]];
+ const point& end = points[meshPoints[e[1]]];
+ point eMid = 0.5*(start+end);
+ pointSum[e[0]] += eMid;
+ nPointSum[e[0]]++;
+ pointSum[e[1]] += eMid;
+ nPointSum[e[1]]++;
+ }
+ }
+
+ forAll(pointSum, pointI)
+ {
+ if
+ (
+ isSmoothPoint[pointI]
+ && isFeaturePoint[pointI]
+ && nPointSum[pointI] >= 2
+ )
+ {
+ newPoints[meshPoints[pointI]] =
+ pointSum[pointI]/nPointSum[pointI];
+ }
+ }
+
+
+ s.movePoints(newPoints);
+
+
+ // Extend isSmoothPoint
+ {
+ PackedBoolList newIsSmoothPoint(isSmoothPoint);
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+ if (isSmoothPoint[e[0]])
+ {
+ newIsSmoothPoint[e[1]] = true;
+ }
+ if (isSmoothPoint[e[1]])
+ {
+ newIsSmoothPoint[e[0]] = true;
+ }
+ }
+ Info<< "From points-to-smooth " << isSmoothPoint.count()
+ << " to " << newIsSmoothPoint.count()
+ << endl;
+ isSmoothPoint.transfer(newIsSmoothPoint);
+ }
+ }
+}
+
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+ argList::addNote("Inflates surface according to point normals.");
+
+ argList::noParallel();
+ argList::addNote
+ (
+ "Creates inflated version of surface using point normals."
+ " Takes surface, distance to inflate and additional safety factor"
+ );
+ argList::addBoolOption
+ (
+ "checkSelfIntersection",
+ "also check for self-intersection"
+ );
+ argList::addOption
+ (
+ "nSmooth",
+ "integer",
+ "number of smoothing iterations (default 20)"
+ );
+ argList::addOption
+ (
+ "featureAngle",
+ "scalar",
+ "feature angle"
+ );
+ argList::addBoolOption
+ (
+ "debug",
+ "switch on additional debug information"
+ );
+
+ argList::validArgs.append("inputFile");
+ argList::validArgs.append("distance");
+ argList::validArgs.append("safety factor [1..]");
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+ runTime.functionObjects().off();
+
+ const word inputName(args.args()[1]);
+ const scalar distance(args.argRead<scalar>(2));
+ const scalar extendFactor(args.argRead<scalar>(3));
+ const bool checkSelfIntersect = args.optionFound("checkSelfIntersection");
+ const label nSmooth = args.optionLookupOrDefault("nSmooth", 10);
+ const scalar featureAngle = args.optionLookupOrDefault<scalar>
+ (
+ "featureAngle",
+ 180
+ );
+ const bool debug = args.optionFound("debug");
+
+
+ Info<< "Inflating surface " << inputName
+ << " according to point normals" << nl
+ << " distance : " << distance << nl
+ << " safety factor : " << extendFactor << nl
+ << " self intersection test : " << checkSelfIntersect << nl
+ << " smoothing iterations : " << nSmooth << nl
+ << " feature angle : " << featureAngle << nl
+ << endl;
+
+
+ if (extendFactor < 1 || extendFactor > 10)
+ {
+ FatalErrorIn(args.executable())
+ << "Illegal safety factor " << extendFactor
+ << ". It is usually 1..2"
+ << exit(FatalError);
+ }
+
+
+
+ // Load triSurfaceMesh
+ triSurfaceMesh s
+ (
+ IOobject
+ (
+ inputName, // name
+ runTime.constant(), // instance
+ "triSurface", // local
+ runTime, // registry
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ )
+ );
+
+ // Mark features
+ const surfaceFeatures features(s, 180.0-featureAngle);
+
+ Info<< "Detected features:" << nl
+ << " feature points : " << features.featurePoints().size()
+ << " out of " << s.nPoints() << nl
+ << " feature edges : " << features.featureEdges().size()
+ << " out of " << s.nEdges() << nl
+ << endl;
+
+ PackedBoolList isFeaturePoint(s.nPoints());
+ forAll(features.featurePoints(), i)
+ {
+ label pointI = features.featurePoints()[i];
+ isFeaturePoint[pointI] = true;
+ }
+ const List<surfaceFeatures::edgeStatus> edgeStat(features.toStatus());
+
+
+
+
+ const pointField initialPoints(s.points());
+
+
+ // Construct scale
+ Info<< "Constructing field scale\n" << endl;
+ triSurfacePointScalarField scale
+ (
+ IOobject
+ (
+ "scale", // name
+ runTime.timeName(), // instance
+ s, // registry
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ s,
+ dimensionedScalar("scale", dimLength, 1.0)
+ );
+
+
+ // Construct unit normals
+
+ Info<< "Calculating vertex normals\n" << endl;
+ const pointField pointNormals
+ (
+ calcPointNormals
+ (
+ s,
+ isFeaturePoint,
+ edgeStat
+ )
+ );
+
+
+ // Construct pointDisplacement
+ Info<< "Constructing field pointDisplacement\n" << endl;
+ triSurfacePointVectorField pointDisplacement
+ (
+ IOobject
+ (
+ "pointDisplacement", // name
+ runTime.timeName(), // instance
+ s, // registry
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ s,
+ dimLength,
+ distance*scale*pointNormals
+ );
+
+
+ const labelList& meshPoints = s.meshPoints();
+
+
+ // Any point on any intersected edge in any of the iterations
+ PackedBoolList isScaledPoint(s.nPoints());
+
+
+ // Detect any self intersections on initial mesh
+ PackedBoolList initialIsEdgeIntersecting;
+ if (checkSelfIntersect)
+ {
+ detectSelfIntersections(s, initialIsEdgeIntersecting);
+ }
+ else
+ {
+ // Mark all edges as already self intersecting so avoid detecting any
+ // new ones
+ initialIsEdgeIntersecting.setSize(s.nEdges(), true);
+ }
+
+
+ // Inflate
+ while (runTime.loop())
+ {
+ Info<< "Time = " << runTime.timeName() << nl << endl;
+
+ // Move to new location
+ pointField newPoints(initialPoints);
+ forAll(meshPoints, i)
+ {
+ newPoints[meshPoints[i]] += pointDisplacement[i];
+ }
+ s.movePoints(newPoints);
+ Info<< "Bounding box : " << s.bounds() << endl;
+
+
+ // Work out scale from pointDisplacement
+ forAll(scale, pointI)
+ {
+ if (s.pointFaces()[pointI].size())
+ {
+ scale[pointI] = mag(pointDisplacement[pointI])/distance;
+ }
+ else
+ {
+ scale[pointI] = 1.0;
+ }
+ }
+
+
+ Info<< "Scale : " << gAverage(scale) << endl;
+ Info<< "Displacement : " << gAverage(pointDisplacement) << endl;
+
+
+
+ // Detect any intersections and scale back
+ PackedBoolList isAffectedPoint;
+ label nIntersections = detectIntersectionPoints
+ (
+ 1e-9, // intersection tolerance
+ s,
+ extendFactor,
+ initialPoints,
+ pointDisplacement,
+
+ checkSelfIntersect,
+ initialIsEdgeIntersecting,
+
+ isAffectedPoint,
+ scale
+ );
+ Info<< "Detected " << nIntersections << " intersections" << nl
+ << endl;
+
+ if (nIntersections == 0)
+ {
+ runTime.write();
+ break;
+ }
+
+
+ // Accumulate all affected points
+ forAll(isAffectedPoint, pointI)
+ {
+ if (isAffectedPoint[pointI])
+ {
+ isScaledPoint[pointI] = true;
+ }
+ }
+
+ // Smear out lowering of scale so any edges not found are
+ // still included
+ for (label i = 0; i < nSmooth; i++)
+ {
+ triSurfacePointScalarField oldScale(scale);
+ oldScale.rename("oldScale");
+ minSmooth
+ (
+ s,
+ PackedBoolList(s.nPoints(), true),
+ oldScale,
+ scale
+ );
+ }
+
+
+ // From scale update the pointDisplacement
+ pointDisplacement *= distance*scale/mag(pointDisplacement);
+
+
+ // Do some smoothing (Lloyds algorithm)
+ lloydsSmoothing(nSmooth, s, isFeaturePoint, edgeStat, isAffectedPoint);
+
+
+ // Update pointDisplacement
+ const pointField& pts = s.points();
+ forAll(meshPoints, i)
+ {
+ label meshPointI = meshPoints[i];
+ pointDisplacement[i] = pts[meshPointI]-initialPoints[meshPointI];
+ }
+
+
+ // Write
+ runTime.write();
+
+ Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+ << " ClockTime = " << runTime.elapsedClockTime() << " s"
+ << nl << endl;
+ }
+
+
+ Info<< "Detected overall intersecting points " << isScaledPoint.count()
+ << " out of " << s.nPoints() << nl << endl;
+
+
+ if (debug)
+ {
+ OBJstream str(runTime.path()/"isScaledPoint.obj");
+ forAll(isScaledPoint, pointI)
+ {
+ if (isScaledPoint[pointI])
+ {
+ str.write(initialPoints[meshPoints[pointI]]);
+ }
+ }
+ }
+
+
+ Info << "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C b/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
index d1e725faee25dd0a8ced157eb14fcbfab91ae64b..93f29a2283c94fbf9e7bd3f529d538449ae4a900 100644
--- a/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
+++ b/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -189,8 +189,6 @@ int main(int argc, char *argv[])
<< endl;
}
- pointField newPoints(surf1.localPoints());
-
for (label iter = 0; iter < iters; iter++)
{
// Lambda
diff --git a/applications/utilities/surface/surfacePatch/Make/files b/applications/utilities/surface/surfacePatch/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..d983714ab3da30dcdf0353c1868ea2b255623ce2
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/Make/files
@@ -0,0 +1,7 @@
+searchableSurfaceModifier/searchableSurfaceModifier.C
+searchableSurfaceModifier/autoPatch.C
+searchableSurfaceModifier/cut.C
+
+surfacePatch.C
+
+EXE = $(FOAM_APPBIN)/surfacePatch
diff --git a/applications/utilities/surface/surfacePatch/Make/options b/applications/utilities/surface/surfacePatch/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..a8173dc04565dcd2835eafe4e6cc2610ff729c71
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/Make/options
@@ -0,0 +1,10 @@
+EXE_INC = \
+ -IsearchableSurfaceModifier \
+ -I$(LIB_SRC)/triSurface/lnInclude \
+ -I$(LIB_SRC)/fileFormats/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude
+
+
+EXE_LIBS = \
+ -ltriSurface \
+ -lmeshTools
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.C b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.C
new file mode 100644
index 0000000000000000000000000000000000000000..e8ca9a0125c5906bb12e58053eedab8fd8f2cc2e
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.C
@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "autoPatch.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceFeatures.H"
+#include "triSurfaceMesh.H"
+#include "PatchTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace searchableSurfaceModifiers
+{
+ defineTypeNameAndDebug(autoPatch, 0);
+ addToRunTimeSelectionTable
+ (
+ searchableSurfaceModifier,
+ autoPatch,
+ dictionary
+ );
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifiers::autoPatch::autoPatch
+(
+ const searchableSurfaces& geometry,
+ const dictionary& dict
+)
+:
+ searchableSurfaceModifier(geometry, dict),
+ featureAngle_(readScalar(dict.lookup("featureAngle")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifiers::autoPatch::~autoPatch()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+bool Foam::searchableSurfaceModifiers::autoPatch::modify
+(
+ const labelList& regions,
+ searchableSurface& geom
+) const
+{
+ triSurface& surf = refCast<triSurfaceMesh>(geom);
+
+ surfaceFeatures set(surf, 180.0-featureAngle_);
+
+ Info<< nl
+ << "Feature set:" << nl
+ << " feature points : " << set.featurePoints().size() << nl
+ << " feature edges : " << set.featureEdges().size() << nl
+ << " of which" << nl
+ << " region edges : " << set.nRegionEdges() << nl
+ << " external edges : " << set.nExternalEdges() << nl
+ << " internal edges : " << set.nInternalEdges() << nl
+ << endl;
+
+
+ // Get per-edge status.
+ boolList borderEdge(surf.nEdges(), false);
+ forAll(set.featureEdges(), i)
+ {
+ borderEdge[set.featureEdges()[i]] = true;
+ }
+
+
+ bool changed = false;
+
+ forAll(regions, i)
+ {
+ label regionI = regions[i];
+
+ labelList subRegion(surf.size(), -1);
+ label nSub = 0;
+
+ forAll(surf, faceI)
+ {
+ if (surf[faceI].region() == regionI && subRegion[faceI] == -1)
+ {
+ PatchTools::markZone(surf, borderEdge, faceI, nSub, subRegion);
+ nSub++;
+ }
+ }
+
+ Info<< "For region " << surf.patches()[regionI].name()
+ << " detected sub regions " << nSub << endl;
+
+
+ if (nSub > 1)
+ {
+ label nOldPatches = surf.patches().size();
+
+ for (label subI = 0; subI < nSub; subI++)
+ {
+ geometricSurfacePatch patch
+ (
+ surf.patches()[regionI].geometricType(),
+ surf.patches()[regionI].name() + Foam::name(subI),
+ surf.patches().size()
+ );
+
+
+ Info<< "For region " << surf.patches()[regionI].name()
+ << " creating region " << patch.name() << endl;
+
+ surf.patches().append(patch);
+ }
+
+ forAll(subRegion, faceI)
+ {
+ if (subRegion[faceI] != -1)
+ {
+ surf[faceI].region() = nOldPatches + subRegion[faceI];
+ }
+ }
+ changed = true;
+ }
+ }
+
+ return changed;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.H b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.H
new file mode 100644
index 0000000000000000000000000000000000000000..4a7e0459b3447cca87c50c3ecc45b3291fd9a694
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/autoPatch.H
@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::searchableSurfaceModifier::autoPatch
+
+Description
+ Patchify triangles based on a feature angle. Replacement of surfaceAutoPatch
+
+SourceFiles
+ autoPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef autoPatch_H
+#define autoPatch_H
+
+#include "searchableSurfaceModifier.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace searchableSurfaceModifiers
+{
+
+/*---------------------------------------------------------------------------*\
+ Class autoPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class autoPatch
+:
+ public searchableSurfaceModifier
+{
+ // Private data
+
+ //- Feature angle
+ const scalar featureAngle_;
+
+
+ // Private Member Functions
+
+public:
+
+ //- Runtime type information
+ TypeName("autoPatch");
+
+
+ // Constructors
+
+ //- Construct from dictionary
+ autoPatch(const searchableSurfaces&, const dictionary&);
+
+ //- Clone
+ autoPtr<searchableSurfaceModifier> clone() const
+ {
+ notImplemented("autoPtr<searchableSurfaceModifier> clone() const");
+ return autoPtr<searchableSurfaceModifier>(NULL);
+ }
+
+
+ //- Destructor
+ virtual ~autoPatch();
+
+
+ // Member Functions
+
+ //- Apply this selector
+ virtual bool modify(const labelList& regions, searchableSurface&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace searchableSurfaceModifiers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.C b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.C
new file mode 100644
index 0000000000000000000000000000000000000000..c04c47da0b5c173e6da11ffb4c2b8f9f5691d012
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.C
@@ -0,0 +1,390 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cut.H"
+#include "addToRunTimeSelectionTable.H"
+#include "searchableSurfaces.H"
+#include "triSurfaceMesh.H"
+#include "searchableBox.H"
+#include "searchableRotatedBox.H"
+#include "surfaceIntersection.H"
+#include "intersectedSurface.H"
+#include "edgeIntersections.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace searchableSurfaceModifiers
+{
+ defineTypeNameAndDebug(cut, 0);
+ addToRunTimeSelectionTable(searchableSurfaceModifier, cut, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::searchableSurfaceModifiers::cut::triangulate
+(
+ const faceList& fcs,
+ pointField& pts,
+ triSurface& cutSurf
+) const
+{
+ label nTris = 0;
+ forAll(fcs, i)
+ {
+ nTris += fcs[i].size()-2;
+ }
+
+ DynamicList<labelledTri> tris(nTris);
+
+ forAll(fcs, i)
+ {
+ const face& f = fcs[i];
+ // Triangulate around vertex 0
+ for (label fp = 1; fp < f.size()-1; fp++)
+ {
+ tris.append(labelledTri(f[0], f[fp], f[f.fcIndex(fp)], i));
+ }
+ }
+ geometricSurfacePatchList patches(fcs.size());
+ forAll(patches, patchI)
+ {
+ patches[patchI] = geometricSurfacePatch
+ (
+ "",
+ "patch" + Foam::name(patchI),
+ patchI
+ );
+ }
+ cutSurf = triSurface(tris.xfer(), patches, pts.xfer());
+}
+
+
+Foam::triSurface& Foam::searchableSurfaceModifiers::cut::triangulate
+(
+ const searchableSurface& cutter,
+ triSurface& cutSurf
+) const
+{
+ if (isA<searchableBox>(cutter))
+ {
+ const searchableBox& bb = refCast<const searchableBox>(cutter);
+
+ pointField pts(bb.points());
+ triangulate(treeBoundBox::faces, pts, cutSurf);
+
+ return cutSurf;
+ }
+ else if (isA<searchableRotatedBox>(cutter))
+ {
+ const searchableRotatedBox& bb =
+ refCast<const searchableRotatedBox>(cutter);
+
+ pointField pts(bb.points());
+ triangulate(treeBoundBox::faces, pts, cutSurf);
+
+ return cutSurf;
+ }
+ else if (isA<triSurfaceMesh>(cutter))
+ {
+ return const_cast<triSurfaceMesh&>
+ (
+ refCast<const triSurfaceMesh>(cutter)
+ );
+ }
+ else
+ {
+ FatalErrorIn
+ (
+ "searchableSurfaceModifiers::cut::triangulate\n"
+ "(\n"
+ " const searchableSurface&,\n"
+ " triSurface&\n"
+ ")"
+ ) << "Triangulation only supported for triSurfaceMesh, searchableBox"
+ << ", not for surface " << cutter.name()
+ << " of type " << cutter.type()
+ << exit(FatalError);
+ return const_cast<triSurfaceMesh&>
+ (
+ refCast<const triSurfaceMesh>(cutter)
+ );
+ }
+}
+
+
+// Keep on shuffling surface points until no more degenerate intersections.
+// Moves both surfaces and updates set of edge cuts.
+bool Foam::searchableSurfaceModifiers::cut::intersectSurfaces
+(
+ triSurface& surf1,
+ edgeIntersections& edgeCuts1,
+ triSurface& surf2,
+ edgeIntersections& edgeCuts2
+) const
+{
+ bool hasMoved1 = false;
+ bool hasMoved2 = false;
+
+ for (label iter = 0; iter < 10; iter++)
+ {
+ Info<< "Determining intersections of surf1 edges with surf2"
+ << " faces" << endl;
+
+ // Determine surface1 edge intersections. Allow surface to be moved.
+
+ // Number of iterations needed to resolve degenerates
+ label nIters1 = 0;
+ {
+ triSurfaceSearch querySurf2(surf2);
+
+ scalarField surf1PointTol
+ (
+ 1E-6*edgeIntersections::minEdgeLength(surf1)
+ );
+
+ // Determine raw intersections
+ edgeCuts1 = edgeIntersections
+ (
+ surf1,
+ querySurf2,
+ surf1PointTol
+ );
+
+ // Shuffle a bit to resolve degenerate edge-face hits
+ {
+ pointField points1(surf1.points());
+
+ nIters1 =
+ edgeCuts1.removeDegenerates
+ (
+ 5, // max iterations
+ surf1,
+ querySurf2,
+ surf1PointTol,
+ points1 // work array
+ );
+
+ if (nIters1 != 0)
+ {
+ // Update geometric quantities
+ surf1.movePoints(points1);
+ hasMoved1 = true;
+ }
+ }
+ }
+
+ Info<< "Determining intersections of surf2 edges with surf1"
+ << " faces" << endl;
+
+ label nIters2 = 0;
+ {
+ triSurfaceSearch querySurf1(surf1);
+
+ scalarField surf2PointTol
+ (
+ 1E-6*edgeIntersections::minEdgeLength(surf2)
+ );
+
+ // Determine raw intersections
+ edgeCuts2 = edgeIntersections
+ (
+ surf2,
+ querySurf1,
+ surf2PointTol
+ );
+
+ // Shuffle a bit to resolve degenerate edge-face hits
+ {
+ pointField points2(surf2.points());
+
+ nIters2 =
+ edgeCuts2.removeDegenerates
+ (
+ 5, // max iterations
+ surf2,
+ querySurf1,
+ surf2PointTol,
+ points2 // work array
+ );
+
+ if (nIters2 != 0)
+ {
+ // Update geometric quantities
+ surf2.movePoints(points2);
+ hasMoved2 = true;
+ }
+ }
+ }
+
+
+ if (nIters1 == 0 && nIters2 == 0)
+ {
+ //Info<< "** Resolved all intersections to be proper"
+ // << "edge-face pierce" << endl;
+ break;
+ }
+ }
+
+ //if (hasMoved1)
+ //{
+ // fileName newFile("surf1.ftr");
+ // Info<< "Surface 1 has been moved. Writing to " << newFile
+ // << endl;
+ // surf1.write(newFile);
+ //}
+ //
+ //if (hasMoved2)
+ //{
+ // fileName newFile("surf2.ftr");
+ // Info<< "Surface 2 has been moved. Writing to " << newFile
+ // << endl;
+ // surf2.write(newFile);
+ //}
+
+ return hasMoved1 || hasMoved2;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifiers::cut::cut
+(
+ const searchableSurfaces& geometry,
+ const dictionary& dict
+)
+:
+ searchableSurfaceModifier(geometry, dict),
+ cutterNames_(dict_.lookup("cutters"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifiers::cut::~cut()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+bool Foam::searchableSurfaceModifiers::cut::modify
+(
+ const labelList& regions,
+ searchableSurface& geom
+) const
+{
+ triSurface& surf = refCast<triSurfaceMesh>(geom);
+
+ bool changed = false;
+
+ // Find the surfaces to cut with
+ forAll(cutterNames_, cutNameI)
+ {
+ labelList geomIDs =
+ findStrings(cutterNames_[cutNameI], geometry_.names());
+
+ forAll(geomIDs, j)
+ {
+ label geomI = geomIDs[j];
+ const searchableSurface& cutter = geometry_[geomI];
+
+ // Triangulate
+ triSurface work;
+ triSurface& cutSurf = triangulate(cutter, work);
+
+ // Determine intersection (with perturbation)
+ edgeIntersections edge1Cuts;
+ edgeIntersections edge2Cuts;
+ intersectSurfaces
+ (
+ surf,
+ edge1Cuts,
+ cutSurf,
+ edge2Cuts
+ );
+
+
+ // Determine intersection edges
+ surfaceIntersection inter(surf, edge1Cuts, cutSurf, edge2Cuts);
+
+
+ // Use intersection edges to cut up faces. (does all the hard work)
+ intersectedSurface surf3(surf, true, inter);
+
+
+ // Mark triangles based on whether they are inside or outside
+ List<volumeType> volTypes;
+ cutter.getVolumeType(surf3.faceCentres(), volTypes);
+
+ label nInside = 0;
+ forAll(volTypes, i)
+ {
+ if (volTypes[i] == volumeType::INSIDE)
+ {
+ nInside++;
+ }
+ }
+
+ // Add a patch for inside the box
+ if (nInside > 0 && surf3.patches().size() > 0)
+ {
+ geometricSurfacePatchList newPatches(surf3.patches());
+ label sz = newPatches.size();
+ newPatches.setSize(sz+1);
+ newPatches[sz] = geometricSurfacePatch
+ (
+ newPatches[sz-1].geometricType(),
+ newPatches[sz-1].name() + "_inside",
+ newPatches[sz-1].index()
+ );
+
+ Info<< "Moving " << nInside << " out of " << surf3.size()
+ << " triangles to region "
+ << newPatches[sz].name() << endl;
+
+
+ List<labelledTri> newTris(surf3);
+ forAll(volTypes, i)
+ {
+ if (volTypes[i] == volumeType::INSIDE)
+ {
+ newTris[i].region() = sz;
+ }
+ }
+ pointField newPoints(surf3.points());
+ surf = triSurface(newTris.xfer(), newPatches, newPoints.xfer());
+
+ changed = true;
+ }
+ }
+ }
+
+ return changed;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.H b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.H
new file mode 100644
index 0000000000000000000000000000000000000000..4081b753b3598837295c1eca243c1807e46d892d
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/cut.H
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::searchableSurfaceModifier::cut
+
+Description
+ Patchify triangles based on orientation w.r.t other (triangulated or
+ triangulatable) surfaces
+
+SourceFiles
+ cut.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cut_H
+#define cut_H
+
+#include "searchableSurfaceModifier.H"
+#include "wordReList.H"
+#include "faceList.H"
+#include "pointField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class edgeIntersections;
+
+namespace searchableSurfaceModifiers
+{
+
+/*---------------------------------------------------------------------------*\
+ Class cut Declaration
+\*---------------------------------------------------------------------------*/
+
+class cut
+:
+ public searchableSurfaceModifier
+{
+ // Private data
+
+ //- Name of surfaces to cut with
+ const wordReList cutterNames_;
+
+
+ // Private Member Functions
+
+ //- Triangulate faces around 0th vertex
+ void triangulate(const faceList&, pointField&, triSurface&) const;
+
+ //- Triangulate searchableSurface (currently only supported for
+ // searchableBox and triSurfaceMesh)
+ triSurface& triangulate(const searchableSurface&, triSurface&) const;
+
+ //- Intersect surfaces. Perturb to avoid degenerates.
+ bool intersectSurfaces
+ (
+ triSurface& surf1,
+ edgeIntersections& edgeCuts1,
+ triSurface& surf2,
+ edgeIntersections& edgeCuts2
+ ) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("cut");
+
+
+ // Constructors
+
+ //- Construct from dictionary
+ cut(const searchableSurfaces&, const dictionary&);
+
+ //- Clone
+ autoPtr<searchableSurfaceModifier> clone() const
+ {
+ notImplemented("autoPtr<searchableSurfaceModifier> clone() const");
+ return autoPtr<searchableSurfaceModifier>(NULL);
+ }
+
+
+ //- Destructor
+ virtual ~cut();
+
+
+ // Member Functions
+
+ //- Apply this selector
+ virtual bool modify(const labelList& regions, searchableSurface&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace searchableSurfaceModifiers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.C b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.C
new file mode 100644
index 0000000000000000000000000000000000000000..e28036012c237a233db761d9ac440074a1eb5ecf
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.C
@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "searchableSurfaceModifier.H"
+#include "triSurface.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(searchableSurfaceModifier, 0);
+ defineRunTimeSelectionTable(searchableSurfaceModifier, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifier::searchableSurfaceModifier
+(
+ const searchableSurfaces& geometry,
+ const dictionary& dict
+)
+:
+ geometry_(geometry),
+ dict_(dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceModifier::~searchableSurfaceModifier()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::searchableSurfaceModifier>
+Foam::searchableSurfaceModifier::New
+(
+ const word& type,
+ const searchableSurfaces& geometry,
+ const dictionary& dict
+)
+{
+ dictionaryConstructorTable::iterator cstrIter =
+ dictionaryConstructorTablePtr_->find(type);
+
+ if (cstrIter == dictionaryConstructorTablePtr_->end())
+ {
+ FatalErrorIn
+ (
+ "searchableSurfaceModifier::New"
+ "(const word&, const searchableSurfaces&, const dictionary&)"
+ ) << "Unknown searchableSurfaceModifier type "
+ << type << nl << nl
+ << "Valid searchableSurfaceModifier types : " << endl
+ << dictionaryConstructorTablePtr_->sortedToc()
+ << exit(FatalError);
+ }
+
+ return autoPtr<searchableSurfaceModifier>(cstrIter()(geometry, dict));
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.H b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.H
new file mode 100644
index 0000000000000000000000000000000000000000..8955dea7416cdf715b8b80c46f8172d3b52978dd
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/searchableSurfaceModifier/searchableSurfaceModifier.H
@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::searchableSurfaceModifier
+
+Description
+ Changing a surface
+
+SourceFiles
+ searchableSurfaceModifier.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableSurfaceModifier_H
+#define searchableSurfaceModifier_H
+
+#include "dictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class triSurface;
+class searchableSurface;
+class searchableSurfaces;
+
+/*---------------------------------------------------------------------------*\
+ Class searchableSurfaceModifier Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableSurfaceModifier
+{
+protected:
+
+ // Protected data
+
+ const searchableSurfaces& geometry_;
+
+ //- Input dictionary
+ const dictionary dict_;
+
+public:
+
+ //- Runtime type information
+ TypeName("searchableSurfaceModifier");
+
+
+ // Declare run-time constructor selection table
+
+ declareRunTimeSelectionTable
+ (
+ autoPtr,
+ searchableSurfaceModifier,
+ dictionary,
+ (
+ const searchableSurfaces& geometry,
+ const dictionary& dict
+ ),
+ (geometry, dict)
+ );
+
+
+ // Constructors
+
+ //- Construct from dictionary
+ searchableSurfaceModifier(const searchableSurfaces&, const dictionary&);
+
+ //- Clone
+ autoPtr<searchableSurfaceModifier> clone() const
+ {
+ notImplemented("autoPtr<searchableSurfaceModifier> clone() const");
+ return autoPtr<searchableSurfaceModifier>(NULL);
+ }
+
+
+ // Selectors
+
+ //- Return a reference to the selected searchableSurfaceModifier
+ static autoPtr<searchableSurfaceModifier> New
+ (
+ const word& type,
+ const searchableSurfaces&,
+ const dictionary& dict
+ );
+
+
+ //- Destructor
+ virtual ~searchableSurfaceModifier();
+
+
+ // Member Functions
+
+ //- Do any operation on surface. Return true if anything changed
+ virtual bool modify(const labelList&, searchableSurface&) const = 0;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/surfacePatch.C b/applications/utilities/surface/surfacePatch/surfacePatch.C
new file mode 100644
index 0000000000000000000000000000000000000000..3df36f6ec652ed3cb79d8a7ea3c3640d33fb50f4
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/surfacePatch.C
@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 <http://www.gnu.org/licenses/>.
+
+Application
+ surfacePatch
+
+Description
+ Patches (regionises) a surface using a user-selectable method.
+
+ The current methods are either based on a geometric feature angle
+ (a replacement for the surfaceAutoPatch functionality) or on intersecting
+ a set of searchableSurfaces - this will re-triangulate the intersecting
+ triangles and regonise the resulting triangles according to being
+ inside/outside the surface.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triSurface.H"
+#include "argList.H"
+#include "Time.H"
+#include "triSurfaceMesh.H"
+#include "searchableSurfaces.H"
+#include "searchableSurfaceModifier.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+ argList::validArgs.clear();
+ argList::noParallel();
+ #include "addDictOption.H"
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+
+ const word dictName("surfacePatchDict");
+ #include "setSystemRunTimeDictionaryIO.H"
+ const IOdictionary meshDict(dictIO);
+
+
+ // Read geometry
+ // ~~~~~~~~~~~~~
+
+ searchableSurfaces allGeometry
+ (
+ IOobject
+ (
+ "abc", // dummy name
+ runTime.constant(), // instance
+ triSurfaceMesh::meshSubDir, // local
+ runTime, // registry
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ ),
+ meshDict.subDict("geometry"),
+ meshDict.lookupOrDefault("singleRegionName", true)
+ );
+
+
+ const dictionary& surfacesDict = meshDict.subDict("surfaces");
+
+ forAllConstIter(dictionary, surfacesDict, surfacesIter)
+ {
+ if (surfacesIter().isDict())
+ {
+ const word& surfName = surfacesIter().keyword();
+ const dictionary& surfDict = surfacesIter().dict();
+
+ // Look up surface
+ searchableSurface& surf = allGeometry[surfName];
+
+
+ bool changed = false;
+
+ // Check for modifier on global level
+ if (surfDict.found("type"))
+ {
+ autoPtr<searchableSurfaceModifier> modifier
+ (
+ searchableSurfaceModifier::New
+ (
+ surfDict.lookup("type"),
+ allGeometry,
+ surfDict
+ )
+ );
+
+ if (modifier().modify(identity(surf.regions().size()), surf))
+ {
+ changed = true;
+ }
+ }
+
+ // Check for modifier on a per-region level
+ if (surfDict.found("regions"))
+ {
+ const dictionary& regionsDict = surfDict.subDict("regions");
+ forAllConstIter(dictionary, regionsDict, regionsIter)
+ {
+ const dictionary& regionDict = regionsIter().dict();
+ const keyType& regionName = regionsIter().keyword();
+
+ autoPtr<searchableSurfaceModifier> modifier
+ (
+ searchableSurfaceModifier::New
+ (
+ regionDict.lookup("type"),
+ allGeometry,
+ regionDict
+ )
+ );
+
+ labelList regionIDs =
+ findStrings(regionName, surf.regions());
+ if (modifier().modify(regionIDs, surf))
+ {
+ changed = true;
+ }
+ }
+ }
+
+
+ if (changed)
+ {
+ const fileName name(surf.name());
+ surf.rename(name.lessExt() + "_patched." + name.ext());
+
+ Info<< "Writing repatched surface " << name << " to "
+ << surf.name() << nl << endl;
+
+ surf.write();
+ }
+ else
+ {
+ Info<< "Surface " << surf.name() << " unchanged" << nl << endl;
+ }
+ }
+ }
+
+ Info<< "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfacePatch/surfacePatchDict b/applications/utilities/surface/surfacePatch/surfacePatchDict
new file mode 100644
index 0000000000000000000000000000000000000000..5e0f5f19a5b0326643acddc80f8f30a9e95f8819
--- /dev/null
+++ b/applications/utilities/surface/surfacePatch/surfacePatchDict
@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object createPatchDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+geometry
+{
+ box
+ {
+ type searchableRotatedBox;
+
+ // box rotated 45 degrees around z axis
+ span (1 1 1);
+ origin (0 0 0);
+ e1 (1 1 0);
+ e3 (0 0 1);
+ }
+
+ singleTri.obj
+ {
+ type triSurfaceMesh;
+ }
+
+ unitCube_me.stl
+ {
+ type triSurfaceMesh;
+ }
+};
+
+
+
+surfaces
+{
+ unitCube_me.stl
+ {
+ regions
+ {
+ // Divide region 'maxZ' into multiple regions according to
+ // a geometric feature angle
+ maxZ
+ {
+ type autoPatch;
+ featureAngle 60;
+ }
+ }
+ }
+
+ singleTri.obj
+ {
+ // Regionise surface according to triangles (after cutting) being
+ // inside or outside the 'box' geometry
+ type cut;
+ cutters (box);
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C b/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
index 6e4db35a0bdce4829099f02569402298cd78aebf..001745fbca4a24ae2c2e36fc51fe32edcdba5105 100644
--- a/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
+++ b/applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -31,21 +31,17 @@ Description
Note
- best decomposition option is hierarchGeomDecomp since
- guarantees square decompositions.
+ guarantees square decompositions.
- triangles might be present on multiple processors.
- merging uses geometric tolerance so take care with writing precision.
\*---------------------------------------------------------------------------*/
-#include "treeBoundBox.H"
-#include "FixedList.H"
#include "argList.H"
#include "Time.H"
#include "polyMesh.H"
#include "distributedTriSurfaceMesh.H"
#include "mapDistribute.H"
-#include "triSurfaceFields.H"
-#include "Pair.H"
using namespace Foam;
@@ -172,7 +168,7 @@ int main(int argc, char *argv[])
"triSurface", // local
runTime, // registry
IOobject::MUST_READ,
- IOobject::NO_WRITE
+ IOobject::AUTO_WRITE
);
const fileName actualPath(io.filePath());
@@ -215,6 +211,7 @@ int main(int argc, char *argv[])
Info<< "Writing dummy bounds dictionary to " << ioDict.name()
<< nl << endl;
+ // Force writing in ascii
ioDict.regIOobject::writeObject
(
IOstream::ASCII,
@@ -239,23 +236,18 @@ int main(int argc, char *argv[])
(
IOobject
(
- surfMesh.searchableSurface::name(), // name
- surfMesh.searchableSurface::instance(), // instance
+ "faceCentres", // name
+ surfMesh.searchableSurface::time().timeName(), // instance
surfMesh.searchableSurface::local(), // local
surfMesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
surfMesh,
- dimLength
+ dimLength,
+ s.faceCentres()
)
);
- triSurfaceVectorField& fc = fcPtr();
-
- forAll(fc, triI)
- {
- fc[triI] = s[triI].centre(s.points());
- }
// Steal pointer and store object on surfMesh
fcPtr.ptr()->store();
@@ -290,7 +282,7 @@ int main(int argc, char *argv[])
Info<< "Writing surface." << nl << endl;
- surfMesh.searchableSurface::write();
+ surfMesh.objectRegistry::write();
Info<< "End\n" << endl;
diff --git a/applications/utilities/surface/surfaceSubset/surfaceSubset.C b/applications/utilities/surface/surfaceSubset/surfaceSubset.C
index c0ae6bd2418295b4d5dd0e407fb62be4cadff521..c2a168ba68d4ab723c863a464a57acd287483590 100644
--- a/applications/utilities/surface/surfaceSubset/surfaceSubset.C
+++ b/applications/utilities/surface/surfaceSubset/surfaceSubset.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -387,6 +387,8 @@ int main(int argc, char *argv[])
surf2.write(outFileName);
+ Info<< "\nEnd\n" << endl;
+
return 0;
}
diff --git a/src/edgeMesh/edgeMesh.C b/src/edgeMesh/edgeMesh.C
index 0e8cb2acada1f5d7473369d2de1e05ce89bc4c40..947472af78913a592cfdbaeccc402f48105b1a3b 100644
--- a/src/edgeMesh/edgeMesh.C
+++ b/src/edgeMesh/edgeMesh.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -286,7 +286,11 @@ void Foam::edgeMesh::scalePoints(const scalar scaleFactor)
}
-void Foam::edgeMesh::mergePoints(const scalar mergeDist)
+void Foam::edgeMesh::mergePoints
+(
+ const scalar mergeDist,
+ labelList& reversePointMap
+)
{
pointField newPoints;
labelList pointMap;
@@ -307,6 +311,9 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
points_.transfer(newPoints);
+ // connectivity changed
+ pointEdgesPtr_.clear();
+
// Renumber and make sure e[0] < e[1] (not really necessary)
forAll(edges_, edgeI)
{
@@ -383,6 +390,9 @@ void Foam::edgeMesh::mergeEdges()
{
edges_[iter()] = iter.key();
}
+
+ // connectivity changed
+ pointEdgesPtr_.clear();
}
diff --git a/src/edgeMesh/edgeMesh.H b/src/edgeMesh/edgeMesh.H
index d3397a290b0816f233172ec779caef371f01b30c..9962ca9176d2c2d84406cee3350138e59d2901ac 100644
--- a/src/edgeMesh/edgeMesh.H
+++ b/src/edgeMesh/edgeMesh.H
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -247,11 +247,12 @@ public:
//- Scale points. A non-positive factor is ignored
virtual void scalePoints(const scalar);
- //- Merge common points (points within mergeDist)
- void mergePoints(const scalar mergeDist);
+ //- Merge common points (points within mergeDist). Return map from
+ // old to new points.
+ virtual void mergePoints(const scalar mergeDist, labelList&);
- //- Merge similar edges
- void mergeEdges();
+ //- Merge duplicate edges
+ virtual void mergeEdges();
// Write
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
index d7bfa4b2f666a880de031415a151d13eef9e3b9a..f2e379b49376bb7aeb6cc4cf5cc0237e2b9fea9c 100644
--- a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -32,6 +32,8 @@ License
#include "DynamicField.H"
#include "edgeMeshFormatsCore.H"
#include "IOmanip.H"
+#include "searchableSurface.H"
+#include "triSurfaceMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -118,7 +120,6 @@ Foam::wordHashSet Foam::extendedEdgeMesh::writeTypes()
}
-
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
bool Foam::extendedEdgeMesh::canReadType
@@ -176,7 +177,7 @@ Foam::extendedEdgeMesh::classifyFeaturePoint
label ptI
) const
{
- labelList ptEds(pointEdges()[ptI]);
+ const labelList& ptEds(pointEdges()[ptI]);
label nPtEds = ptEds.size();
label nExternal = 0;
@@ -217,46 +218,182 @@ Foam::extendedEdgeMesh::classifyFeaturePoint
}
-Foam::extendedEdgeMesh::edgeStatus
-Foam::extendedEdgeMesh::classifyEdge
+void Foam::extendedEdgeMesh::cut
(
- const List<vector>& norms,
- const labelList& edNorms,
- const vector& fC0tofC1
+ const searchableSurface& surf,
+
+ labelList& pointMap,
+ labelList& edgeMap,
+ labelList& pointsFromEdge,
+ labelList& oldEdge,
+ labelList& surfTri
)
{
- label nEdNorms = edNorms.size();
+ const edgeList& edges = this->edges();
+ const pointField& points = this->points();
- if (nEdNorms == 1)
+
+ List<List<pointIndexHit> > edgeHits(edges.size());
{
- return OPEN;
+ pointField start(edges.size());
+ pointField end(edges.size());
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+ start[edgeI] = points[e[0]];
+ end[edgeI] = points[e[1]];
+ }
+ surf.findLineAll(start, end, edgeHits);
}
- else if (nEdNorms == 2)
+
+ // Count number of hits
+ label nHits = 0;
+ forAll(edgeHits, edgeI)
{
- const vector& n0(norms[edNorms[0]]);
- const vector& n1(norms[edNorms[1]]);
+ nHits += edgeHits[edgeI].size();
+ }
- if ((n0 & n1) > cosNormalAngleTol_)
+ DynamicField<point> newPoints(points);
+ DynamicList<label> newToOldPoint(identity(points.size()));
+
+ newPoints.setCapacity(newPoints.size()+nHits);
+ newToOldPoint.setCapacity(newPoints.capacity());
+
+ DynamicList<edge> newEdges(edges);
+ DynamicList<label> newToOldEdge(identity(edges.size()));
+
+ newEdges.setCapacity(newEdges.size()+nHits);
+ newToOldEdge.setCapacity(newEdges.capacity());
+
+ // Information on additional points
+ DynamicList<label> dynPointsFromEdge(nHits);
+ DynamicList<label> dynOldEdge(nHits);
+ DynamicList<label> dynSurfTri(nHits);
+
+ forAll(edgeHits, edgeI)
+ {
+ const List<pointIndexHit>& eHits = edgeHits[edgeI];
+
+ if (eHits.size())
{
- return FLAT;
+ label prevPtI = edges[edgeI][0];
+ forAll(eHits, eHitI)
+ {
+ label newPtI = newPoints.size();
+
+ newPoints.append(eHits[eHitI].hitPoint());
+ newToOldPoint.append(edges[edgeI][0]); // map from start point
+ dynPointsFromEdge.append(newPtI);
+ dynOldEdge.append(edgeI);
+ dynSurfTri.append(eHits[eHitI].index());
+
+ if (eHitI == 0)
+ {
+ newEdges[edgeI] = edge(prevPtI, newPtI);
+ }
+ else
+ {
+ newEdges.append(edge(prevPtI, newPtI));
+ newToOldEdge.append(edgeI);
+ }
+ prevPtI = newPtI;
+ }
+ newEdges.append(edge(prevPtI, edges[edgeI][1]));
+ newToOldEdge.append(edgeI);
}
- else if ((fC0tofC1 & n0) > 0.0)
+ }
+
+ pointField allPoints;
+ allPoints.transfer(newPoints);
+ pointMap.transfer(newToOldPoint);
+
+ edgeList allEdges;
+ allEdges.transfer(newEdges);
+ edgeMap.transfer(newToOldEdge);
+
+ pointsFromEdge.transfer(dynPointsFromEdge);
+ oldEdge.transfer(dynOldEdge);
+ surfTri.transfer(dynSurfTri);
+
+ // Update local information
+ autoMap(allPoints, allEdges, pointMap, edgeMap);
+}
+
+
+void Foam::extendedEdgeMesh::select
+(
+ const searchableSurface& surf,
+ const volumeType volType, // side to keep
+ labelList& pointMap, // sub to old points
+ labelList& edgeMap // sub to old edges
+)
+{
+ const edgeList& edges = this->edges();
+ const pointField& points = this->points();
+
+ // Test edge centres for inside/outside
+ if (volType == volumeType::INSIDE || volType == volumeType::OUTSIDE)
+ {
+ pointField edgeCentres(edges.size());
+ forAll(edgeCentres, edgeI)
{
- return INTERNAL;
+ const edge& e = edges[edgeI];
+ edgeCentres[edgeI] = e.centre(points);
}
- else
+ List<volumeType> volTypes;
+ surf.getVolumeType(edgeCentres, volTypes);
+
+ // Extract edges on correct side
+ edgeMap.setSize(edges.size());
+ label compactEdgeI = 0;
+
+ forAll(volTypes, edgeI)
{
- return EXTERNAL;
+ if (volTypes[edgeI] == volType)
+ {
+ edgeMap[compactEdgeI++] = edgeI;
+ }
}
- }
- else if (nEdNorms > 2)
- {
- return MULTIPLE;
+ edgeMap.setSize(compactEdgeI);
+
+ // Extract used points
+ labelList pointToCompact(points.size(), -1);
+ forAll(edgeMap, i)
+ {
+ const edge& e = edges[edgeMap[i]];
+ pointToCompact[e[0]] = labelMax; // tag with a value
+ pointToCompact[e[1]] = labelMax;
+ }
+
+ pointMap.setSize(points.size());
+ label compactPointI = 0;
+ forAll(pointToCompact, pointI)
+ {
+ if (pointToCompact[pointI] != -1)
+ {
+ pointToCompact[pointI] = compactPointI;
+ pointMap[compactPointI++] = pointI;
+ }
+ }
+ pointMap.setSize(compactPointI);
+ pointField subPoints(points, pointMap);
+
+ // Renumber edges
+ edgeList subEdges(edgeMap.size());
+ forAll(edgeMap, i)
+ {
+ const edge& e = edges[edgeMap[i]];
+ subEdges[i][0] = pointToCompact[e[0]];
+ subEdges[i][1] = pointToCompact[e[1]];
+ }
+
+ // Reset primitives and map other quantities
+ autoMap(subPoints, subEdges, pointMap, edgeMap);
}
else
{
- // There is a problem - the edge has no normals
- return NONE;
+ pointMap = identity(points.size());
+ edgeMap = identity(edges.size());
}
}
@@ -287,10 +424,7 @@ Foam::extendedEdgeMesh::extendedEdgeMesh()
{}
-Foam::extendedEdgeMesh::extendedEdgeMesh
-(
- const extendedEdgeMesh& fem
-)
+Foam::extendedEdgeMesh::extendedEdgeMesh(const extendedEdgeMesh& fem)
:
edgeMesh(fem),
concaveStart_(fem.concaveStart()),
@@ -1332,126 +1466,649 @@ void Foam::extendedEdgeMesh::flipNormals()
}
-void Foam::extendedEdgeMesh::writeObj
+void Foam::extendedEdgeMesh::autoMap
(
- const fileName& prefix
-) const
+ const pointField& subPoints,
+ const edgeList& subEdges,
+ const labelList& pointMap,
+ const labelList& edgeMap
+)
{
- Info<< nl << "Writing extendedEdgeMesh components to " << prefix
- << endl;
+ // Determine slicing for subEdges
+ label subIntStart = edgeMap.size();
+ label subFlatStart = edgeMap.size();
+ label subOpenStart = edgeMap.size();
+ label subMultipleStart = edgeMap.size();
+
+ forAll(edgeMap, subEdgeI)
+ {
+ label edgeI = edgeMap[subEdgeI];
+ if (edgeI >= internalStart() && subIntStart == edgeMap.size())
+ {
+ subIntStart = subEdgeI;
+ }
+ if (edgeI >= flatStart() && subFlatStart == edgeMap.size())
+ {
+ subFlatStart = subEdgeI;
+ }
+ if (edgeI >= openStart() && subOpenStart == edgeMap.size())
+ {
+ subOpenStart = subEdgeI;
+ }
+ if (edgeI >= multipleStart() && subMultipleStart == edgeMap.size())
+ {
+ subMultipleStart = subEdgeI;
+ }
+ }
+
+
+ // Determine slicing for subPoints
+
+ label subConcaveStart = pointMap.size();
+ label subMixedStart = pointMap.size();
+ label subNonFeatStart = pointMap.size();
+
+ forAll(pointMap, subPointI)
+ {
+ label pointI = pointMap[subPointI];
+ if (pointI >= concaveStart() && subConcaveStart == pointMap.size())
+ {
+ subConcaveStart = subPointI;
+ }
+ if (pointI >= mixedStart() && subMixedStart == pointMap.size())
+ {
+ subMixedStart = subPointI;
+ }
+ if
+ (
+ pointI >= nonFeatureStart()
+ && subNonFeatStart == pointMap.size()
+ )
+ {
+ subNonFeatStart = subPointI;
+ }
+ }
- edgeMesh::write(prefix + "_edgeMesh.obj");
- OBJstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
- Info<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
- for(label i = 0; i < concaveStart_; i++)
+ // Compact region edges
+ labelList subRegionEdges;
{
- convexFtPtStr.write(points()[i]);
+ PackedBoolList isRegionEdge(edges().size());
+ forAll(regionEdges(), i)
+ {
+ label edgeI = regionEdges()[i];
+ isRegionEdge[edgeI] = true;
+ }
+
+ DynamicList<label> newRegionEdges(regionEdges().size());
+ forAll(edgeMap, subEdgeI)
+ {
+ label edgeI = edgeMap[subEdgeI];
+ if (isRegionEdge[edgeI])
+ {
+ newRegionEdges.append(subEdgeI);
+ }
+ }
+ subRegionEdges.transfer(newRegionEdges);
}
- OBJstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
- Info<< "Writing concave feature points to "
- << concaveFtPtStr.name() << endl;
- for(label i = concaveStart_; i < mixedStart_; i++)
+ labelListList subFeaturePointEdges;
+ if (featurePointEdges().size())
{
- convexFtPtStr.write(points()[i]);
+ subFeaturePointEdges.setSize(subNonFeatStart);
+ for (label subPointI = 0; subPointI < subNonFeatStart; subPointI++)
+ {
+ label pointI = pointMap[subPointI];
+ const labelList& pEdges = featurePointEdges()[pointI];
+
+ labelList& subPEdges = subFeaturePointEdges[subPointI];
+ subPEdges.setSize(pEdges.size());
+
+ if (pEdges.size())
+ {
+ forAll(pEdges, i)
+ {
+ subPEdges[i] = edgeMap[pEdges[i]];
+ }
+ }
+ }
}
- OBJstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
- Info<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
- for(label i = mixedStart_; i < nonFeatureStart_; i++)
+ vectorField subEdgeDirections(edgeDirections(), edgeMap);
+
+ // Find used normals
+ labelList reverseNormalMap(normals().size(), -1);
+ DynamicList<label> normalMap(normals().size());
+
{
- mixedFtPtStr.write(points()[i]);
+ PackedBoolList isSubNormal(normals().size());
+ for (label subPointI = 0; subPointI < subNonFeatStart; subPointI++)
+ {
+ label pointI = pointMap[subPointI];
+ const labelList& pNormals = featurePointNormals()[pointI];
+
+ forAll(pNormals, i)
+ {
+ isSubNormal[pNormals[i]] = true;
+ }
+ }
+ forAll(edgeMap, subEdgeI)
+ {
+ label edgeI = edgeMap[subEdgeI];
+ const labelList& eNormals = edgeNormals()[edgeI];
+
+ forAll(eNormals, i)
+ {
+ isSubNormal[eNormals[i]] = true;
+ }
+ }
+
+ forAll(isSubNormal, normalI)
+ {
+ if (isSubNormal[normalI])
+ {
+ label subNormalI = normalMap.size();
+ reverseNormalMap[normalI] = subNormalI;
+ normalMap.append(subNormalI);
+ }
+ }
}
- OBJstream mixedFtPtStructureStr(prefix + "_mixedFeaturePtsStructure.obj");
- Info<< "Writing mixed feature point structure to "
- << mixedFtPtStructureStr.name() << endl;
- for(label i = mixedStart_; i < nonFeatureStart_; i++)
+ // Use compaction map on data referencing normals
+ labelListList subNormalDirections;
+
+ if (normalDirections().size())
{
- const labelList& ptEds = pointEdges()[i];
+ subNormalDirections.setSize(edgeMap.size());
- forAll(ptEds, j)
+ forAll(edgeMap, subEdgeI)
{
- const edge& e = edges()[ptEds[j]];
- mixedFtPtStructureStr.write
+ label edgeI = edgeMap[subEdgeI];
+ const labelList& eNormals = normalDirections()[edgeI];
+
+ labelList& subNormals = subNormalDirections[subEdgeI];
+ subNormals.setSize(eNormals.size());
+ forAll(eNormals, i)
+ {
+ if (eNormals[i] >= 0)
+ {
+ subNormals[i] = reverseNormalMap[eNormals[i]];
+ }
+ else
+ {
+ subNormals[i] = -1;
+ }
+ }
+ }
+ }
+
+ labelListList subEdgeNormals(edgeMap.size());
+ forAll(edgeMap, subEdgeI)
+ {
+ label edgeI = edgeMap[subEdgeI];
+ const labelList& eNormals = edgeNormals()[edgeI];
+ labelList& subNormals = subEdgeNormals[subEdgeI];
+
+ subNormals = UIndirectList<label>(reverseNormalMap, eNormals);
+ }
+
+ labelListList subPointNormals(pointMap.size());
+ for (label subPointI = 0; subPointI < subNonFeatStart; subPointI++)
+ {
+ label pointI = pointMap[subPointI];
+ const labelList& pNormals = featurePointNormals()[pointI];
+ labelList& subNormals = subPointNormals[subPointI];
+
+ subNormals = UIndirectList<label>(reverseNormalMap, pNormals);
+ }
+
+ // Use compaction map to compact normal data
+ vectorField subNormals(normals(), normalMap);
+
+ List<extendedEdgeMesh::sideVolumeType> subNormalVolumeTypes;
+ if (normalVolumeTypes().size())
+ {
+ subNormalVolumeTypes =
+ UIndirectList<extendedEdgeMesh::sideVolumeType>
(
- linePointRef(points()[e[0]],
- points()[e[1]])
+ normalVolumeTypes(),
+ normalMap
);
+ }
+
+ extendedEdgeMesh subMesh
+ (
+ subPoints,
+ subEdges,
+
+ // Feature points slices
+ subConcaveStart,
+ subMixedStart,
+ subNonFeatStart,
+
+ // Feature edges slices
+ subIntStart,
+ subFlatStart,
+ subOpenStart,
+ subMultipleStart,
+
+ // All normals
+ subNormals,
+ subNormalVolumeTypes,
+
+ // Per edge edge vector
+ subEdgeDirections,
+
+ // Per edge list of normal indices
+ subNormalDirections,
+ // Per edge list of normal indices
+ subEdgeNormals,
+
+ // Per point list of normal indices
+ subPointNormals,
+ subFeaturePointEdges,
+ subRegionEdges
+ );
+
+ transfer(subMesh);
+}
+
+
+void Foam::extendedEdgeMesh::trim
+(
+ const searchableSurface& surf,
+ const volumeType volType,
+ labelList& pointMap,
+ labelList& edgeMap
+)
+{
+ // Cut edges with the other surfaces
+
+ labelList allPointMap; // from all to original point
+ labelList allEdgeMap; // from all to original edge
+
+ labelList pointsFromEdge; // list of new points created by cutting
+ labelList oldEdge; // for each of these points the orginal edge
+ labelList surfTri; // for each of these points the surface triangle
+ cut
+ (
+ surf,
+ allPointMap,
+ allEdgeMap,
+ pointsFromEdge,
+ oldEdge,
+ surfTri
+ );
+
+ const label nOldPoints = points().size();
+
+ // Remove outside edges and compact
+
+ labelList subPointMap; // sub to old points
+ labelList subEdgeMap; // sub to old edges
+ select(surf, volType, subPointMap, subEdgeMap);
+
+ // Update overall point maps
+ pointMap = UIndirectList<label>(allPointMap, subPointMap);
+ edgeMap = UIndirectList<label>(allEdgeMap, subEdgeMap);
+
+ // Extract current point and edge status
+ List<edgeStatus> edgeStat(edges().size());
+ List<pointStatus> pointStat(points().size());
+ forAll(edgeStat, edgeI)
+ {
+ edgeStat[edgeI] = getEdgeStatus(edgeI);
+ }
+ forAll(pointStat, pointI)
+ {
+ pointStat[pointI] = getPointStatus(pointI);
+ }
+
+ // Re-classify exposed points (from cutting)
+ labelList oldPointToIndex(nOldPoints, -1);
+ forAll(pointsFromEdge, i)
+ {
+ oldPointToIndex[pointsFromEdge[i]] = i;
+ }
+ forAll(subPointMap, pointI)
+ {
+ label oldPointI = subPointMap[pointI];
+ label index = oldPointToIndex[oldPointI];
+ if (index != -1)
+ {
+ pointStat[pointI] = classifyFeaturePoint(pointI);
}
}
- OBJstream externalStr(prefix + "_externalEdges.obj");
- Info<< "Writing external edges to " << externalStr.name() << endl;
+ // Reset based on new point and edge status
+ labelList sortedToOriginalPoint;
+ labelList sortedToOriginalEdge;
+ setFromStatus
+ (
+ pointStat,
+ edgeStat,
+ sortedToOriginalPoint,
+ sortedToOriginalEdge
+ );
+
+ // Update the overall pointMap, edgeMap
+ pointMap = UIndirectList<label>(pointMap, sortedToOriginalPoint)();
+ edgeMap = UIndirectList<label>(edgeMap, sortedToOriginalEdge)();
+}
+
+
+void Foam::extendedEdgeMesh::setFromStatus
+(
+ const List<extendedEdgeMesh::pointStatus>& pointStat,
+ const List<extendedEdgeMesh::edgeStatus>& edgeStat,
+ labelList& sortedToOriginalPoint,
+ labelList& sortedToOriginalEdge
+)
+{
+ // Use pointStatus and edgeStatus to determine new ordering
+ label pointConcaveStart;
+ label pointMixedStart;
+ label pointNonFeatStart;
+
+ label edgeInternalStart;
+ label edgeFlatStart;
+ label edgeOpenStart;
+ label edgeMultipleStart;
+ sortedOrder
+ (
+ pointStat,
+ edgeStat,
+ sortedToOriginalPoint,
+ sortedToOriginalEdge,
+
+ pointConcaveStart,
+ pointMixedStart,
+ pointNonFeatStart,
+
+ edgeInternalStart,
+ edgeFlatStart,
+ edgeOpenStart,
+ edgeMultipleStart
+ );
+
- for (label i = externalStart_; i < internalStart_; i++)
+ // Order points and edges
+ labelList reversePointMap(points().size(), -1);
+ forAll(sortedToOriginalPoint, sortedI)
{
- const edge& e = edges()[i];
- externalStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ reversePointMap[sortedToOriginalPoint[sortedI]] = sortedI;
}
- OBJstream internalStr(prefix + "_internalEdges.obj");
- Info<< "Writing internal edges to " << internalStr.name() << endl;
+ edgeList sortedEdges(UIndirectList<edge>(edges(), sortedToOriginalEdge)());
+ forAll(sortedEdges, sortedI)
+ {
+ inplaceRenumber(reversePointMap, sortedEdges[sortedI]);
+ }
+
+ // Update local data
+ autoMap
+ (
+ pointField(points(), sortedToOriginalPoint),
+ sortedEdges,
+ sortedToOriginalPoint,
+ sortedToOriginalEdge
+ );
- for (label i = internalStart_; i < flatStart_; i++)
+ // Reset the slice starts
+ concaveStart_ = pointConcaveStart;
+ mixedStart_ = pointMixedStart;
+ nonFeatureStart_ = pointNonFeatStart;
+ internalStart_ = edgeInternalStart;
+ flatStart_ = edgeFlatStart;
+ openStart_ = edgeOpenStart;
+ multipleStart_ = edgeMultipleStart;
+}
+
+
+bool Foam::extendedEdgeMesh::mergePointsAndSort
+(
+ const scalar mergeDist,
+ labelList& pointMap,
+ labelList& edgeMap
+)
+{
+ const label nOldPoints = points().size();
+
+ // Detect and merge collocated feature points
+ labelList oldToMerged;
+ label nNewPoints = ::Foam::mergePoints
+ (
+ points(),
+ SMALL,
+ false,
+ oldToMerged
+ );
+
+ pointMap.setSize(nNewPoints);
+ pointMap = -1;
+ forAll(oldToMerged, oldI)
{
- const edge& e = edges()[i];
- internalStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ label newI = oldToMerged[oldI];
+ if (pointMap[newI] == -1)
+ {
+ pointMap[newI] = oldI;
+ }
}
- OBJstream flatStr(prefix + "_flatEdges.obj");
- Info<< "Writing flat edges to " << flatStr.name() << endl;
+ // Renumber edges
+ edgeList newEdges(edges().size());
+ forAll(edges(), edgeI)
+ {
+ const edge& oldE = edges()[edgeI];
+ newEdges[edgeI] = edge(oldToMerged[oldE[0]], oldToMerged[oldE[1]]);
+ }
+
+ // Shuffle basic information (reorders point data)
+ autoMap
+ (
+ pointField(points(), pointMap),
+ newEdges,
+ pointMap,
+ identity(newEdges.size())
+ );
- for (label i = flatStart_; i < openStart_; i++)
+ // Re-classify the merged points
+ List<edgeStatus> edgeStat(edges().size());
+ forAll(edgeStat, edgeI)
{
- const edge& e = edges()[i];
- flatStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ edgeStat[edgeI] = getEdgeStatus(edgeI);
}
- OBJstream openStr(prefix + "_openEdges.obj");
- Info<< "Writing open edges to " << openStr.name() << endl;
+ List<pointStatus> pointStat(points().size());
+ forAll(pointStat, pointI)
+ {
+ pointStat[pointI] = getPointStatus(pointI);
+ }
- for (label i = openStart_; i < multipleStart_; i++)
+ // Re-classify merged points
+ labelList nPoints(nNewPoints, 0);
+ forAll(oldToMerged, oldPointI)
{
- const edge& e = edges()[i];
- openStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ nPoints[oldToMerged[oldPointI]]++;
}
- OBJstream multipleStr(prefix + "_multipleEdges.obj");
- Info<< "Writing multiple edges to " << multipleStr.name() << endl;
+ forAll(nPoints, pointI)
+ {
+ if (nPoints[pointI] != 1)
+ {
+ pointStat[pointI] = classifyFeaturePoint(pointI);
+ }
+ }
+
+ labelList sortedToOriginalPoint;
+ setFromStatus
+ (
+ pointStat,
+ edgeStat,
+ sortedToOriginalPoint,
+ edgeMap // point merging above did not affect edge order
+ );
+ pointMap = UIndirectList<label>(pointMap, sortedToOriginalPoint)();
+
+ return nNewPoints != nOldPoints;
+}
+
+
+void Foam::extendedEdgeMesh::writeObj(const fileName& prefix) const
+{
+ Info<< nl << "Writing extendedEdgeMesh components to " << prefix
+ << endl;
+
+ edgeMesh::write(prefix + "_edgeMesh.obj");
- for (label i = multipleStart_; i < edges().size(); i++)
{
- const edge& e = edges()[i];
- multipleStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ OBJstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
+ Info<< "Writing " << concaveStart_
+ << " convex feature points to " << convexFtPtStr.name() << endl;
+
+ for(label i = 0; i < concaveStart_; i++)
+ {
+ convexFtPtStr.write(points()[i]);
+ }
+ }
+
+ {
+ OBJstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
+ Info<< "Writing " << mixedStart_-concaveStart_
+ << " concave feature points to "
+ << concaveFtPtStr.name() << endl;
+
+ for(label i = concaveStart_; i < mixedStart_; i++)
+ {
+ concaveFtPtStr.write(points()[i]);
+ }
+ }
+
+ {
+ OBJstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
+ Info<< "Writing " << nonFeatureStart_-mixedStart_
+ << " mixed feature points to " << mixedFtPtStr.name() << endl;
+
+ for(label i = mixedStart_; i < nonFeatureStart_; i++)
+ {
+ mixedFtPtStr.write(points()[i]);
+ }
}
- OBJstream regionStr(prefix + "_regionEdges.obj");
- Info<< "Writing region edges to " << regionStr.name() << endl;
+ {
+ OBJstream mixedFtPtStructureStr(prefix+"_mixedFeaturePtsStructure.obj");
+ Info<< "Writing "
+ << nonFeatureStart_-mixedStart_
+ << " mixed feature point structure to "
+ << mixedFtPtStructureStr.name() << endl;
+
+ for(label i = mixedStart_; i < nonFeatureStart_; i++)
+ {
+ const labelList& ptEds = pointEdges()[i];
+
+ forAll(ptEds, j)
+ {
+ const edge& e = edges()[ptEds[j]];
+ mixedFtPtStructureStr.write
+ (
+ linePointRef(points()[e[0]],
+ points()[e[1]])
+ );
+ }
+ }
+ }
- forAll(regionEdges_, i)
{
- const edge& e = edges()[regionEdges_[i]];
- regionStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ OBJstream externalStr(prefix + "_externalEdges.obj");
+ Info<< "Writing " << internalStart_-externalStart_
+ << " external edges to " << externalStr.name() << endl;
+
+ for (label i = externalStart_; i < internalStart_; i++)
+ {
+ const edge& e = edges()[i];
+ externalStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ }
}
- OBJstream edgeDirsStr(prefix + "_edgeDirections.obj");
- Info<< "Writing edge directions to " << edgeDirsStr.name() << endl;
+ {
+ OBJstream internalStr(prefix + "_internalEdges.obj");
+ Info<< "Writing " << flatStart_-internalStart_
+ << " internal edges to " << internalStr.name() << endl;
+
+ for (label i = internalStart_; i < flatStart_; i++)
+ {
+ const edge& e = edges()[i];
+ internalStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ }
+ }
- forAll(edgeDirections_, i)
{
- const vector& eVec = edgeDirections_[i];
- const edge& e = edges()[i];
+ OBJstream flatStr(prefix + "_flatEdges.obj");
+ Info<< "Writing " << openStart_-flatStart_
+ << " flat edges to " << flatStr.name() << endl;
- edgeDirsStr.write
- (
- linePointRef(points()[e.start()], eVec + points()[e.start()])
- );
+ for (label i = flatStart_; i < openStart_; i++)
+ {
+ const edge& e = edges()[i];
+ flatStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ }
+ }
+
+ {
+ OBJstream openStr(prefix + "_openEdges.obj");
+ Info<< "Writing " << multipleStart_-openStart_
+ << " open edges to " << openStr.name() << endl;
+
+ for (label i = openStart_; i < multipleStart_; i++)
+ {
+ const edge& e = edges()[i];
+ openStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ }
+ }
+
+ {
+ OBJstream multipleStr(prefix + "_multipleEdges.obj");
+ Info<< "Writing " << edges().size()-multipleStart_
+ << " multiple edges to " << multipleStr.name() << endl;
+
+ for (label i = multipleStart_; i < edges().size(); i++)
+ {
+ const edge& e = edges()[i];
+ multipleStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ }
+ }
+
+ {
+ OBJstream regionStr(prefix + "_regionEdges.obj");
+ Info<< "Writing " << regionEdges_.size()
+ << " region edges to " << regionStr.name() << endl;
+
+ forAll(regionEdges_, i)
+ {
+ const edge& e = edges()[regionEdges_[i]];
+ regionStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+ }
+ }
+
+ {
+ OBJstream edgeDirsStr(prefix + "_edgeDirections.obj");
+ Info<< "Writing " << edgeDirections_.size()
+ << " edge directions to " << edgeDirsStr.name() << endl;
+
+ forAll(edgeDirections_, i)
+ {
+ const vector& eVec = edgeDirections_[i];
+ const edge& e = edges()[i];
+
+ edgeDirsStr.write
+ (
+ linePointRef(points()[e.start()], eVec + points()[e.start()])
+ );
+ }
}
}
@@ -1506,6 +2163,228 @@ void Foam::extendedEdgeMesh::writeStats(Ostream& os) const
}
+Foam::extendedEdgeMesh::edgeStatus
+Foam::extendedEdgeMesh::classifyEdge
+(
+ const List<vector>& norms,
+ const labelList& edNorms,
+ const vector& fC0tofC1
+)
+{
+ label nEdNorms = edNorms.size();
+
+ if (nEdNorms == 1)
+ {
+ return OPEN;
+ }
+ else if (nEdNorms == 2)
+ {
+ const vector& n0(norms[edNorms[0]]);
+ const vector& n1(norms[edNorms[1]]);
+
+ if ((n0 & n1) > cosNormalAngleTol_)
+ {
+ return FLAT;
+ }
+ else if ((fC0tofC1 & n0) > 0.0)
+ {
+ return INTERNAL;
+ }
+ else
+ {
+ return EXTERNAL;
+ }
+ }
+ else if (nEdNorms > 2)
+ {
+ return MULTIPLE;
+ }
+ else
+ {
+ // There is a problem - the edge has no normals
+ return NONE;
+ }
+}
+
+
+void Foam::extendedEdgeMesh::sortedOrder
+(
+ const List<extendedEdgeMesh::pointStatus>& pointStat,
+ const List<extendedEdgeMesh::edgeStatus>& edgeStat,
+ labelList& sortedToOriginalPoint,
+ labelList& sortedToOriginalEdge,
+
+ label& pointConcaveStart,
+ label& pointMixedStart,
+ label& pointNonFeatStart,
+
+ label& edgeInternalStart,
+ label& edgeFlatStart,
+ label& edgeOpenStart,
+ label& edgeMultipleStart
+)
+{
+ sortedToOriginalPoint.setSize(pointStat.size());
+ sortedToOriginalPoint = -1;
+
+ sortedToOriginalEdge.setSize(edgeStat.size());
+ sortedToOriginalEdge = -1;
+
+
+ // Order edges
+ // ~~~~~~~~~~~
+
+ label nConvex = 0;
+ label nConcave = 0;
+ label nMixed = 0;
+ label nNonFeat = 0;
+
+ forAll(pointStat, pointI)
+ {
+ switch (pointStat[pointI])
+ {
+ case extendedEdgeMesh::CONVEX:
+ nConvex++;
+ break;
+
+ case extendedEdgeMesh::CONCAVE:
+ nConcave++;
+ break;
+
+ case extendedEdgeMesh::MIXED:
+ nMixed++;
+ break;
+
+ case extendedEdgeMesh::NONFEATURE:
+ nNonFeat++;
+ break;
+
+ default:
+ FatalErrorIn("order(..)") << "Problem" << exit(FatalError);
+ break;
+ }
+ }
+
+ label convexStart = 0;
+ label concaveStart = nConvex;
+ label mixedStart = concaveStart+nConcave;
+ label nonFeatStart = mixedStart+nMixed;
+
+
+ // Copy to parameters
+ pointConcaveStart = concaveStart;
+ pointMixedStart = mixedStart;
+ pointNonFeatStart = nonFeatStart;
+
+ forAll(pointStat, pointI)
+ {
+ switch (pointStat[pointI])
+ {
+ case extendedEdgeMesh::CONVEX:
+ sortedToOriginalPoint[convexStart++] = pointI;
+ break;
+
+ case extendedEdgeMesh::CONCAVE:
+ sortedToOriginalPoint[concaveStart++] = pointI;
+ break;
+
+ case extendedEdgeMesh::MIXED:
+ sortedToOriginalPoint[mixedStart++] = pointI;
+ break;
+
+ case extendedEdgeMesh::NONFEATURE:
+ sortedToOriginalPoint[nonFeatStart++] = pointI;
+ break;
+ }
+ }
+
+
+ // Order edges
+ // ~~~~~~~~~~~
+
+ label nExternal = 0;
+ label nInternal = 0;
+ label nFlat = 0;
+ label nOpen = 0;
+ label nMultiple = 0;
+
+ forAll(edgeStat, edgeI)
+ {
+ switch (edgeStat[edgeI])
+ {
+ case extendedEdgeMesh::EXTERNAL:
+ nExternal++;
+ break;
+
+ case extendedEdgeMesh::INTERNAL:
+ nInternal++;
+ break;
+
+ case extendedEdgeMesh::FLAT:
+ nFlat++;
+ break;
+
+ case extendedEdgeMesh::OPEN:
+ nOpen++;
+ break;
+
+ case extendedEdgeMesh::MULTIPLE:
+ nMultiple++;
+ break;
+
+ case extendedEdgeMesh::NONE:
+ default:
+ FatalErrorIn("order(..)") << "Problem" << exit(FatalError);
+ break;
+ }
+ }
+
+ label externalStart = 0;
+ label internalStart = nExternal;
+ label flatStart = internalStart + nInternal;
+ label openStart = flatStart + nFlat;
+ label multipleStart = openStart + nOpen;
+
+
+ // Copy to parameters
+ edgeInternalStart = internalStart;
+ edgeFlatStart = flatStart;
+ edgeOpenStart = openStart;
+ edgeMultipleStart = multipleStart;
+
+ forAll(edgeStat, edgeI)
+ {
+ switch (edgeStat[edgeI])
+ {
+ case extendedEdgeMesh::EXTERNAL:
+ sortedToOriginalEdge[externalStart++] = edgeI;
+ break;
+
+ case extendedEdgeMesh::INTERNAL:
+ sortedToOriginalEdge[internalStart++] = edgeI;
+ break;
+
+ case extendedEdgeMesh::FLAT:
+ sortedToOriginalEdge[flatStart++] = edgeI;
+ break;
+
+ case extendedEdgeMesh::OPEN:
+ sortedToOriginalEdge[openStart++] = edgeI;
+ break;
+
+ case extendedEdgeMesh::MULTIPLE:
+ sortedToOriginalEdge[multipleStart++] = edgeI;
+ break;
+
+ case extendedEdgeMesh::NONE:
+ default:
+ FatalErrorIn("order(..)") << "Problem" << exit(FatalError);
+ break;
+ }
+ }
+}
+
+
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Istream& Foam::operator>>
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
index 58e2ba8cabe922298fa2ef03f171861d60653e34..1dfb921e489f3eb64c6cc969b9bae4d3d4a51c94 100644
--- a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -68,7 +68,7 @@ namespace Foam
{
class surfaceFeatures;
-class objectRegistry;
+class searchableSurface;
// Forward declaration of friend functions and operators
class extendedEdgeMesh;
@@ -138,7 +138,7 @@ protected:
static label externalStart_;
- // Private data
+ // Protected data
//- Index of the start of the concave feature points
label concaveStart_;
@@ -200,12 +200,33 @@ protected:
mutable PtrList<indexedOctree<treeDataEdge> > edgeTreesByType_;
- // Private Member Functions
+ // Protected Member Functions
- //- Classify the type of feature point. Requires valid stored member
+ //- Classify the type of feature point. Requires valid stored member
// data for edges and normals.
pointStatus classifyFeaturePoint(label ptI) const;
+ //- Cut edges with surface. Return map from cut points&edges back
+ // to original
+ void cut
+ (
+ const searchableSurface&,
+ labelList& pMap,
+ labelList& eMap,
+ labelList& pointsFromEdge, // new points created by cutting
+ labelList& oldEdge, // the orginal edge
+ labelList& surfTri // the surface triangle index
+ );
+
+ //- Remove outside/inside edges. volType denotes which side to keep
+ void select
+ (
+ const searchableSurface& surf,
+ const volumeType volType,
+ labelList& pMap,
+ labelList& eMap
+ );
+
template<class Patch>
void sortPointsAndEdges
(
@@ -262,9 +283,8 @@ public:
const Xfer<edgeList>&
);
- //- Construct given a surface with selected edges,point
- // (surfaceFeatures), an objectRegistry and a
- // fileName to write to.
+ //- Construct given a surface with selected edges,points
+ // (surfaceFeatures)
// Extracts, classifies and reorders the data from surfaceFeatures.
extendedEdgeMesh
(
@@ -501,6 +521,43 @@ public:
// etc.
void flipNormals();
+ //- Update with derived geometry
+ void autoMap
+ (
+ const pointField&,
+ const edgeList&,
+ const labelList& pointMap,
+ const labelList& edgeMap
+ );
+
+ //- Trim to surface. Keep volType side. Return map from current back
+ // to original points (-1 for newly introduced points), edges
+ void trim
+ (
+ const searchableSurface& surf,
+ const volumeType volType,
+ labelList& pointMap,
+ labelList& edgeMap
+ );
+
+ //- Order according to point and edge status
+ void setFromStatus
+ (
+ const List<extendedEdgeMesh::pointStatus>& pointStat,
+ const List<extendedEdgeMesh::edgeStatus>& edgeStat,
+ labelList& sortedToOriginalPoint,
+ labelList& sortedToOriginalEdge
+ );
+
+ //- Geometric merge points. Returns true if any points merged.
+ // Return maps from new back to original points/edges.
+ bool mergePointsAndSort
+ (
+ const scalar mergeDist,
+ labelList& pointMap,
+ labelList& edgeMap
+ );
+
// Read
@@ -532,6 +589,24 @@ public:
const vector& fC0tofC1
);
+ //- Determine the ordering
+ static void sortedOrder
+ (
+ const List<extendedEdgeMesh::pointStatus>& pointStat,
+ const List<extendedEdgeMesh::edgeStatus>& edgeStat,
+ labelList& sortedToOriginalPoint,
+ labelList& sortedToOriginalEdge,
+
+ label& pointConcaveStart,
+ label& pointMixedStart,
+ label& pointNonFeatStart,
+
+ label& edgeInternalStart,
+ label& edgeFlatStart,
+ label& edgeOpenStart,
+ label& edgeMultipleStart
+ );
+
// Ostream Operator