zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-12-28 03:37:59 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'cvMesh'

Browse Source
This commit is contained in:
laurence
2012-02-07 12:12:12 +00:00
parent 0dbdea4781 3f23f5f40f
commit 35fdacd810
42 changed files with 3748 additions and 1957 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/files
View File

@ -4,6 +4,7 @@ conformalVoronoiMesh/conformalVoronoiMesh.C
conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
conformalVoronoiMesh/conformalVoronoiMeshIO.C
conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
cvControls/cvControls.C

3
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
View File

@ -14,7 +14,8 @@ EXE_INC = \
-I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
-I$(LIB_SRC)/edgeMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude
-I$(LIB_SRC)/triSurface/lnInclude \
-I../vectorTools
EXE_LIBS = \
-lmeshTools \

42
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3Ddefs.H
View File

@ -72,6 +72,48 @@ typedef Delaunay::Cell_handle Cell_handle;
typedef Delaunay::Point Point;
//- Spatial sort traits to use with a pair of point pointers and an integer.
// Taken from a post on the CGAL lists: 2010-01/msg00004.html by
// Sebastien Loriot (Geometry Factory).
template<class Triangulation>
struct Traits_for_spatial_sort
:
public Triangulation::Geom_traits
{
typedef typename Triangulation::Geom_traits Gt;
typedef std::pair<const typename Triangulation::Point*, int> Point_3;
struct Less_x_3
{
bool operator()(const Point_3& p, const Point_3& q) const
{
return typename Gt::Less_x_3()(*(p.first), *(q.first));
}
};
struct Less_y_3
{
bool operator()(const Point_3& p, const Point_3& q) const
{
return typename Gt::Less_y_3()(*(p.first), *(q.first));
}
};
struct Less_z_3
{
bool operator()(const Point_3& p, const Point_3& q) const
{
return typename Gt::Less_z_3()(*(p.first), *(q.first));
}
};
Less_x_3 less_x_3_object () const {return Less_x_3();}
Less_y_3 less_y_3_object () const {return Less_y_3();}
Less_z_3 less_z_3_object () const {return Less_z_3();}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

955
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
View File

File diff suppressed because it is too large Load Diff

239
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
View File

@ -35,6 +35,10 @@ SourceFiles
conformalVoronoiMeshI.H
conformalVoronoiMesh.C
conformalVoronoiMeshIO.C
conformalVoronoiMeshConformToSurface.C
conformalVoronoiMeshFeaturePoints.C
conformalVoronoiMeshFeaturePointSpecialisations.C
conformalVoronoiMeshCalcDualMesh.C
\*---------------------------------------------------------------------------*/
@ -44,6 +48,7 @@ SourceFiles
#define CGAL_INEXACT
#include "CGALTriangulation3Ddefs.H"
#include <CGAL/Spatial_sort_traits_adapter_3.h>
#include "uint.H"
#include "ulong.H"
#include "searchableSurfaces.H"
@ -57,6 +62,8 @@ SourceFiles
#include "plane.H"
#include "SortableList.H"
#include "meshTools.H"
#include "dynamicIndexedOctree.H"
#include "dynamicTreeDataPoint.H"
#include "indexedOctree.H"
#include "treeDataPoint.H"
#include "unitConversion.H"
@ -73,6 +80,7 @@ SourceFiles
#include "processorPolyPatch.H"
#include "zeroGradientFvPatchFields.H"
#include "globalIndex.H"
#include "pointFeatureEdgesTypes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -116,6 +124,15 @@ public:
private:
// Static data
static const scalar tolParallel;
static const scalar searchConeAngle;
static const scalar searchAngleOppositeSurface;
// Private data
//- The time registry of the application
@ -146,7 +163,7 @@ private:
//- Store the feature constraining points to be reinserted after a
// triangulation clear. Maintained with relative types and indices.
std::list<Vb> featureVertices_;
List<Vb> featureVertices_;
//- Storing the locations of all of the features to be conformed to.
// Single pointField required by the featurePointTree.
@ -155,6 +172,14 @@ private:
//- Search tree for feature point locations
mutable autoPtr<indexedOctree<treeDataPoint> > featurePointTreePtr_;
//- Search tree for edge point locations
mutable autoPtr<dynamicIndexedOctree<dynamicTreeDataPoint> >
edgeLocationTreePtr_;
//- Search tree for surface point locations
mutable autoPtr<dynamicIndexedOctree<dynamicTreeDataPoint> >
surfacePtLocationTreePtr_;
//- Store locations where the cell size and alignments will be
// pre-calculated and looked up
pointField sizeAndAlignmentLocations_;
@ -170,7 +195,7 @@ private:
//- Store the surface and feature edge conformation locations to be
// reinserted
std::list<Vb> surfaceConformationVertices_;
List<Vb> surfaceConformationVertices_;
//- Method for inserting initial points. Runtime selectable.
autoPtr<initialPointsMethod> initialPointsMethod_;
@ -235,6 +260,12 @@ private:
const Foam::point& pt
) const;
//- Return the square of the local surface point exclusion distance
inline scalar surfacePtExclusionDistanceSqr
(
const Foam::point& pt
) const;
//- Return the square of the local surface search distance
inline scalar surfaceSearchDistanceSqr(const Foam::point& pt) const;
@ -279,7 +310,7 @@ private:
// processors
void insertPoints
(
std::list<Point>& points,
List<Point>& points,
bool distribute = true
);
@ -396,6 +427,8 @@ private:
//- Determine and insert point groups at the feature points
void insertFeaturePoints();
bool edgesShareNormal(const label e1, const label e2) const;
//- Create point groups at convex feature points
void createConvexFeaturePoints
(
@ -420,12 +453,24 @@ private:
DynamicList<label>& types
);
//- Fill the pointFeatureEdgesType struct with the types of feature
// edges that are attached to the point.
List<extendedFeatureEdgeMesh::edgeStatus> calcPointFeatureEdgesTypes
(
const extendedFeatureEdgeMesh& feMesh,
const labelList& pEds,
pointFeatureEdgesTypes& pFEdgeTypes
);
//- Create feature point groups if a specialisation exists for the
// structure
bool createSpecialisedFeaturePoint
(
const extendedFeatureEdgeMesh& feMesh,
label ptI,
const labelList& pEds,
const pointFeatureEdgesTypes& pFEdgeTypes,
const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
const label ptI,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
@ -434,6 +479,11 @@ private:
//- Store the locations of all of the features to be conformed to
void constructFeaturePointLocations();
List<pointIndexHit> findSurfacePtLocationsNearFeaturePoint
(
const Foam::point& featurePoint
) const;
//- Reinsert stored feature point defining points
void reinsertFeaturePoints(bool distribute = false);
@ -443,13 +493,18 @@ private:
//- Check if a location is in exclusion range around a feature point
bool nearFeaturePt(const Foam::point& pt) const;
//- Clear the entire tesselation and insert bounding points
void reset();
//- Clear the entire tesselation
// Reinsert bounding points, feature points and recalculate
// startOfInternalPoints_
void reset(const bool distribute = false);
//- Insert far points in a large bounding box to avoid dual edges
// spanning huge distances
void insertBoundingPoints();
//- Reinsert the bounding points
void reinsertBoundingPoints();
//- Insert the initial points into the triangulation, based on the
// initialPointsMethod
void insertInitialPoints();
@ -468,10 +523,10 @@ private:
//- Store data for sizeAndAlignmentLocations_, storedSizes_ and
// storedAlignments_ and initialise the sizeAndAlignmentTreePtr_
void storeSizesAndAlignments(const std::list<Point>& storePts);
void storeSizesAndAlignments(const List<Point>& storePts);
//- Restore the sizes and alignments if required
void updateSizesAndAlignments(const std::list<Point>& storePts);
void updateSizesAndAlignments(const List<Point>& storePts);
//- Demand driven construction of octree for and alignment points
const indexedOctree<treeDataPoint>& sizeAndAlignmentTree() const;
@ -526,6 +581,14 @@ private:
const Delaunay::Finite_vertices_iterator& vit
) const;
//- Return all intersections
bool dualCellSurfaceAllIntersections
(
const Delaunay::Finite_vertices_iterator& vit,
DynamicList<pointIndexHit>& info,
DynamicList<label>& hitSurface
) const;
//- Return false if the line is entirely outside the current processor
// domain, true is either point is inside, or the processor domain
// bounadry is intersected (i.e. the points are box outside but the
@ -533,6 +596,7 @@ private:
// intersect.
bool clipLineToProc
(
const Foam::point& pt,
Foam::point& a,
Foam::point& b
) const;
@ -624,21 +688,67 @@ private:
label callCount = 0
) const;
//- Find angle between the normals of two close surface points.
scalar angleBetweenSurfacePoints(Foam::point pA, Foam::point pB) const;
//- Check if a surface point is near another.
bool nearSurfacePoint
(
pointIndexHit& pHit,
label& surfaceHit,
DynamicList<Foam::point>& existingSurfacePtLocations
) const;
//- Append a point to the surface point tree and the existing list
bool appendToSurfacePtTree
(
const Foam::point& pt,
DynamicList<Foam::point>& existingSurfacePtLocations
) const;
//- Append a point to the edge location tree and the existing list
bool appendToEdgeLocationTree
(
const Foam::point& pt,
DynamicList<Foam::point>& existingEdgeLocations
) const;
//- Check if a point is near any feature edge points.
bool pointIsNearFeatureEdgeLocation(const Foam::point& pt) const;
bool pointIsNearFeatureEdgeLocation
(
const Foam::point& pt,
pointIndexHit& info
) const;
//- Check if a point is near any surface conformation points.
bool pointIsNearSurfaceLocation(const Foam::point& pt) const;
bool pointIsNearSurfaceLocation
(
const Foam::point& pt,
pointIndexHit& info
) const;
//- Check if a location is in the exclusion range of an existing feature
//- edge conformation location
bool nearFeatureEdgeLocation
(
const Foam::point& pt,
DynamicList<Foam::point>& newEdgeLocations,
pointField& existingEdgeLocations,
autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
pointIndexHit& pHit,
DynamicList<Foam::point>& existingEdgeLocations
) const;
//- Build or rebuild the edgeLocationTree
//- Build or rebuild the edge location tree
void buildEdgeLocationTree
(
autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree,
const pointField& existingEdgeLocations
const DynamicList<Foam::point>& existingEdgeLocations
) const;
//- Build or rebuild the surface point location tree
void buildSurfacePtLocationTree
(
const DynamicList<Foam::point>& existingSurfacePtLocations
) const;
//- Build or rebuild the sizeAndAlignmentTree
@ -650,8 +760,8 @@ private:
(
const Delaunay::Finite_vertices_iterator& vit,
const Foam::point& vert,
const pointIndexHit& surfHit,
label hitSurface,
const DynamicList<pointIndexHit>& surfHit,
const DynamicList<label>& hitSurface,
scalar surfacePtReplaceDistCoeffSqr,
scalar edgeSearchDistCoeffSqr,
DynamicList<pointIndexHit>& surfaceHits,
@ -659,8 +769,8 @@ private:
DynamicList<pointIndexHit>& featureEdgeHits,
DynamicList<label>& featureEdgeFeaturesHit,
DynamicList<Foam::point>& newEdgeLocations,
pointField& existingEdgeLocations,
autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
DynamicList<Foam::point>& existingEdgeLocations,
DynamicList<Foam::point>& existingSurfacePtLocations
) const;
//- Store the surface conformation with the indices offset to be
@ -996,6 +1106,9 @@ public:
//- Write Delaunay points to .obj file
void writePoints(const fileName& fName, bool internalOnly) const;
//- Write the boundary Delaunay points to .obj file
void writeBoundaryPoints(const fileName& fName) const;
//- Write list of points to file
void writePoints
(
@ -1045,9 +1158,93 @@ public:
// actual cell size (cbrt(actual cell volume)).
void writeCellSizes(const fvMesh& mesh) const;
//- Calculate and write the cell centres.
void writeCellCentres(const fvMesh& mesh) const;
//- Find the cellSet of the boundary cells which have points that
// protrude out of the surface beyond a tolerance.
void findRemainingProtrusionSet(const fvMesh& mesh) const;
//- Function inserting points into a triangulation and setting the
// index and type data of the point in the correct order. This is
// faster than inserting points individually.
//
// Adapted from a post on the CGAL lists: 2010-01/msg00004.html by
// Sebastien Loriot (Geometry Factory).
//
// @todo problems putting declaration in the .C file. Function
// prototype is not recognised.
template<class Triangulation, class Point_iterator>
void rangeInsertWithInfo
(
Point_iterator begin,
Point_iterator end,
Triangulation& T,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
typedef std::vector
<
std::pair<const typename Triangulation::Point*, label>
> vectorPairPointIndex;
vectorPairPointIndex points;
label index = 0;
for (Point_iterator it = begin; it != end; ++it)
{
points.push_back
(
std::make_pair(&(toPoint(*it)), index++)
);
}
std::random_shuffle(points.begin(), points.end());
spatial_sort
(
points.begin(),
points.end(),
Traits_for_spatial_sort<Triangulation>()
);
typename Triangulation::Cell_handle hint;
for
(
typename vectorPairPointIndex::const_iterator
p = points.begin();
p != points.end();
++p
)
{
typename Triangulation::Locate_type lt;
typename Triangulation::Cell_handle c;
label li, lj;
c = T.locate(*(p->first), lt, li, lj, hint);
typename Triangulation::Vertex_handle v
= T.insert(*(p->first), lt, c, li, lj);
label oldIndex = p->second;
label type = types[oldIndex];
if (type > Vb::vtFar)
{
// This is a member of a point pair, don't use the type
// directly (note that this routine never gets called
// for referredPoints so type will never be -procI)
type += T.number_of_vertices() - 1;
}
v->index() = indices[oldIndex] + T.number_of_vertices() - 1;
v->type() = type;
}
}
};
@ -1059,6 +1256,10 @@ public:
#include "conformalVoronoiMeshI.H"
//#ifdef NoRepository
//# include "conformalVoronoiMeshTemplates.C"
//#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

106
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
View File

@ -79,53 +79,59 @@ void Foam::conformalVoronoiMesh::calcDualMesh
}
}
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
{
std::list<Cell_handle> cells;
incident_cells(vit, std::back_inserter(cells));
bool hasProcPt = false;
for
(
std::list<Cell_handle>::iterator cit=cells.begin();
cit != cells.end();
++cit
)
{
if
(
!(*cit)->vertex(0)->real()
|| !(*cit)->vertex(1)->real()
|| !(*cit)->vertex(2)->real()
|| !(*cit)->vertex(3)->real()
)
{
hasProcPt = true;
break;
}
}
if (hasProcPt)
{
for
(
std::list<Cell_handle>::iterator cit=cells.begin();
cit != cells.end();
++cit
)
{
(*cit)->filterCount() =
cvMeshControls().filterCountSkipThreshold() + 1;
}
}
}
// REMOVED BECAUSE THIS CODE STOPS ALL FACES NEAR ANY BOUNDARY (PROC OR REAL)
// FROM BEING FILTERED.
//
// for
// (
// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
// vit != finite_vertices_end();
// vit++
// )
// {
// std::list<Cell_handle> cells;
// incident_cells(vit, std::back_inserter(cells));
//
// bool hasProcPt = false;
//
// for
// (
// std::list<Cell_handle>::iterator cit = cells.begin();
// cit != cells.end();
// ++cit
// )
// {
// // Allow filtering if any vertices are far points. Otherwise faces
// // with boundary points attached to a cell with a far point will not
// // be filtered.
// if
// (
// ( !(*cit)->vertex(0)->real() && !(*cit)->vertex(0)->farPoint() )
// || ( !(*cit)->vertex(1)->real() && !(*cit)->vertex(1)->farPoint() )
// || ( !(*cit)->vertex(2)->real() && !(*cit)->vertex(2)->farPoint() )
// || ( !(*cit)->vertex(3)->real() && !(*cit)->vertex(3)->farPoint() )
// )
// {
// hasProcPt = true;
//
// break;
// }
// }
//
// if (hasProcPt)
// {
// for
// (
// std::list<Cell_handle>::iterator cit = cells.begin();
// cit != cells.end();
// ++cit
// )
// {
// (*cit)->filterCount() =
// cvMeshControls().filterCountSkipThreshold() + 1;
// }
// }
// }
PackedBoolList boundaryPts(number_of_cells(), false);
@ -250,7 +256,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
timeCheck("End of filtering iteration");
} while (nBadQualityFaces > nInitialBadQualityFaces);
} while (nBadQualityFaces > 0); //nInitialBadQualityFaces);
}
}
@ -987,7 +993,6 @@ Foam::label Foam::conformalVoronoiMesh::collapseFaces
if (dualFace.size() < 3)
{
// This face has been collapsed already
continue;
}
@ -995,8 +1000,6 @@ Foam::label Foam::conformalVoronoiMesh::collapseFaces
if (maxFC > cvMeshControls().filterCountSkipThreshold())
{
// A vertex on this face has been limited too many
// times, skip
continue;
}
@ -1900,7 +1903,6 @@ void Foam::conformalVoronoiMesh::indexDualVertices
)
{
// This is a boundary dual vertex
boundaryPts[dualVertI] = true;
}

894
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
View File

File diff suppressed because it is too large Load Diff

642
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
View File

@ -24,33 +24,20 @@ License
\*---------------------------------------------------------------------------*/
#include "conformalVoronoiMesh.H"
#include "vectorTools.H"
using namespace Foam::vectorTools;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
Foam::List<Foam::extendedFeatureEdgeMesh::edgeStatus>
Foam::conformalVoronoiMesh::calcPointFeatureEdgesTypes
(
const extendedFeatureEdgeMesh& feMesh,
label ptI,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
const labelList& pEds,
pointFeatureEdgesTypes& pFEdgeTypes
)
{
const labelList& pEds(feMesh.pointEdges()[ptI]);
if (pEds.size() != 3)
{
// Only three edge specialisations available
return false;
}
label nExternal = 0;
label nInternal = 0;
label nFlat = 0;
label nOpen = 0;
label nMultiple = 0;
List<extendedFeatureEdgeMesh::edgeStatus> allEdStat(pEds.size());
forAll(pEds, i)
@ -65,293 +52,388 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
{
case extendedFeatureEdgeMesh::EXTERNAL:
{
nExternal++;
pFEdgeTypes.nExternal++;
break;
}
case extendedFeatureEdgeMesh::INTERNAL:
{
nInternal++;
pFEdgeTypes.nInternal++;
break;
}
case extendedFeatureEdgeMesh::FLAT:
{
nFlat++;
pFEdgeTypes.nFlat++;
break;
}
case extendedFeatureEdgeMesh::OPEN:
{
nOpen++;
pFEdgeTypes.nOpen++;
break;
}
case extendedFeatureEdgeMesh::MULTIPLE:
{
nMultiple++;
pFEdgeTypes.nMultiple++;
break;
}
case extendedFeatureEdgeMesh::NONE:
{
pFEdgeTypes.nNonFeature++;
break;
}
}
}
if (nExternal == 2 && nInternal == 1)
if (debug)
{
// if (!positionOnThisProc(pt))
// {
// continue;
// }
// Info<< "nExternal == 2 && nInternal == 1" << endl;
// const Foam::point& featPt = feMesh.points()[ptI];
// scalar ppDist = pointPairDistance(featPt);
// const vectorField& normals = feMesh.normals();
// const labelListList& edgeNormals = feMesh.edgeNormals();
// label concaveEdgeI = pEds
// [
// findIndex(allEdStat, extendedFeatureEdgeMesh::INTERNAL)
// ];
// // // Find which planes are joined to the concave edge
// // List<label> concaveEdgePlanes(2,label(-1));
// // label concaveEdgeI = concaveEdges[0];
// // // Pick up the two faces adjacent to the concave feature edge
// // const labelList& eFaces = qSurf_.edgeFaces()[concaveEdgeI];
// // label faceA = eFaces[0];
// // vector nA = qSurf_.faceNormals()[faceA];
// // scalar maxNormalDotProduct = -SMALL;
// // forAll(uniquePlaneNormals, uPN)
// // {
// // scalar normalDotProduct = nA & uniquePlaneNormals[uPN];
// // if (normalDotProduct > maxNormalDotProduct)
// // {
// // maxNormalDotProduct = normalDotProduct;
// // concaveEdgePlanes[0] = uPN;
// // }
// // }
// // label faceB = eFaces[1];
// // vector nB = qSurf_.faceNormals()[faceB];
// // maxNormalDotProduct = -SMALL;
// // forAll(uniquePlaneNormals, uPN)
// // {
// // scalar normalDotProduct = nB & uniquePlaneNormals[uPN];
// // if (normalDotProduct > maxNormalDotProduct)
// // {
// // maxNormalDotProduct = normalDotProduct;
// // concaveEdgePlanes[1] = uPN;
// // }
// // }
// // const vector& concaveEdgePlaneANormal =
// // uniquePlaneNormals[concaveEdgePlanes[0]];
// // const vector& concaveEdgePlaneBNormal =
// // uniquePlaneNormals[concaveEdgePlanes[1]];
// // label convexEdgesPlaneI;
// // if (findIndex(concaveEdgePlanes, 0) == -1)
// // {
// // convexEdgesPlaneI = 0;
// // }
// // else if (findIndex(concaveEdgePlanes, 1) == -1)
// // {
// // convexEdgesPlaneI = 1;
// // }
// // else
// // {
// // convexEdgesPlaneI = 2;
// // }
// // const vector& convexEdgesPlaneNormal =
// // uniquePlaneNormals[convexEdgesPlaneI];
// // const edge& concaveEdge = edges[concaveEdgeI];
// // // Check direction of edge, if the feature point is at the end()
// // // the reverse direction.
// // scalar edgeDirection = 1.0;
// // if (ptI == concaveEdge.end())
// // {
// // edgeDirection *= -1.0;
// // }
// const vector& concaveEdgePlaneANormal =
// normals[edgeNormals[concaveEdgeI][0]];
// const vector& concaveEdgePlaneBNormal =
// normals[edgeNormals[concaveEdgeI][1]];
// // Intersect planes parallel to the concave edge planes offset
// // by ppDist and the plane defined by featPt and the edge
// // vector.
// plane planeA
// (
// featPt + ppDist*concaveEdgePlaneANormal,
// concaveEdgePlaneANormal
// );
// plane planeB
// (
// featPt + ppDist*concaveEdgePlaneBNormal,
// concaveEdgePlaneBNormal
// );
// const vector& concaveEdgeDir = feMesh.edgeDirection
// (
// concaveEdgeI,
// ptI
// );
// Foam::point concaveEdgeLocalFeatPt = featPt + ppDist*concaveEdgeDir;
// // Finding the nearest point on the intersecting line to the
// // edge point. Floating point errors often encountered using
// // planePlaneIntersect
// plane planeF(concaveEdgeLocalFeatPt, concaveEdgeDir);
// Foam::point concaveEdgeExternalPt = planeF.planePlaneIntersect
// (
// planeA,
// planeB
// );
// label concaveEdgeExternalPtI = number_of_vertices() + 4;
// // Redefine planes to be on the feature surfaces to project
// // through
// planeA = plane(featPt, concaveEdgePlaneANormal);
// planeB = plane(featPt, concaveEdgePlaneBNormal);
// Foam::point internalPtA =
// concaveEdgeExternalPt
// - 2*planeA.distance(concaveEdgeExternalPt)
// *concaveEdgePlaneANormal;
// label internalPtAI = insertPoint
// (
// internalPtA,
// concaveEdgeExternalPtI
// );
// Foam::point internalPtB =
// concaveEdgeExternalPt
// - 2*planeB.distance(concaveEdgeExternalPt)
// *concaveEdgePlaneBNormal;
// label internalPtBI = insertPoint
// (
// internalPtB,
// concaveEdgeExternalPtI
// );
// // TEMPORARY VARIABLE TO TEST
// vector convexEdgesPlaneNormal = -concaveEdgeDir;
// plane planeC(featPt, convexEdgesPlaneNormal);
// Foam::point externalPtD =
// internalPtA
// + 2*planeC.distance(internalPtA)*convexEdgesPlaneNormal;
// insertPoint(externalPtD, internalPtAI);
// Foam::point externalPtE =
// internalPtB
// + 2*planeC.distance(internalPtB)*convexEdgesPlaneNormal;
// insertPoint(externalPtE, internalPtBI);
// insertPoint(concaveEdgeExternalPt, internalPtAI);
// Info<< nl << "# featPt " << endl;
// meshTools::writeOBJ(Info, featPt);
// Info<< "# internalPtA" << endl;
// meshTools::writeOBJ(Info, internalPtA);
// Info<< "# internalPtB" << endl;
// meshTools::writeOBJ(Info, internalPtB);
// Info<< "# externalPtD" << endl;
// meshTools::writeOBJ(Info, externalPtD);
// Info<< "# externalPtE" << endl;
// meshTools::writeOBJ(Info, externalPtE);
// scalar totalAngle = radToDeg
// (
// constant::mathematical::pi +
// acos(mag(concaveEdgePlaneANormal & concaveEdgePlaneBNormal))
// );
// if (totalAngle > cvMeshControls().maxQuadAngle())
// {
// // Add additional mitering points
// scalar angleSign = 1.0;
// if
// (
// geometryToConformTo_.outside
// (
// featPt - convexEdgesPlaneNormal*ppDist
// )
// )
// {
// angleSign = -1.0;
// }
// scalar phi =
// angleSign*acos(concaveEdgeDir & -convexEdgesPlaneNormal);
// scalar guard =
// (
// 1 + sin(phi)*ppDist/mag
// (
// concaveEdgeLocalFeatPt - concaveEdgeExternalPt
// )
// )/cos(phi) - 1;
// Foam::point internalPtF =
// concaveEdgeExternalPt
// + (2 + guard)*(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
// label internalPtFI =
// insertPoint(internalPtF, number_of_vertices() + 1);
// Foam::point externalPtG =
// internalPtF
// + 2*planeC.distance(internalPtF) * convexEdgesPlaneNormal;
// insertPoint(externalPtG, internalPtFI);
// }
// return true;
return false;
Info<< pFEdgeTypes << endl;
}
else if (nExternal == 1 && nInternal == 2)
return allEdStat;
}
bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
(
const extendedFeatureEdgeMesh& feMesh,
const labelList& pEds,
const pointFeatureEdgesTypes& pFEdgesTypes,
const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
const label ptI,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
if
(
pFEdgesTypes.nExternal == 2
&& pFEdgesTypes.nInternal == 1
&& pEds.size() == 3
)
{
Info<< "nExternal == 2 && nInternal == 1" << endl;
const Foam::point& featPt = feMesh.points()[ptI];
if (!positionOnThisProc(featPt))
{
return false;
}
const label initialNumOfPoints = pts.size();
const scalar ppDist = pointPairDistance(featPt);
const vectorField& normals = feMesh.normals();
const labelListList& edgeNormals = feMesh.edgeNormals();
label concaveEdgeI = -1;
labelList convexEdgesI(2, -1);
label nConvex = 0;
forAll(pEds, i)
{
const extendedFeatureEdgeMesh::edgeStatus& eS = allEdStat[i];
if (eS == extendedFeatureEdgeMesh::INTERNAL)
{
concaveEdgeI = pEds[i];
}
else if (eS == extendedFeatureEdgeMesh::EXTERNAL)
{
convexEdgesI[nConvex++] = pEds[i];
}
else if (eS == extendedFeatureEdgeMesh::FLAT)
{
WarningIn("Foam::conformalVoronoiMesh::"
"createSpecialisedFeaturePoint")
<< "Edge " << eS << " is flat"
<< endl;
}
else
{
FatalErrorIn("Foam::conformalVoronoiMesh::"
"createSpecialisedFeaturePoint")
<< "Edge " << eS << " not concave/convex"
<< exit(FatalError);
}
}
const vector& concaveEdgePlaneANormal =
normals[edgeNormals[concaveEdgeI][0]];
const vector& concaveEdgePlaneBNormal =
normals[edgeNormals[concaveEdgeI][1]];
// Intersect planes parallel to the concave edge planes offset
// by ppDist and the plane defined by featPt and the edge vector.
plane planeA
(
featPt + ppDist*concaveEdgePlaneANormal,
concaveEdgePlaneANormal
);
plane planeB
(
featPt + ppDist*concaveEdgePlaneBNormal,
concaveEdgePlaneBNormal
);
const vector& concaveEdgeDir = feMesh.edgeDirection
(
concaveEdgeI,
ptI
);
// Todo,needed later but want to get rid of this.
const Foam::point concaveEdgeLocalFeatPt = featPt + ppDist*concaveEdgeDir;
// Finding the nearest point on the intersecting line to the edge
// point. Floating point errors often occur using planePlaneIntersect
plane planeF(concaveEdgeLocalFeatPt, concaveEdgeDir);
const Foam::point concaveEdgeExternalPt = planeF.planePlaneIntersect
(
planeA,
planeB
);
// Redefine planes to be on the feature surfaces to project through
planeA = plane(featPt, concaveEdgePlaneANormal);
planeB = plane(featPt, concaveEdgePlaneBNormal);
const Foam::point internalPtA =
concaveEdgeExternalPt
- 2.0*planeA.distance(concaveEdgeExternalPt)
*concaveEdgePlaneANormal;
pts.append(internalPtA);
indices.append(0);
types.append(4);
const Foam::point internalPtB =
concaveEdgeExternalPt
- 2.0*planeB.distance(concaveEdgeExternalPt)
*concaveEdgePlaneBNormal;
pts.append(internalPtB);
indices.append(0);
types.append(3);
// Add the external points
Foam::point externalPtD;
Foam::point externalPtE;
vector convexEdgePlaneCNormal(0,0,0);
vector convexEdgePlaneDNormal(0,0,0);
const labelList& concaveEdgeNormals = edgeNormals[concaveEdgeI];
const labelList& convexEdgeANormals = edgeNormals[convexEdgesI[0]];
const labelList& convexEdgeBNormals = edgeNormals[convexEdgesI[1]];
forAll(concaveEdgeNormals, edgeNormalI)
{
bool convexEdgeA = false;
bool convexEdgeB = false;
forAll(convexEdgeANormals, edgeAnormalI)
{
const vector& concaveNormal
= normals[concaveEdgeNormals[edgeNormalI]];
const vector& convexNormal
= normals[convexEdgeANormals[edgeAnormalI]];
Info<< "Angle between vectors = "
<< degAngleBetween(concaveNormal, convexNormal) << endl;
// Need a looser tolerance, because sometimes adjacent triangles
// on the same surface will be slightly out of alignment.
if (areParallel(concaveNormal, convexNormal, tolParallel))
{
convexEdgeA = true;
}
}
forAll(convexEdgeBNormals, edgeBnormalI)
{
const vector& concaveNormal
= normals[concaveEdgeNormals[edgeNormalI]];
const vector& convexNormal
= normals[convexEdgeBNormals[edgeBnormalI]];
Info<< "Angle between vectors = "
<< degAngleBetween(concaveNormal, convexNormal) << endl;
// Need a looser tolerance, because sometimes adjacent triangles
// on the same surface will be slightly out of alignment.
if (areParallel(concaveNormal, convexNormal, tolParallel))
{
convexEdgeB = true;
}
}
if ((convexEdgeA && convexEdgeB) || (!convexEdgeA && !convexEdgeB))
{
WarningIn
(
"Foam::conformalVoronoiMesh"
"::createSpecialisedFeaturePoint"
)
<< "Both or neither of the convex edges share the concave "
<< "edge's normal."
<< " convexEdgeA = " << convexEdgeA
<< " convexEdgeB = " << convexEdgeB
<< endl;
// Remove points that have just been added before returning
for (label i = 0; i < 2; ++i)
{
pts.remove();
indices.remove();
types.remove();
}
return false;
}
if (convexEdgeA)
{
forAll(convexEdgeANormals, edgeAnormalI)
{
const vector& concaveNormal
= normals[concaveEdgeNormals[edgeNormalI]];
const vector& convexNormal
= normals[convexEdgeANormals[edgeAnormalI]];
if
(
!areParallel(concaveNormal, convexNormal, tolParallel)
)
{
convexEdgePlaneCNormal = convexNormal;
plane planeC(featPt, convexEdgePlaneCNormal);
externalPtD =
internalPtA
+ 2.0*planeC.distance(internalPtA)
*convexEdgePlaneCNormal;
pts.append(externalPtD);
indices.append(0);
types.append(-2);
}
}
}
if (convexEdgeB)
{
forAll(convexEdgeBNormals, edgeBnormalI)
{
const vector& concaveNormal
= normals[concaveEdgeNormals[edgeNormalI]];
const vector& convexNormal
= normals[convexEdgeBNormals[edgeBnormalI]];
if
(
!areParallel(concaveNormal, convexNormal, tolParallel)
)
{
convexEdgePlaneDNormal = convexNormal;
plane planeD(featPt, convexEdgePlaneDNormal);
externalPtE =
internalPtB
+ 2.0*planeD.distance(internalPtB)
*convexEdgePlaneDNormal;
pts.append(externalPtE);
indices.append(0);
types.append(-2);
}
}
}
}
pts.append(concaveEdgeExternalPt);
indices.append(0);
types.append(-4);
const scalar totalAngle = radToDeg
(
constant::mathematical::pi
+ radAngleBetween(concaveEdgePlaneANormal, concaveEdgePlaneBNormal)
);
if (totalAngle > cvMeshControls().maxQuadAngle())
{
// Add additional mitering points
//scalar angleSign = 1.0;
Info<<convexEdgePlaneCNormal << " " <<convexEdgePlaneDNormal << endl;
vector convexEdgesPlaneNormal = 0.5*(convexEdgePlaneCNormal + convexEdgePlaneDNormal);
plane planeM(featPt, convexEdgesPlaneNormal);
// if
// (
// geometryToConformTo_.outside
// (
// featPt - convexEdgesPlaneNormal*ppDist
// )
// )
// {
// angleSign = -1.0;
// }
// scalar phi =
// angleSign*acos(concaveEdgeDir & -convexEdgesPlaneNormal);
//
// scalar guard =
// (
// 1.0 + sin(phi)*ppDist/mag
// (
// concaveEdgeLocalFeatPt - concaveEdgeExternalPt
// )
// )/cos(phi) - 1.0;
const Foam::point internalPtF =
concaveEdgeExternalPt
//+ (2.0 + guard)*(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
+ 2.0*(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
pts.append(internalPtF);
indices.append(0);
types.append(1);
const Foam::point externalPtG =
internalPtF
+ 2.0*planeM.distance(internalPtF)*convexEdgesPlaneNormal;
pts.append(externalPtG);
indices.append(0);
types.append(-1);
}
if (debug)
{
for (label ptI = initialNumOfPoints; ptI < pts.size(); ++ptI)
{
Info<< "Point " << ptI << " : ";
meshTools::writeOBJ(Info, pts[ptI]);
}
}
return true;
}
else if (pFEdgesTypes.nExternal == 1 && pFEdgesTypes.nInternal == 2)
{
// Info<< "nExternal == 1 && nInternal == 2" << endl;

923
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C Normal file
View File

@ -0,0 +1,923 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
\*---------------------------------------------------------------------------*/
#include "conformalVoronoiMesh.H"
#include "vectorTools.H"
using namespace Foam::vectorTools;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void Foam::conformalVoronoiMesh::insertBoundingPoints()
{
pointField farPts = geometryToConformTo_.globalBounds().points();
// Shift corners of bounds relative to origin
farPts -= geometryToConformTo_.globalBounds().midpoint();
// Scale the box up
farPts *= 10.0;
// Shift corners of bounds back to be relative to midpoint
farPts += geometryToConformTo_.globalBounds().midpoint();
limitBounds_ = treeBoundBox(farPts);
forAll(farPts, fPI)
{
insertPoint(farPts[fPI], Vb::vtFar);
}
}
void Foam::conformalVoronoiMesh::createEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
label edgeI = edHit.index();
extendedFeatureEdgeMesh::edgeStatus edStatus = feMesh.getEdgeStatus(edgeI);
switch (edStatus)
{
case extendedFeatureEdgeMesh::EXTERNAL:
{
createExternalEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::INTERNAL:
{
createInternalEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::FLAT:
{
createFlatEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::OPEN:
{
//createOpenEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::MULTIPLE:
{
//createMultipleEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::NONE:
{
break;
}
}
}
void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const Foam::point& edgePt = edHit.hitPoint();
scalar ppDist = pointPairDistance(edgePt);
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
// As this is an external edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
if (areParallel(nA, nB))
{
// The normals are nearly parallel, so this is too sharp a feature to
// conform to.
return;
}
// Normalised distance of reference point from edge point
vector refVec((nA + nB)/(1 + (nA & nB)));
if (magSqr(refVec) > sqr(5.0))
{
// Limit the size of the conformation
ppDist *= 5.0/mag(refVec);
// Pout<< nl << "createExternalEdgePointGroup limit "
// << "edgePt " << edgePt << nl
// << "refVec " << refVec << nl
// << "mag(refVec) " << mag(refVec) << nl
// << "ppDist " << ppDist << nl
// << "nA " << nA << nl
// << "nB " << nB << nl
// << "(nA & nB) " << (nA & nB) << nl
// << endl;
}
// Convex. So refPt will be inside domain and hence a master point
Foam::point refPt = edgePt - ppDist*refVec;
// Result when the points are eventually inserted.
// Add number_of_vertices() at insertion of first vertex to all numbers:
// pt index type
// refPt 0 1
// reflectedA 1 0
// reflectedB 2 0
// Insert the master point pairing the the first slave
pts.append(refPt);
indices.append(0);
types.append(1);
// Insert the slave points by reflecting refPt in both faces.
// with each slave refering to the master
Foam::point reflectedA = refPt + 2*ppDist*nA;
pts.append(reflectedA);
indices.append(0);
types.append(-1);
Foam::point reflectedB = refPt + 2*ppDist*nB;
pts.append(reflectedB);
indices.append(0);
types.append(-2);
}
void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const Foam::point& edgePt = edHit.hitPoint();
scalar ppDist = pointPairDistance(edgePt);
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
// As this is an external edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
if (areParallel(nA, nB))
{
// The normals are nearly parallel, so this is too sharp a feature to
// conform to.
return;
}
// Normalised distance of reference point from edge point
vector refVec((nA + nB)/(1 + (nA & nB)));
if (magSqr(refVec) > sqr(5.0))
{
// Limit the size of the conformation
ppDist *= 5.0/mag(refVec);
// Pout<< nl << "createInternalEdgePointGroup limit "
// << "edgePt " << edgePt << nl
// << "refVec " << refVec << nl
// << "mag(refVec) " << mag(refVec) << nl
// << "ppDist " << ppDist << nl
// << "nA " << nA << nl
// << "nB " << nB << nl
// << "(nA & nB) " << (nA & nB) << nl
// << endl;
}
// Concave. master and reflected points inside the domain.
Foam::point refPt = edgePt - ppDist*refVec;
// Generate reflected master to be outside.
Foam::point reflMasterPt = refPt + 2*(edgePt - refPt);
// Reflect reflMasterPt in both faces.
Foam::point reflectedA = reflMasterPt - 2*ppDist*nA;
Foam::point reflectedB = reflMasterPt - 2*ppDist*nB;
scalar totalAngle =
radToDeg(constant::mathematical::pi + radAngleBetween(nA, nB));
// Number of quadrants the angle should be split into
int nQuads = int(totalAngle/cvMeshControls().maxQuadAngle()) + 1;
// The number of additional master points needed to obtain the
// required number of quadrants.
int nAddPoints = min(max(nQuads - 2, 0), 2);
// index of reflMaster
label reflectedMaster = 2 + nAddPoints;
// Add number_of_vertices() at insertion of first vertex to all numbers:
// Result for nAddPoints 1 when the points are eventually inserted
// pt index type
// reflectedA 0 2
// reflectedB 1 2
// reflMasterPt 2 0
// Result for nAddPoints 1 when the points are eventually inserted
// pt index type
// reflectedA 0 3
// reflectedB 1 3
// refPt 2 3
// reflMasterPt 3 0
// Result for nAddPoints 2 when the points are eventually inserted
// pt index type
// reflectedA 0 4
// reflectedB 1 4
// reflectedAa 2 4
// reflectedBb 3 4
// reflMasterPt 4 0
// Master A is inside.
pts.append(reflectedA);
indices.append(0);
types.append(reflectedMaster--);
// Master B is inside.
pts.append(reflectedB);
indices.append(0);
types.append(reflectedMaster--);
if (nAddPoints == 1)
{
// One additinal point is the reflection of the slave point,
// i.e. the original reference point
pts.append(refPt);
indices.append(0);
types.append(reflectedMaster--);
}
else if (nAddPoints == 2)
{
Foam::point reflectedAa = refPt + ppDist*nB;
pts.append(reflectedAa);
indices.append(0);
types.append(reflectedMaster--);
Foam::point reflectedBb = refPt + ppDist*nA;
pts.append(reflectedBb);
indices.append(0);
types.append(reflectedMaster--);
}
// Slave is outside.
pts.append(reflMasterPt);
indices.append(0);
types.append(-(nAddPoints + 2));
}
void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const Foam::point& edgePt = edHit.hitPoint();
const scalar ppDist = pointPairDistance(edgePt);
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
// As this is a flat edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
// Average normal to remove any bias to one side, although as this
// is a flat edge, the normals should be essentially the same
const vector n = 0.5*(nA + nB);
// Direction along the surface to the control point, sense of edge
// direction not important, as +s and -s can be used because this
// is a flat edge
vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
createPointPair(ppDist, edgePt + s, n, pts, indices, types);
createPointPair(ppDist, edgePt - s, n, pts, indices, types);
}
void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
Info<< "NOT INSERTING OPEN EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
}
void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
Info<< "NOT INSERTING MULTIPLE EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
}
void Foam::conformalVoronoiMesh::reinsertBoundingPoints()
{
tmp<pointField> farPts = limitBounds_.points();
forAll(farPts(), fPI)
{
insertPoint(farPts()[fPI], Vb::vtFar);
}
}
void Foam::conformalVoronoiMesh::reinsertFeaturePoints(bool distribute)
{
Info<< nl << "Reinserting stored feature points" << endl;
label preReinsertionSize(number_of_vertices());
if (distribute)
{
DynamicList<Foam::point> pointsToInsert;
DynamicList<label> indices;
DynamicList<label> types;
for
(
List<Vb>::iterator vit=featureVertices_.begin();
vit != featureVertices_.end();
++vit
)
{
pointsToInsert.append(topoint(vit->point()));
indices.append(vit->index());
types.append(vit->type());
}
// Insert distributed points
insertPoints(pointsToInsert, indices, types, true);
// Save points in new distribution
featureVertices_.clear();
featureVertices_.setSize(pointsToInsert.size());
forAll(pointsToInsert, pI)
{
featureVertices_[pI] =
Vb(toPoint(pointsToInsert[pI]), indices[pI], types[pI]);
}
}
else
{
for
(
List<Vb>::iterator vit=featureVertices_.begin();
vit != featureVertices_.end();
++vit
)
{
// Assuming that all of the reinsertions are pair points, and that
// the index and type are relative, i.e. index 0 and type relative
// to it.
insertPoint
(
vit->point(),
vit->index() + number_of_vertices(),
vit->type() + number_of_vertices()
);
}
}
Info<< " Reinserted "
<< returnReduce
(
label(number_of_vertices()) - preReinsertionSize,
sumOp<label>()
)
<< " vertices" << endl;
}
bool Foam::conformalVoronoiMesh::edgesShareNormal
(
const label e1,
const label e2
) const
{
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
const vectorField& e1normals = feMesh.edgeNormals(e1);
const vectorField& e2normals = feMesh.edgeNormals(e2);
forAll(e1normals, nI1)
{
forAll(e2normals, nI2)
{
if (degAngleBetween(e1normals[nI1], e2normals[nI2]) < 1)
{
return true;
}
}
}
}
return false;
}
void Foam::conformalVoronoiMesh::createConvexFeaturePoints
(
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
for
(
label ptI = feMesh.convexStart();
ptI < feMesh.concaveStart();
ptI++
)
{
const Foam::point& featPt = feMesh.points()[ptI];
if (!positionOnThisProc(featPt))
{
continue;
}
const vectorField& featPtNormals = feMesh.featurePointNormals(ptI);
const scalar ppDist = - pointPairDistance(featPt);
vector cornerNormal = sum(featPtNormals);
cornerNormal /= mag(cornerNormal);
Foam::point internalPt = featPt + ppDist*cornerNormal;
// Result when the points are eventually inserted (example n = 4)
// Add number_of_vertices() at insertion of first vertex to all
// numbers:
// pt index type
// internalPt 0 1
// externalPt0 1 0
// externalPt1 2 0
// externalPt2 3 0
// externalPt3 4 0
pts.append(internalPt);
indices.append(0);
types.append(1);
label internalPtIndex = -1;
forAll(featPtNormals, nI)
{
const vector& n = featPtNormals[nI];
plane planeN = plane(featPt, n);
Foam::point externalPt =
internalPt + 2.0 * planeN.distance(internalPt) * n;
pts.append(externalPt);
indices.append(0);
types.append(internalPtIndex--);
}
}
}
}
void Foam::conformalVoronoiMesh::createConcaveFeaturePoints
(
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
for
(
label ptI = feMesh.concaveStart();
ptI < feMesh.mixedStart();
ptI++
)
{
const Foam::point& featPt = feMesh.points()[ptI];
if (!positionOnThisProc(featPt))
{
continue;
}
const vectorField& featPtNormals = feMesh.featurePointNormals(ptI);
const scalar ppDist = pointPairDistance(featPt);
vector cornerNormal = sum(featPtNormals);
cornerNormal /= mag(cornerNormal);
Foam::point externalPt = featPt + ppDist*cornerNormal;
label externalPtIndex = featPtNormals.size();
// Result when the points are eventually inserted (example n = 5)
// Add number_of_vertices() at insertion of first vertex to all
// numbers:
// pt index type
// internalPt0 0 5
// internalPt1 1 5
// internalPt2 2 5
// internalPt3 3 5
// internalPt4 4 5
// externalPt 5 4
forAll(featPtNormals, nI)
{
const vector& n = featPtNormals[nI];
const plane planeN = plane(featPt, n);
const Foam::point internalPt =
externalPt - 2.0 * planeN.distance(externalPt) * n;
pts.append(internalPt);
indices.append(0);
types.append(externalPtIndex--);
}
pts.append(externalPt);
indices.append(0);
types.append(-1);
}
}
}
void Foam::conformalVoronoiMesh::createMixedFeaturePoints
(
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
if (cvMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
{
Info<< nl << "Skipping specialised handling for mixed feature points"
<< endl;
return;
}
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh = feMeshes[i];
const labelListList& pointsEdges = feMesh.pointEdges();
const pointField& points = feMesh.points();
for
(
label ptI = feMesh.mixedStart();
ptI < feMesh.nonFeatureStart();
ptI++
)
{
const labelList& pEds = pointsEdges[ptI];
pointFeatureEdgesTypes pFEdgeTypes(ptI);
const List<extendedFeatureEdgeMesh::edgeStatus> allEdStat
= calcPointFeatureEdgesTypes(feMesh, pEds, pFEdgeTypes);
bool specialisedSuccess = false;
if (cvMeshControls().specialiseFeaturePoints())
{
specialisedSuccess = createSpecialisedFeaturePoint
(
feMesh, pEds, pFEdgeTypes, allEdStat, ptI,
pts, indices, types
);
}
if (!specialisedSuccess)
{
// Specialisations available for some mixed feature points. For
// non-specialised feature points, inserting mixed internal and
// external edge groups at feature point.
// Skipping unsupported mixed feature point types
bool skipEdge = false;
forAll(pEds, e)
{
const label edgeI = pEds[e];
const extendedFeatureEdgeMesh::edgeStatus edStatus
= feMesh.getEdgeStatus(edgeI);
if
(
edStatus == extendedFeatureEdgeMesh::OPEN
|| edStatus == extendedFeatureEdgeMesh::MULTIPLE
)
{
Info<< "Edge type " << edStatus
<< " found for mixed feature point " << ptI
<< ". Not supported."
<< endl;
skipEdge = true;
}
}
if (skipEdge)
{
Info<< "Skipping point " << ptI << nl << endl;
continue;
}
const Foam::point& pt = points[ptI];
if (!positionOnThisProc(pt))
{
continue;
}
const scalar edgeGroupDistance = mixedFeaturePointDistance(pt);
forAll(pEds, e)
{
const label edgeI = pEds[e];
const Foam::point edgePt =
pt + edgeGroupDistance*feMesh.edgeDirection(edgeI, ptI);
const pointIndexHit edgeHit(true, edgePt, edgeI);
createEdgePointGroup(feMesh, edgeHit, pts, indices, types);
}
}
}
}
}
//
//void Foam::conformalVoronoiMesh::createFeaturePoints
//(
// DynamicList<Foam::point>& pts,
// DynamicList<label>& indices,
// DynamicList<label>& types
//)
//{
// const PtrList<extendedFeatureEdgeMesh>& feMeshes
// (
// geometryToConformTo_.features()
// );
//
// forAll(feMeshes, i)
// {
// const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
//
// for
// (
// label ptI = feMesh.convexStart();
// ptI < feMesh.mixedStart();
// ++ptI
// )
// {
// const Foam::point& featPt = feMesh.points()[ptI];
//
// if (!positionOnThisProc(featPt))
// {
// continue;
// }
//
// const scalar searchRadiusSqr = 5.0*targetCellSize(featPt);
//
// labelList indices =
// surfacePtLocationTreePtr_().findSphere(featPt, searchRadiusSqr);
//
// pointField nearestSurfacePoints(indices.size());
//
// forAll(indices, pI)
// {
// nearestSurfacePoints[pI] =
// surfaceConformationVertices_[indices[pI]];
// }
//
// forAll()
//
// // Now find the nearest points within the edge cones.
//
// // Calculate preliminary surface point locations
//
//
// }
// }
//}
void Foam::conformalVoronoiMesh::insertFeaturePoints()
{
Info<< nl << "Conforming to feature points" << endl;
DynamicList<Foam::point> pts;
DynamicList<label> indices;
DynamicList<label> types;
const label preFeaturePointSize = number_of_vertices();
createConvexFeaturePoints(pts, indices, types);
createConcaveFeaturePoints(pts, indices, types);
// createFeaturePoints(pts, indices, types);
createMixedFeaturePoints(pts, indices, types);
// Insert the created points, distributing to the appropriate processor
insertPoints(pts, indices, types, true);
if (cvMeshControls().objOutput())
{
writePoints("featureVertices.obj", pts);
}
label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
if (Pstream::parRun())
{
reduce(nFeatureVertices, sumOp<label>());
}
if (nFeatureVertices > 0)
{
Info<< " Inserted " << nFeatureVertices
<< " feature vertices" << endl;
}
featureVertices_.clear();
featureVertices_.setSize(pts.size());
forAll(pts, pI)
{
featureVertices_[pI] = Vb(toPoint(pts[pI]), indices[pI], types[pI]);
}
constructFeaturePointLocations();
}
void Foam::conformalVoronoiMesh::constructFeaturePointLocations()
{
DynamicList<Foam::point> ftPtLocs;
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
if (cvMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
{
// Ignoring mixed feature points
for
(
label ptI = feMesh.convexStart();
ptI < feMesh.mixedStart();
ptI++
)
{
ftPtLocs.append(feMesh.points()[ptI]);
}
}
else
{
for
(
label ptI = feMesh.convexStart();
ptI < feMesh.nonFeatureStart();
ptI++
)
{
ftPtLocs.append(feMesh.points()[ptI]);
}
}
}
featurePointLocations_.transfer(ftPtLocs);
}
Foam::List<Foam::pointIndexHit>
Foam::conformalVoronoiMesh::findSurfacePtLocationsNearFeaturePoint
(
const Foam::point& featurePoint
) const
{
DynamicList<pointIndexHit> dynPointList;
const scalar searchRadiusSqr = 3*targetCellSize(featurePoint);
labelList surfacePtList = surfacePtLocationTreePtr_().findSphere
(
featurePoint,
searchRadiusSqr
);
forAll(surfacePtList, elemI)
{
label index = surfacePtList[elemI];
const Foam::point& p
= surfacePtLocationTreePtr_().shapes().shapePoints()[index];
pointIndexHit nearHit(true, p, index);
dynPointList.append(nearHit);
}
return dynPointList.shrink();
}

14
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
View File

@ -186,6 +186,20 @@ inline Foam::scalar Foam::conformalVoronoiMesh::featureEdgeExclusionDistanceSqr
}
inline Foam::scalar Foam::conformalVoronoiMesh::surfacePtExclusionDistanceSqr
(
const Foam::point& pt
) const
{
return
sqr
(
targetCellSize(pt)
*cvMeshControls().surfacePtExclusionDistanceCoeff()
);
}
inline Foam::scalar Foam::conformalVoronoiMesh::surfaceSearchDistanceSqr
(
const Foam::point& pt

75
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
View File

@ -175,6 +175,30 @@ void Foam::conformalVoronoiMesh::writePoints
}
void Foam::conformalVoronoiMesh::writeBoundaryPoints
(
const fileName& fName
) const
{
OFstream str(runTime_.path()/fName);
Pout<< nl << "Writing boundary points to " << str.name() << endl;
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
++vit
)
{
if (!vit->internalPoint())
{
meshTools::writeOBJ(str, topoint(vit->point()));
}
}
}
void Foam::conformalVoronoiMesh::writePoints
(
const fileName& fName,
@ -597,29 +621,10 @@ void Foam::conformalVoronoiMesh::writeMesh
<< exit(FatalError);
}
// pointIOField cellCs
// (
// IOobject
// (
// "cellCentres",
// mesh.pointsInstance(),
// polyMesh::meshSubDir,
// mesh,
// IOobject::NO_READ,
// IOobject::AUTO_WRITE
// ),
// cellCentres
// );
// Info<< nl
// << "Writing " << cellCs.name()
// << " to " << cellCs.instance()
// << endl;
// cellCs.write();
writeCellSizes(mesh);
writeCellCentres(mesh);
findRemainingProtrusionSet(mesh);
}
@ -805,6 +810,34 @@ void Foam::conformalVoronoiMesh::writeCellSizes
}
void Foam::conformalVoronoiMesh::writeCellCentres
(
const fvMesh& mesh
) const
{
Info<< "Writing components of cellCentre positions to volScalarFields"
<< " ccx, ccy, ccz in " << runTime_.timeName() << endl;
for (direction i=0; i<vector::nComponents; i++)
{
volScalarField cci
(
IOobject
(
"cc" + word(vector::componentNames[i]),
runTime_.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh.C().component(i)
);
cci.write();
}
}
void Foam::conformalVoronoiMesh::findRemainingProtrusionSet
(
const fvMesh& mesh

290
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
View File

@ -33,6 +33,10 @@ Description
An indexed form of CGAL::Triangulation_cell_base_3<K> used to keep
track of the Delaunay cells (tets) in the tessellation.
SourceFiles
indexedCellI.H
indexedCell.C
\*---------------------------------------------------------------------------*/
#ifndef indexedCell_H
@ -119,248 +123,72 @@ public:
};
indexedCell()
:
Cb(),
index_(ctFar),
filterCount_(0)
{}
// Constructors
inline indexedCell();
indexedCell
(
Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
)
:
Cb(v0, v1, v2, v3),
index_(ctFar),
filterCount_(0)
{}
indexedCell
(
Vertex_handle v0,
Vertex_handle v1,
Vertex_handle v2,
Vertex_handle v3,
Cell_handle n0,
Cell_handle n1,
Cell_handle n2,
Cell_handle n3
)
:
Cb(v0, v1, v2, v3, n0, n1, n2, n3),
index_(ctFar),
filterCount_(0)
{}
int& cellIndex()
{
return index_;
}
int cellIndex() const
{
return index_;
}
inline bool farCell() const
{
return index_ == ctFar;
}
inline int& filterCount()
{
return filterCount_;
}
inline int filterCount() const
{
return filterCount_;
}
//- Is the Delaunay cell real, i.e. any real vertex
inline bool real() const
{
return
inline indexedCell
(
this->vertex(0)->real()
|| this->vertex(1)->real()
|| this->vertex(2)->real()
|| this->vertex(3)->real()
);
}
//- Does the Dual vertex form part of a processor patch
inline bool parallelDualVertex() const
{
return
(
this->vertex(0)->referred()
|| this->vertex(1)->referred()
|| this->vertex(2)->referred()
|| this->vertex(3)->referred()
);
}
//- Does the Dual vertex form part of a processor patch
inline Foam::label dualVertexMasterProc() const
{
if (!parallelDualVertex())
{
return -1;
}
// The master processor is the lowest numbered of the four on this tet.
int masterProc = Foam::Pstream::nProcs() + 1;
for (int i = 0; i < 4; i++)
{
if (this->vertex(i)->referred())
{
masterProc = min(masterProc, this->vertex(i)->procIndex());
}
else
{
masterProc = min(masterProc, Foam::Pstream::myProcNo());
}
}
return masterProc;
}
inline Foam::FixedList<Foam::label, 4> processorsAttached() const
{
if (!parallelDualVertex())
{
return Foam::FixedList<Foam::label, 4>(Foam::Pstream::myProcNo());
}
Foam::FixedList<Foam::label, 4> procsAttached
(
Foam::Pstream::myProcNo()
Vertex_handle v0,
Vertex_handle v1,
Vertex_handle v2,
Vertex_handle v3
);
for (int i = 0; i < 4; i++)
{
if (this->vertex(i)->referred())
{
procsAttached[i] = this->vertex(i)->procIndex();
}
}
return procsAttached;
}
// Using the globalIndex object, return a list of four (sorted) global
// Delaunay vertex indices that uniquely identify this tet in parallel
inline Foam::tetCell vertexGlobalIndices
(
const Foam::globalIndex& globalDelaunayVertexIndices
) const
{
// tetVertexGlobalIndices
Foam::tetCell tVGI;
for (int i = 0; i < 4; i++)
{
Vertex_handle v = this->vertex(i);
// Finding the global index of each Delaunay vertex
if (v->referred())
{
// Referred vertices may have negative indices
tVGI[i] = globalDelaunayVertexIndices.toGlobal
(
v->procIndex(),
Foam::mag(v->index())
);
}
else
{
tVGI[i] = globalDelaunayVertexIndices.toGlobal
(
Foam::Pstream::myProcNo(),
v->index()
);
}
}
// bubble sort
for (int i = 0; i < tVGI.size(); i++)
{
for (int j = tVGI.size() - 1 ; j > i; j--)
{
if (tVGI[j - 1] > tVGI[j])
{
Foam::Swap(tVGI[j - 1], tVGI[j]);
}
}
}
return tVGI;
}
// Is the Delaunay cell part of the final dual mesh, i.e. any vertex form
// part of the internal or boundary definition
inline bool internalOrBoundaryDualVertex() const
{
return
inline indexedCell
(
this->vertex(0)->internalOrBoundaryPoint()
|| this->vertex(1)->internalOrBoundaryPoint()
|| this->vertex(2)->internalOrBoundaryPoint()
|| this->vertex(3)->internalOrBoundaryPoint()
Vertex_handle v0,
Vertex_handle v1,
Vertex_handle v2,
Vertex_handle v3,
Cell_handle n0,
Cell_handle n1,
Cell_handle n2,
Cell_handle n3
);
}
// Is the Delaunay cell real or referred (or mixed), i.e. all vertices form
// part of the real or referred internal or boundary definition
inline bool anyInternalOrBoundaryDualVertex() const
{
return
// Member Functions
inline int& cellIndex();
inline int cellIndex() const;
inline bool farCell() const;
inline int& filterCount();
inline int filterCount() const;
//- Is the Delaunay cell real, i.e. any real vertex
inline bool real() const;
//- Does the Dual vertex form part of a processor patch
inline bool parallelDualVertex() const;
//- Does the Dual vertex form part of a processor patch
inline Foam::label dualVertexMasterProc() const;
inline Foam::FixedList<Foam::label, 4> processorsAttached() const;
//- Using the globalIndex object, return a list of four (sorted) global
// Delaunay vertex indices that uniquely identify this tet in parallel
inline Foam::tetCell vertexGlobalIndices
(
this->vertex(0)->anyInternalOrBoundaryPoint()
|| this->vertex(1)->anyInternalOrBoundaryPoint()
|| this->vertex(2)->anyInternalOrBoundaryPoint()
|| this->vertex(3)->anyInternalOrBoundaryPoint()
);
}
const Foam::globalIndex& globalDelaunayVertexIndices
) const;
//- Is the Delaunay cell part of the final dual mesh, i.e. any vertex
// form part of the internal or boundary definition
inline bool internalOrBoundaryDualVertex() const;
// A dual vertex on the boundary will result from a Delaunay cell with
// least one Delaunay vertex outside and at least one inside
inline bool boundaryDualVertex() const
{
return
(
(
this->vertex(0)->internalOrBoundaryPoint()
|| this->vertex(1)->internalOrBoundaryPoint()
|| this->vertex(2)->internalOrBoundaryPoint()
|| this->vertex(3)->internalOrBoundaryPoint()
)
&& (
!this->vertex(0)->internalOrBoundaryPoint()
|| !this->vertex(1)->internalOrBoundaryPoint()
|| !this->vertex(2)->internalOrBoundaryPoint()
|| !this->vertex(3)->internalOrBoundaryPoint()
)
);
}
//- Is the Delaunay cell real or referred (or mixed), i.e. all vertices
// form part of the real or referred internal or boundary definition
inline bool anyInternalOrBoundaryDualVertex() const;
//- A dual vertex on the boundary will result from a Delaunay cell with
// least one Delaunay vertex outside and at least one inside
inline bool boundaryDualVertex() const;
// Info
@ -387,6 +215,8 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "indexedCellI.H"
#ifdef NoRepository
# include "indexedCell.C"
#endif

288
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H Normal file
View File

@ -0,0 +1,288 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
As a special exception, you have permission to link this program with the
CGAL library and distribute executables, as long as you follow the
requirements of the GNU GPL in regard to all of the software in the
executable aside from CGAL.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Gt, class Cb>
CGAL::indexedCell<Gt, Cb>::indexedCell()
:
Cb(),
index_(ctFar),
filterCount_(0)
{}
template<class Gt, class Cb>
CGAL::indexedCell<Gt, Cb>::indexedCell
(
Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
)
:
Cb(v0, v1, v2, v3),
index_(ctFar),
filterCount_(0)
{}
template<class Gt, class Cb>
CGAL::indexedCell<Gt, Cb>::indexedCell
(
Vertex_handle v0,
Vertex_handle v1,
Vertex_handle v2,
Vertex_handle v3,
Cell_handle n0,
Cell_handle n1,
Cell_handle n2,
Cell_handle n3
)
:
Cb(v0, v1, v2, v3, n0, n1, n2, n3),
index_(ctFar),
filterCount_(0)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Gt, class Cb>
int& CGAL::indexedCell<Gt, Cb>::cellIndex()
{
return index_;
}
template<class Gt, class Cb>
int CGAL::indexedCell<Gt, Cb>::cellIndex() const
{
return index_;
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::farCell() const
{
return index_ == ctFar;
}
template<class Gt, class Cb>
inline int& CGAL::indexedCell<Gt, Cb>::filterCount()
{
return filterCount_;
}
template<class Gt, class Cb>
inline int CGAL::indexedCell<Gt, Cb>::filterCount() const
{
return filterCount_;
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::real() const
{
return
(
this->vertex(0)->real()
|| this->vertex(1)->real()
|| this->vertex(2)->real()
|| this->vertex(3)->real()
);
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::parallelDualVertex() const
{
return
(
this->vertex(0)->referred()
|| this->vertex(1)->referred()
|| this->vertex(2)->referred()
|| this->vertex(3)->referred()
);
}
template<class Gt, class Cb>
inline Foam::label CGAL::indexedCell<Gt, Cb>::dualVertexMasterProc() const
{
if (!parallelDualVertex())
{
return -1;
}
// The master processor is the lowest numbered of the four on this tet.
int masterProc = Foam::Pstream::nProcs() + 1;
for (int i = 0; i < 4; i++)
{
if (this->vertex(i)->referred())
{
masterProc = min(masterProc, this->vertex(i)->procIndex());
}
else
{
masterProc = min(masterProc, Foam::Pstream::myProcNo());
}
}
return masterProc;
}
template<class Gt, class Cb>
inline Foam::FixedList<Foam::label, 4>
CGAL::indexedCell<Gt, Cb>::processorsAttached() const
{
if (!parallelDualVertex())
{
return Foam::FixedList<Foam::label, 4>(Foam::Pstream::myProcNo());
}
Foam::FixedList<Foam::label, 4> procsAttached
(
Foam::Pstream::myProcNo()
);
for (int i = 0; i < 4; i++)
{
if (this->vertex(i)->referred())
{
procsAttached[i] = this->vertex(i)->procIndex();
}
}
return procsAttached;
}
template<class Gt, class Cb>
inline Foam::tetCell CGAL::indexedCell<Gt, Cb>::vertexGlobalIndices
(
const Foam::globalIndex& globalDelaunayVertexIndices
) const
{
// tetVertexGlobalIndices
Foam::tetCell tVGI;
for (int i = 0; i < 4; i++)
{
Vertex_handle v = this->vertex(i);
// Finding the global index of each Delaunay vertex
if (v->referred())
{
// Referred vertices may have negative indices
tVGI[i] = globalDelaunayVertexIndices.toGlobal
(
v->procIndex(),
Foam::mag(v->index())
);
}
else
{
tVGI[i] = globalDelaunayVertexIndices.toGlobal
(
Foam::Pstream::myProcNo(),
v->index()
);
}
}
// bubble sort
for (int i = 0; i < tVGI.size(); i++)
{
for (int j = tVGI.size() - 1 ; j > i; j--)
{
if (tVGI[j - 1] > tVGI[j])
{
Foam::Swap(tVGI[j - 1], tVGI[j]);
}
}
}
return tVGI;
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::internalOrBoundaryDualVertex() const
{
return
(
this->vertex(0)->internalOrBoundaryPoint()
|| this->vertex(1)->internalOrBoundaryPoint()
|| this->vertex(2)->internalOrBoundaryPoint()
|| this->vertex(3)->internalOrBoundaryPoint()
);
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::anyInternalOrBoundaryDualVertex() const
{
return
(
this->vertex(0)->anyInternalOrBoundaryPoint()
|| this->vertex(1)->anyInternalOrBoundaryPoint()
|| this->vertex(2)->anyInternalOrBoundaryPoint()
|| this->vertex(3)->anyInternalOrBoundaryPoint()
);
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::boundaryDualVertex() const
{
return
(
(
this->vertex(0)->internalOrBoundaryPoint()
|| this->vertex(1)->internalOrBoundaryPoint()
|| this->vertex(2)->internalOrBoundaryPoint()
|| this->vertex(3)->internalOrBoundaryPoint()
)
&& (
!this->vertex(0)->internalOrBoundaryPoint()
|| !this->vertex(1)->internalOrBoundaryPoint()
|| !this->vertex(2)->internalOrBoundaryPoint()
|| !this->vertex(3)->internalOrBoundaryPoint()
)
);
}
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //

304
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
View File

@ -33,6 +33,10 @@ Description
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep
track of the Delaunay vertices in the tessellation.
SourceFiles
indexedVertexI.H
indexedVertex.C
\*---------------------------------------------------------------------------*/
#ifndef indexedVertex_H
@ -88,12 +92,15 @@ class indexedVertex
// Lowest numbered is inside one (master).
int type_;
// Required alignment of the dual cell of this vertex
//- Required alignment of the dual cell of this vertex
Foam::tensor alignment_;
// Target size of the dual cell of this vertex
//- Target size of the dual cell of this vertex
Foam::scalar targetCellSize_;
//- Specify whether the vertex is fixed or movable.
bool vertexFixed_;
public:
@ -104,6 +111,12 @@ public:
vtFar = INT_MIN + 2
};
enum vertexMotion
{
fixed = 0,
movable = 1
};
typedef typename Vb::Vertex_handle Vertex_handle;
typedef typename Vb::Cell_handle Cell_handle;
typedef typename Vb::Point Point;
@ -116,268 +129,105 @@ public:
};
indexedVertex()
:
Vb(),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
// Constructors
inline indexedVertex();
indexedVertex(const Point& p)
:
Vb(p),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
inline indexedVertex(const Point& p);
inline indexedVertex(const Point& p, int index, int type);
indexedVertex(const Point& p, int index, int type)
:
Vb(p),
index_(index),
type_(type),
alignment_(),
targetCellSize_(0.0)
{}
inline indexedVertex(const Point& p, Cell_handle f);
inline indexedVertex(Cell_handle f);
indexedVertex(const Point& p, Cell_handle f)
:
Vb(f, p),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
// Member Functions
indexedVertex(Cell_handle f)
:
Vb(f),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
inline int& index();
inline int index() const;
int& index()
{
return index_;
}
inline int& type();
inline int type() const;
int index() const
{
return index_;
}
inline Foam::tensor& alignment();
inline const Foam::tensor& alignment() const;
int& type()
{
return type_;
}
inline Foam::scalar& targetCellSize();
inline Foam::scalar targetCellSize() const;
int type() const
{
return type_;
}
inline bool uninitialised() const;
//- Is point a far-point
inline bool farPoint() const;
inline Foam::tensor& alignment()
{
return alignment_;
}
//- Is point internal, i.e. not on boundary
inline bool internalPoint() const;
//- Is point internal, i.e. not on boundary, external query.
inline static bool internalPoint(int type);
inline const Foam::tensor& alignment() const
{
return alignment_;
}
// is this a referred vertex
inline bool referred() const;
// is this a referred internal or boundary vertex
inline bool referredInternalOrBoundaryPoint() const;
inline Foam::scalar& targetCellSize()
{
return targetCellSize_;
}
// is this a referred external (pair slave) vertex
inline bool referredExternal() const;
// Is this a "real" point on this processor, i.e. is it internal or part
// of the boundary description, and not a "far" or "referred" point
inline bool real() const;
inline Foam::scalar targetCellSize() const
{
return targetCellSize_;
}
// For referred vertices, what is the original processor index
inline int procIndex() const;
inline static int encodeProcIndex(int procI);
inline bool uninitialised() const
{
return type_ == vtInternal && index_ == vtInternal;
}
//- Set the point to be internal
inline void setInternal();
//- Is point internal and near the boundary
inline bool nearBoundary() const;
//- Is point a far-point
inline bool farPoint() const
{
return type_ == vtFar;
}
//- Set the point to be near the boundary
inline void setNearBoundary();
//- Either master or slave of pointPair.
inline bool pairPoint() const;
//- Is point internal, i.e. not on boundary
inline bool internalPoint() const
{
return internalPoint(type_);
}
//- Master of a pointPair is the lowest numbered one.
inline bool ppMaster() const;
//- Master of a pointPair is the lowest numbered one, external query.
inline static bool ppMaster(int index, int type);
//- Is point internal, i.e. not on boundary, external query.
inline static bool internalPoint(int type)
{
return type <= vtInternal;
}
//- Slave of a pointPair is the highest numbered one.
inline bool ppSlave() const;
//- Either original internal point or master of pointPair.
inline bool internalOrBoundaryPoint() const;
// is this a referred vertex
inline bool referred() const
{
return (type_ < 0 && type_ > vtFar);
}
//- Either original internal point or master of pointPair.
// External query.
inline static bool internalOrBoundaryPoint(int index, int type);
//- Is point near the boundary or part of the boundary definition
inline bool nearOrOnBoundary() const;
// is this a referred internal or boundary vertex
inline bool referredInternalOrBoundaryPoint() const
{
return referred() && index_ >= 0;
}
//- Either a real or referred internal or boundary point
inline bool anyInternalOrBoundaryPoint() const;
//- Is the vertex fixed or movable
inline bool isVertexFixed() const;
// is this a referred external (pair slave) vertex
inline bool referredExternal() const
{
return referred() && index_ < 0;
}
// is this a "real" point on this processor, i.e. is it internal or part of
// the boundary description, and not a "far" or "referred" point
inline bool real() const
{
return internalPoint() || pairPoint();
}
// For referred vertices, what is the original processor index
inline int procIndex() const
{
if (referred())
{
return -(type_ + 1);
}
else
{
return -1;
}
}
inline static int encodeProcIndex(int procI)
{
return -(procI + 1);
}
//- Set the point to be internal
inline void setInternal()
{
type_ = vtInternal;
}
//- Is point internal and near the boundary
inline bool nearBoundary() const
{
return type_ == vtNearBoundary;
}
//- Set the point to be near the boundary
inline void setNearBoundary()
{
type_ = vtNearBoundary;
}
//- Either master or slave of pointPair.
inline bool pairPoint() const
{
return type_ >= 0;
}
//- Master of a pointPair is the lowest numbered one.
inline bool ppMaster() const
{
return ppMaster(index_, type_);
}
//- Master of a pointPair is the lowest numbered one, external query.
inline static bool ppMaster(int index, int type)
{
if (index >= 0 && type > index)
{
return true;
}
return false;
}
//- Slave of a pointPair is the highest numbered one.
inline bool ppSlave() const
{
if (type_ >= 0 && type_ < index_)
{
return true;
}
else
{
return false;
}
}
//- Either original internal point or master of pointPair.
inline bool internalOrBoundaryPoint() const
{
return internalOrBoundaryPoint(index_, type_);
}
//- Either original internal point or master of pointPair, external query.
inline static bool internalOrBoundaryPoint(int index, int type)
{
return internalPoint(type) || ppMaster(index, type);
}
//- Is point near the boundary or part of the boundary definition
inline bool nearOrOnBoundary() const
{
return pairPoint() || nearBoundary();
}
//- Either a real or referred internal or boundary point
inline bool anyInternalOrBoundaryPoint() const
{
return internalOrBoundaryPoint() || referredInternalOrBoundaryPoint();
}
//- Fix the vertex so that it can't be moved
inline void setVertexFixed() const;
// inline void operator=(const Delaunay::Finite_vertices_iterator vit)
// {
@ -414,6 +264,8 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "indexedVertexI.H"
#ifdef NoRepository
# include "indexedVertex.C"
#endif

340
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H Normal file
View File

@ -0,0 +1,340 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
As a special exception, you have permission to link this program with the
CGAL library and distribute executables, as long as you follow the
requirements of the GNU GPL in regard to all of the software in the
executable aside from CGAL.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Gt, class Vb>
inline CGAL::indexedVertex<Gt, Vb>::indexedVertex()
:
Vb(),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0),
vertexFixed_(false)
{}
template<class Gt, class Vb>
inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p)
:
Vb(p),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0),
vertexFixed_(false)
{}
template<class Gt, class Vb>
inline CGAL::indexedVertex<Gt, Vb>::indexedVertex
(
const Point& p,
int index,
int type
)
:
Vb(p),
index_(index),
type_(type),
alignment_(),
targetCellSize_(0.0),
vertexFixed_(false)
{}
template<class Gt, class Vb>
inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p, Cell_handle f)
:
Vb(f, p),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0),
vertexFixed_(false)
{}
template<class Gt, class Vb>
inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(Cell_handle f)
:
Vb(f),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0),
vertexFixed_(false)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Gt, class Vb>
inline int& CGAL::indexedVertex<Gt, Vb>::index()
{
return index_;
}
template<class Gt, class Vb>
inline int CGAL::indexedVertex<Gt, Vb>::index() const
{
return index_;
}
template<class Gt, class Vb>
inline int& CGAL::indexedVertex<Gt, Vb>::type()
{
return type_;
}
template<class Gt, class Vb>
inline int CGAL::indexedVertex<Gt, Vb>::type() const
{
return type_;
}
template<class Gt, class Vb>
inline Foam::tensor& CGAL::indexedVertex<Gt, Vb>::alignment()
{
return alignment_;
}
template<class Gt, class Vb>
inline const Foam::tensor& CGAL::indexedVertex<Gt, Vb>::alignment() const
{
return alignment_;
}
template<class Gt, class Vb>
inline Foam::scalar& CGAL::indexedVertex<Gt, Vb>::targetCellSize()
{
return targetCellSize_;
}
template<class Gt, class Vb>
inline Foam::scalar CGAL::indexedVertex<Gt, Vb>::targetCellSize() const
{
return targetCellSize_;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::uninitialised() const
{
return type_ == vtInternal && index_ == vtInternal;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::farPoint() const
{
return type_ == vtFar;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::internalPoint() const
{
return internalPoint(type_);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::internalPoint(int type)
{
return type <= vtInternal;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::referred() const
{
return (type_ < 0 && type_ > vtFar);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::referredInternalOrBoundaryPoint() const
{
return referred() && index_ >= 0;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::referredExternal() const
{
return referred() && index_ < 0;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::real() const
{
return internalPoint() || pairPoint();
}
template<class Gt, class Vb>
inline int CGAL::indexedVertex<Gt, Vb>::procIndex() const
{
if (referred())
{
return -(type_ + 1);
}
else
{
return -1;
}
}
template<class Gt, class Vb>
inline int CGAL::indexedVertex<Gt, Vb>::encodeProcIndex(int procI)
{
return -(procI + 1);
}
template<class Gt, class Vb>
inline void CGAL::indexedVertex<Gt, Vb>::setInternal()
{
type_ = vtInternal;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::nearBoundary() const
{
return type_ == vtNearBoundary;
}
template<class Gt, class Vb>
inline void CGAL::indexedVertex<Gt, Vb>::setNearBoundary()
{
type_ = vtNearBoundary;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::pairPoint() const
{
return type_ >= 0;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::ppMaster() const
{
return ppMaster(index_, type_);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::ppMaster(int index, int type)
{
if (index >= 0 && type > index)
{
return true;
}
return false;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::ppSlave() const
{
if (type_ >= 0 && type_ < index_)
{
return true;
}
else
{
return false;
}
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::internalOrBoundaryPoint() const
{
return internalOrBoundaryPoint(index_, type_);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::internalOrBoundaryPoint
(
int index,
int type
)
{
return internalPoint(type) || ppMaster(index, type);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::nearOrOnBoundary() const
{
return pairPoint() || nearBoundary();
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::anyInternalOrBoundaryPoint() const
{
return internalOrBoundaryPoint() || referredInternalOrBoundaryPoint();
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::isVertexFixed() const
{
return vertexFixed_;
}
template<class Gt, class Vb>
inline void CGAL::indexedVertex<Gt, Vb>::setVertexFixed() const
{
vertexFixed_ = true;
}
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //

87
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointFeatureEdgesTypes.H Normal file
View File

@ -0,0 +1,87 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
As a special exception, you have permission to link this program with the
CGAL library and distribute executables, as long as you follow the
requirements of the GNU GPL in regard to all of the software in the
executable aside from CGAL.
Class
pointFeatureEdgesTypes
Description
struct for holding information on the types of feature edges attached to
feature points
\*---------------------------------------------------------------------------*/
#ifndef pointFeatureEdgesTypes_H
#define pointFeatureEdgesTypes_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class pointFeatureEdgesTypes Declaration
\*---------------------------------------------------------------------------*/
//- Hold the types of feature edges attached to the point.
struct pointFeatureEdgesTypes
{
label ptI;
label nExternal, nInternal;
label nFlat, nOpen, nMultiple, nNonFeature;
pointFeatureEdgesTypes(const label pointLabel)
{
ptI = pointLabel;
nExternal = nInternal = 0;
nFlat = nOpen = nMultiple = nNonFeature = 0;
}
friend Ostream& operator<<(Ostream& os, const pointFeatureEdgesTypes& p)
{
os << "E " << "I " << "F " << "O " << "M " << "N "
<< "Pt" << nl
<< p.nExternal << " " << p.nInternal << " "
<< p.nFlat << " " << p.nOpen << " "
<< p.nMultiple << " " << p.nNonFeature << " "
<< p.ptI
<< endl;
return os;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

43
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
View File

@ -730,6 +730,49 @@ void Foam::conformationSurfaces::findEdgeNearestByType
}
void Foam::conformationSurfaces::findAllNearestEdges
(
const point& sample,
const scalar searchRadiusSqr,
List<List<pointIndexHit> >& edgeHitsByFeature,
List<label>& featuresHit
) const
{
// Initialise
//featuresHit.setSize(features_.size());
//featuresHit = -1;
//edgeHitsByFeature.setSize(features_.size());
// Work arrays
List<pointIndexHit> hitInfo(extendedFeatureEdgeMesh::nEdgeTypes);
forAll(features_, testI)
{
features_[testI].allNearestFeatureEdges
(
sample,
searchRadiusSqr,
hitInfo
);
bool anyHit = false;
forAll(hitInfo, hitI)
{
if (hitInfo[hitI].hit())
{
anyHit = true;
}
}
if (anyHit)
{
edgeHitsByFeature.append(hitInfo);
featuresHit.append(testI);
}
}
}
void Foam::conformationSurfaces::writeFeatureObj(const fileName& prefix) const
{
OFstream ftStr(runTime_.time().path()/prefix + "_allFeatures.obj");

16
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
View File

@ -166,7 +166,7 @@ public:
//- Check if point is closer to the surfaces to conform to than
// testDistSqr, in which case return false, otherwise assess in or
// outside and erturn a result depending on the testForInside flag
// outside and return a result depending on the testForInside flag
Field<bool> wellInOutSide
(
const pointField& samplePts,
@ -228,7 +228,7 @@ public:
label& hitSurface
) const;
//- Find the nearest point to the sample and return it to the
//- Find the nearest point to the sample and return it to the
// pointIndexHit
void findSurfaceNearest
(
@ -270,7 +270,17 @@ public:
const point& sample,
scalar nearestDistSqr,
List<pointIndexHit>& edgeHit,
List<label>& featureHit
List<label>& featuresHit
) const;
//- Find the nearest points on each feature edge that is within
// a given distance from the sample point
void findAllNearestEdges
(
const point& sample,
const scalar searchRadiusSqr,
List<List<pointIndexHit> >& edgeHitsByFeature,
List<label>& featuresHit
) const;
//- Find which patch is intersected by the line from one point to

12
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
View File

@ -30,11 +30,9 @@ License
Foam::cvControls::cvControls
(
const conformalVoronoiMesh& cvMesh,
const dictionary& cvMeshDict
)
:
cvMesh_(cvMesh),
cvMeshDict_(cvMeshDict)
{
// Surface conformation controls
@ -61,6 +59,11 @@ Foam::cvControls::cvControls
surfDict.lookup("featureEdgeExclusionDistanceCoeff")
);
specialiseFeaturePoints_ = Switch
(
surfDict.lookupOrDefault<Switch>("specialiseFeaturePoints", false)
);
surfaceSearchDistanceCoeff_ = readScalar
(
surfDict.lookup("surfaceSearchDistanceCoeff")
@ -96,6 +99,11 @@ Foam::cvControls::cvControls
coarseDict.subDict("iteration")
);
surfacePtExclusionDistanceCoeff_ = readScalar
(
coarseInitialDict.lookup("surfacePtExclusionDistanceCoeff")
);
edgeSearchDistCoeffSqr_coarse_initial_ = sqr
(
readScalar

18
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
View File

@ -54,9 +54,6 @@ class cvControls
{
// Private data
//- Reference to the conformalVoronoiMesh holding this cvControls object
const conformalVoronoiMesh& cvMesh_;
//- Reference to the cvMeshDict
const dictionary& cvMeshDict_;
@ -80,6 +77,9 @@ class cvControls
// fraction of the local target cell size
scalar featureEdgeExclusionDistanceCoeff_;
//- Switch for using specialised feature points
Switch specialiseFeaturePoints_;
//- Surface search distance coefficient - fraction of the local
// target cell size
scalar surfaceSearchDistanceCoeff_;
@ -99,6 +99,11 @@ class cvControls
// Controls for coarse surface conformation
//- Distance to an existing surface conformation point location
// within which other surface point locations are excluded
// - fraction of the local target cell size
scalar surfacePtExclusionDistanceCoeff_;
//- Distance to search for feature edges near to
// surface protrusions - fraction of the local target
// cell size. Coarse conformation, initial protrusion
@ -291,7 +296,6 @@ public:
//- Construct from references to conformalVoronoiMesh and dictionary
cvControls
(
const conformalVoronoiMesh& cvMesh,
const dictionary& cvMeshDict
);
@ -319,6 +323,12 @@ public:
//- Return the featureEdgeExclusionDistanceCoeff
inline scalar featureEdgeExclusionDistanceCoeff() const;
//- Return the surfacePtExclusionDistanceCoeff
inline scalar surfacePtExclusionDistanceCoeff() const;
//- Return whether to use specialised feature points
inline Switch specialiseFeaturePoints() const;
//- Return the surfaceSearchDistanceCoeff
inline scalar surfaceSearchDistanceCoeff() const;

9
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
View File

@ -54,6 +54,15 @@ inline Foam::scalar Foam::cvControls::featureEdgeExclusionDistanceCoeff() const
return featureEdgeExclusionDistanceCoeff_;
}
inline Foam::scalar Foam::cvControls::surfacePtExclusionDistanceCoeff() const
{
return surfacePtExclusionDistanceCoeff_;
}
inline Foam::Switch Foam::cvControls::specialiseFeaturePoints() const
{
return specialiseFeaturePoints_;
}
inline Foam::scalar Foam::cvControls::surfaceSearchDistanceCoeff() const
{

9
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.C
View File

@ -42,12 +42,10 @@ defineRunTimeSelectionTable(faceAreaWeightModel, dictionary);
faceAreaWeightModel::faceAreaWeightModel
(
const word& type,
const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh
const dictionary& relaxationDict
)
:
dictionary(relaxationDict),
cvMesh_(cvMesh),
coeffDict_(subDict(type + "Coeffs"))
{}
@ -56,8 +54,7 @@ faceAreaWeightModel::faceAreaWeightModel
autoPtr<faceAreaWeightModel> faceAreaWeightModel::New
(
const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh
const dictionary& relaxationDict
)
{
word faceAreaWeightModelTypeName
@ -85,7 +82,7 @@ autoPtr<faceAreaWeightModel> faceAreaWeightModel::New
<< exit(FatalError);
}
return autoPtr<faceAreaWeightModel>(cstrIter()(relaxationDict, cvMesh));
return autoPtr<faceAreaWeightModel>(cstrIter()(relaxationDict));
}

15
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.H
View File

@ -39,7 +39,6 @@ SourceFiles
#define faceAreaWeightModel_H
#include "point.H"
#include "conformalVoronoiMesh.H"
#include "dictionary.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
@ -62,9 +61,6 @@ protected:
// Protected data
//- Reference to the conformalVoronoiMesh holding this cvControls object
const conformalVoronoiMesh& cvMesh_;
//- Method coeffs dictionary
dictionary coeffDict_;
@ -93,10 +89,9 @@ public:
faceAreaWeightModel,
dictionary,
(
const dictionary& faceAreaWeightDict,
const conformalVoronoiMesh& cvMesh
const dictionary& faceAreaWeightDict
),
(faceAreaWeightDict, cvMesh)
(faceAreaWeightDict)
);
@ -106,8 +101,7 @@ public:
faceAreaWeightModel
(
const word& type,
const dictionary& faceAreaWeightDict,
const conformalVoronoiMesh& cvMesh
const dictionary& faceAreaWeightDict
);
@ -116,8 +110,7 @@ public:
//- Return a reference to the selected faceAreaWeightModel
static autoPtr<faceAreaWeightModel> New
(
const dictionary& faceAreaWeightDict,
const conformalVoronoiMesh& cvMesh
const dictionary& faceAreaWeightDict
);

5
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
View File

@ -45,11 +45,10 @@ addToRunTimeSelectionTable
piecewiseLinearRamp::piecewiseLinearRamp
(
const dictionary& faceAreaWeightDict,
const conformalVoronoiMesh& cvMesh
const dictionary& faceAreaWeightDict
)
:
faceAreaWeightModel(typeName, faceAreaWeightDict, cvMesh),
faceAreaWeightModel(typeName, faceAreaWeightDict),
lAF_(readScalar(coeffDict().lookup("lowerAreaFraction"))),
uAF_(readScalar(coeffDict().lookup("upperAreaFraction")))
{}

3
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.H
View File

@ -71,8 +71,7 @@ public:
//- Construct from components
piecewiseLinearRamp
(
const dictionary& faceAreaWeightDict,
const conformalVoronoiMesh& cvMesh
const dictionary& faceAreaWeightDict
);

34
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
View File

@ -177,7 +177,7 @@ bool Foam::autoDensity::combinedWellInside
void Foam::autoDensity::recurseAndFill
(
std::list<Vb::Point>& initialPoints,
DynamicList<Vb::Point>& initialPoints,
const treeBoundBox& bb,
label levelLimit,
word recursionName
@ -272,7 +272,7 @@ void Foam::autoDensity::recurseAndFill
bool Foam::autoDensity::fillBox
(
std::list<Vb::Point>& initialPoints,
DynamicList<Vb::Point>& initialPoints,
const treeBoundBox& bb,
bool overlapping
) const
@ -677,7 +677,7 @@ bool Foam::autoDensity::fillBox
// Pout<< "Final volume probability break accept"
// << endl;
initialPoints.push_back
initialPoints.append
(
Vb::Point(p.x(), p.y(), p.z())
);
@ -688,7 +688,7 @@ bool Foam::autoDensity::fillBox
break;
}
initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
volumeAdded += localVolume;
}
@ -800,7 +800,7 @@ bool Foam::autoDensity::fillBox
// Pout<< "Final volume probability break accept"
// << endl;
initialPoints.push_back
initialPoints.append
(
Vb::Point(p.x(), p.y(), p.z())
);
@ -811,7 +811,7 @@ bool Foam::autoDensity::fillBox
break;
}
initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
volumeAdded += localVolume;
}
@ -849,10 +849,7 @@ autoDensity::autoDensity
:
initialPointsMethod(typeName, initialPointsDict, cvMesh),
globalTrialPoints_(0),
minCellSizeLimit_
(
detailsDict().lookupOrDefault<scalar>("minCellSizeLimit", 0.0)
),
minCellSizeLimit_(readScalar(detailsDict().lookup("minCellSizeLimit"))),
minLevels_(readLabel(detailsDict().lookup("minLevels"))),
maxSizeRatio_(readScalar(detailsDict().lookup("maxSizeRatio"))),
volRes_(readLabel(detailsDict().lookup("sampleResolution"))),
@ -882,7 +879,7 @@ autoDensity::autoDensity
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
std::list<Vb::Point> autoDensity::initialPoints() const
List<Vb::Point> autoDensity::initialPoints() const
{
treeBoundBox hierBB;
@ -902,7 +899,18 @@ std::list<Vb::Point> autoDensity::initialPoints() const
);
}
std::list<Vb::Point> initialPoints;
// Initialise size of points list.
const scalar volumeBoundBox = Foam::pow3(hierBB.typDim());
const scalar volumeSmallestCell = Foam::pow3(minCellSizeLimit_);
const int initialPointEstimate
= min
(
static_cast<int>(volumeBoundBox/(volumeSmallestCell + SMALL)/10),
1e6
);
DynamicList<Vb::Point> initialPoints(initialPointEstimate);
Info<< nl << " " << typeName << endl;
@ -919,6 +927,8 @@ std::list<Vb::Point> autoDensity::initialPoints() const
"recursionBox"
);
initialPoints.shrink();
label nInitialPoints = initialPoints.size();
if (Pstream::parRun())

6
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.H
View File

@ -116,7 +116,7 @@ private:
//- Descend into octants of the supplied bound box
void recurseAndFill
(
std::list<Vb::Point>& initialPoints,
DynamicList<Vb::Point>& initialPoints,
const treeBoundBox& bb,
label levelLimit,
word recursionName
@ -127,7 +127,7 @@ private:
// in favour of further recursion.
bool fillBox
(
std::list<Vb::Point>& initialPoints,
DynamicList<Vb::Point>& initialPoints,
const treeBoundBox& bb,
bool overlapping
) const;
@ -156,7 +156,7 @@ public:
// Member Functions
//- Return the initial points for the conformalVoronoiMesh
virtual std::list<Vb::Point> initialPoints() const;
virtual List<Vb::Point> initialPoints() const;
};

8
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C
View File

@ -56,7 +56,7 @@ bodyCentredCubic::bodyCentredCubic
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
std::list<Vb::Point> bodyCentredCubic::initialPoints() const
List<Vb::Point> bodyCentredCubic::initialPoints() const
{
boundBox bb;
@ -91,7 +91,7 @@ std::list<Vb::Point> bodyCentredCubic::initialPoints() const
scalar pert = randomPerturbationCoeff_*cmptMin(delta);
std::list<Vb::Point> initialPoints;
DynamicList<Vb::Point> initialPoints(ni*nj*nk/10);
for (label i = 0; i < ni; i++)
{
@ -177,13 +177,13 @@ std::list<Vb::Point> bodyCentredCubic::initialPoints() const
{
const point& p(points[i]);
initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
}
}
}
}
return initialPoints;
return initialPoints.shrink();
}

2
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.H
View File

@ -89,7 +89,7 @@ public:
// Member Functions
//- Return the initial points for the conformalVoronoiMesh
virtual std::list<Vb::Point> initialPoints() const;
virtual List<Vb::Point> initialPoints() const;
};

8
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.C
View File

@ -56,7 +56,7 @@ faceCentredCubic::faceCentredCubic
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
std::list<Vb::Point> faceCentredCubic::initialPoints() const
List<Vb::Point> faceCentredCubic::initialPoints() const
{
boundBox bb;
@ -91,7 +91,7 @@ std::list<Vb::Point> faceCentredCubic::initialPoints() const
scalar pert = randomPerturbationCoeff_*cmptMin(delta);
std::list<Vb::Point> initialPoints;
DynamicList<Vb::Point> initialPoints(ni*nj*nk/10);
for (label i = 0; i < ni; i++)
{
@ -238,13 +238,13 @@ std::list<Vb::Point> faceCentredCubic::initialPoints() const
{
const point& p(points[i]);
initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
}
}
}
}
return initialPoints;
return initialPoints.shrink();
}

2
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.H
View File

@ -89,7 +89,7 @@ public:
// Member Functions
//- Return the initial points for the conformalVoronoiMesh
virtual std::list<Vb::Point> initialPoints() const;
virtual List<Vb::Point> initialPoints() const;
};

2
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/initialPointsMethod/initialPointsMethod.H
View File

@ -139,7 +139,7 @@ public:
}
//- Return the initial points for the conformalVoronoiMesh
virtual std::list<Vb::Point> initialPoints() const = 0;
virtual List<Vb::Point> initialPoints() const = 0;
};

12
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
View File

@ -51,7 +51,7 @@ pointFile::pointFile
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
std::list<Vb::Point> pointFile::initialPoints() const
List<Vb::Point> pointFile::initialPoints() const
{
pointIOField points
(
@ -69,7 +69,7 @@ std::list<Vb::Point> pointFile::initialPoints() const
if (points.empty())
{
FatalErrorIn("std::list<Vb::Point> pointFile::initialPoints() const")
FatalErrorIn("List<Vb::Point> pointFile::initialPoints() const")
<< "Point file contain no points"
<< exit(FatalError) << endl;
}
@ -126,8 +126,6 @@ std::list<Vb::Point> pointFile::initialPoints() const
}
}
std::list<Vb::Point> initialPoints;
Field<bool> insidePoints = cvMesh_.geometryToConformTo().wellInside
(
points,
@ -138,16 +136,20 @@ std::list<Vb::Point> pointFile::initialPoints() const
)
);
DynamicList<Vb::Point> initialPoints(insidePoints.size()/10);
forAll(insidePoints, i)
{
if (insidePoints[i])
{
const point& p(points[i]);
initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
}
}
initialPoints.shrink();
label nPointsRejected = points.size() - initialPoints.size();
if (Pstream::parRun())

2
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.H
View File

@ -85,7 +85,7 @@ public:
// Member Functions
//- Return the initial points for the conformalVoronoiMesh
virtual std::list<Vb::Point> initialPoints() const;
virtual List<Vb::Point> initialPoints() const;
};

9
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
View File

@ -56,7 +56,7 @@ uniformGrid::uniformGrid
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
std::list<Vb::Point> uniformGrid::initialPoints() const
List<Vb::Point> uniformGrid::initialPoints() const
{
boundBox bb;
@ -91,7 +91,8 @@ std::list<Vb::Point> uniformGrid::initialPoints() const
scalar pert = randomPerturbationCoeff_*cmptMin(delta);
std::list<Vb::Point> initialPoints;
// Initialise points list
DynamicList<Vb::Point> initialPoints(ni*nj*nk/10);
for (label i = 0; i < ni; i++)
{
@ -153,13 +154,13 @@ std::list<Vb::Point> uniformGrid::initialPoints() const
{
const point& p(points[i]);
initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
}
}
}
}
return initialPoints;
return initialPoints.shrink();
}

2
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.H
View File

@ -89,7 +89,7 @@ public:
// Member Functions
//- Return the initial points for the conformalVoronoiMesh
virtual std::list<Vb::Point> initialPoints() const;
virtual List<Vb::Point> initialPoints() const;
};

4
applications/utilities/mesh/generation/cvMesh/cvMesh.C
View File

@ -59,6 +59,8 @@ int main(int argc, char *argv[])
)
);
conformalVoronoiMesh::debug = true;
conformalVoronoiMesh mesh(runTime, cvMeshDict);
while (runTime.loop())
@ -72,7 +74,7 @@ int main(int argc, char *argv[])
<< nl << endl;
}
mesh.writeMesh(runTime.constant());
mesh.writeMesh(runTime.constant(), true);
Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"

157
applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H Normal file
View File

@ -0,0 +1,157 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
Class
Foam::vectorTools
Description
Functions for analysing the relationships between vectors
SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef vectorTools_H
#define vectorTools_H
#include "vector.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Namespace vectorTools Declaration
\*---------------------------------------------------------------------------*/
//- Collection of functions for testing relationships between two vectors.
namespace vectorTools
{
//- Test if a and b are parallel: a.b = 1
// Uses the cross product, so the tolerance is proportional to
// the sine of the angle between a and b in radians
template <typename T>
bool areParallel
(
const Vector<T>& a,
const Vector<T>& b,
const T& tolerance = SMALL
)
{
return (mag(a ^ b) < tolerance) ? true : false;
// return ( mag( mag(a & b)/(mag(a)*mag(b)) - 1.0 ) < tolerance )
// ? true
// : false;
}
//- Test if a and b are orthogonal: a.b = 0
// Uses the dot product, so the tolerance is proportional to
// the cosine of the angle between a and b in radians
template <typename T>
bool areOrthogonal
(
const Vector<T>& a,
const Vector<T>& b,
const T& tolerance = SMALL
)
{
return (mag(a & b) < tolerance) ? true : false;
}
//- Test if angle between a and b is acute: a.b > 0
template <typename T>
bool areAcute
(
const Vector<T>& a,
const Vector<T>& b
)
{
return ((a & b) > 0) ? true : false;
}
//- Test if angle between a and b is obtuse: a.b < 0
template <typename T>
bool areObtuse
(
const Vector<T>& a,
const Vector<T>& b
)
{
return ((a & b) < 0) ? true : false;
}
//- Calculate angle between a and b in radians
template <typename T>
T cosPhi
(
const Vector<T>& a,
const Vector<T>& b,
const T& tolerance = SMALL
)
{
scalar cosPhi = (a & b)/(mag(a)*mag(b) + tolerance);
// Enforce bounding between -1 and 1
return min(max(cosPhi, -1), 1);
}
//- Calculate angle between a and b in radians
template <typename T>
T radAngleBetween
(
const Vector<T>& a,
const Vector<T>& b,
const T& tolerance = SMALL
)
{
scalar cosPhi = (a & b)/(mag(a)*mag(b) + tolerance);
// Enforce bounding between -1 and 1
return acos( min(max(cosPhi, -1), 1) );
}
//- Calculate angle between a and b in degrees
template <typename T>
T degAngleBetween
(
const Vector<T>& a,
const Vector<T>& b,
const T& tolerance = SMALL
)
{
return radToDeg(radAngleBetween(a, b, tolerance));
}
} // End namespace vectorTools
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

49
applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
View File

@ -8,6 +8,7 @@
#include "pyrMatcher.H"
#include "tetWedgeMatcher.H"
#include "tetMatcher.H"
#include "IOmanip.H"
void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
@ -62,22 +63,19 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
<< nInternalEdges-nInternal1Edges << nl;
}
Info<< " faces: "
<< returnReduce(mesh.faces().size(), sumOp<label>()) << nl
<< " internal faces: "
<< returnReduce(mesh.faceNeighbour().size(), sumOp<label>()) << nl
<< " cells: "
<< returnReduce(mesh.cells().size(), sumOp<label>()) << nl
<< " boundary patches: "
<< mesh.boundaryMesh().size() << nl
<< " point zones: "
<< mesh.pointZones().size() << nl
<< " face zones: "
<< mesh.faceZones().size() << nl
<< " cell zones: "
<< mesh.cellZones().size() << nl
<< endl;
label nFaces = returnReduce(mesh.faces().size(), sumOp<label>());
label nIntFaces = returnReduce(mesh.faceNeighbour().size(), sumOp<label>());
label nCells = returnReduce(mesh.cells().size(), sumOp<label>());
Info<< " faces: " << nFaces << nl
<< " internal faces: " << nIntFaces << nl
<< " cells: " << nCells << nl
<< " faces per cell: " << scalar(nFaces + nIntFaces)/nCells << nl
<< " boundary patches: " << mesh.boundaryMesh().size() << nl
<< " point zones: " << mesh.pointZones().size() << nl
<< " face zones: " << mesh.faceZones().size() << nl
<< " cell zones: " << mesh.cellZones().size() << nl
<< endl;
// Construct shape recognizers
hexMatcher hex;
@ -96,6 +94,8 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
label nTetWedge = 0;
label nUnknown = 0;
Map<label> polyhedralFaces;
for (label cellI = 0; cellI < mesh.nCells(); cellI++)
{
if (hex.isA(mesh, cellI))
@ -125,6 +125,7 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
else
{
nUnknown++;
polyhedralFaces(mesh.cells()[cellI].size())++;
}
}
@ -144,5 +145,21 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
<< " tet wedges: " << nTetWedge << nl
<< " tetrahedra: " << nTet << nl
<< " polyhedra: " << nUnknown
<< nl << endl;
<< endl;
if (nUnknown > 0)
{
Pstream::mapCombineGather(polyhedralFaces, plusEqOp<label>());
Info<< " Breakdown of polyhedra by number of faces:" << endl;
Info<< " faces" << " number of cells" << endl;
forAllConstIter(Map<label>, polyhedralFaces, iter)
{
Info<< setf(std::ios::right) << setw(13)
<< iter.key() << " " << iter() << nl;
}
}
Info<< endl;
}