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

Merge branch 'master' of ssh://opencfd:8007/home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
Henry
2013-09-26 22:42:05 +01:00
parent 3004606a5a f663d40545
commit 88f82124a2
73 changed files with 2926 additions and 1506 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
applications/test/PtrList/Test-PtrList.C
View File

@ -32,6 +32,8 @@ Description
#include "scalar.H"
#include "IOstreams.H"
#include "PtrList.H"
#include "plane.H"
#include "DynamicList.H"
using namespace Foam;
@ -76,6 +78,7 @@ int main(int argc, char *argv[])
{
PtrList<Scalar> list1(10);
PtrList<Scalar> list2(15);
PtrList<Scalar> listApp;
forAll(list1, i)
{
@ -87,9 +90,14 @@ int main(int argc, char *argv[])
list2.set(i, new Scalar(10 + 1.3*i));
}
for (label i = 0; i < 5; ++i)
{
listApp.append(new Scalar(1.3*i));
}
Info<<"list1: " << list1 << endl;
Info<<"list2: " << list2 << endl;
Info<<"listApp: " << listApp << endl;
Info<<"indirectly delete some items via set(.., 0) :" << endl;
for (label i = 0; i < 3; i++)
@ -115,6 +123,13 @@ int main(int argc, char *argv[])
Info<<"list2: " << list2 << endl;
Info<<"list3: " << list3 << endl;
PtrList<plane> planes;
planes.append(new plane(vector::one, vector::one));
planes.append(new plane(vector(1,2,3), vector::one));
forAll(planes, p)
Info<< "plane " << planes[p] << endl;
Info<< nl << "Done." << endl;
return 0;
}

7
applications/test/sort/Test-sortList.C
View File

@ -50,11 +50,18 @@ int main(int argc, char *argv[])
labelList a(orig);
sortedOrder(a, order);
SortableList<label> aReverse(a.size());
aReverse = a;
Info<< "unsorted: " << a << endl;
sort(a);
Info<< "sorted: " << a << endl;
Info<< "indices: " << order << endl;
aReverse.reverseSort();
Info<< "reverse sorted: " << aReverse << endl;
Info<< "reverse indices: " << aReverse.indices() << endl;
SortableList<label> b(orig);
Info<< "unsorted: " << orig << endl;
Info<< "sorted: " << b << endl;

45
applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
View File

@ -52,6 +52,7 @@ Usage
#include "polyTopoChange.H"
#include "fvMesh.H"
#include "polyMeshFilter.H"
#include "faceSet.H"
using namespace Foam;
@ -74,9 +75,9 @@ int main(int argc, char *argv[])
argList::addOption
(
"collapseFaceZone",
"zoneName",
"Collapse faces that are in the supplied face zone"
"collapseFaceSet",
"faceSet",
"Collapse faces that are in the supplied face set"
);
# include "addOverwriteOption.H"
@ -93,15 +94,15 @@ int main(int argc, char *argv[])
const bool overwrite = args.optionFound("overwrite");
const bool collapseFaces = args.optionFound("collapseFaces");
const bool collapseFaceZone = args.optionFound("collapseFaceZone");
const bool collapseFaceSet = args.optionFound("collapseFaceSet");
if (collapseFaces && collapseFaceZone)
if (collapseFaces && collapseFaceSet)
{
FatalErrorIn("main(int, char*[])")
<< "Both face zone collapsing and face collapsing have been"
<< "selected. Choose only one of:" << nl
<< " -collapseFaces" << nl
<< " -collapseFaceZone <faceZoneName>"
<< " -collapseFaceSet <faceSet>"
<< abort(FatalError);
}
@ -117,7 +118,6 @@ int main(int argc, char *argv[])
),
labelList(mesh.nPoints(), labelMin)
);
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
@ -128,6 +128,14 @@ int main(int argc, char *argv[])
label nBadFaces = 0;
faceSet indirectPatchFaces
(
mesh,
"indirectPatchFaces",
IOobject::MUST_READ,
IOobject::AUTO_WRITE
);
{
meshFilterPtr.set(new polyMeshFilter(mesh, pointPriority));
polyMeshFilter& meshFilter = meshFilterPtr();
@ -139,19 +147,26 @@ int main(int argc, char *argv[])
// the face filtering is sped up.
nBadFaces = meshFilter.filterEdges(0);
{
polyTopoChange meshMod(newMesh);
polyTopoChange meshMod(newMesh());
meshMod.changeMesh(mesh, false);
polyMeshFilter::copySets(newMesh(), mesh);
}
pointPriority = meshFilter.pointPriority();
}
if (collapseFaceZone)
if (collapseFaceSet)
{
const word faceZoneName = args.optionRead<word>("collapseFaceZone");
const faceZone& fZone = mesh.faceZones()[faceZoneName];
const word faceSetName(args.optionRead<word>("collapseFaceSet"));
faceSet fSet
(
mesh,
faceSetName,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
);
meshFilterPtr.reset(new polyMeshFilter(mesh, pointPriority));
polyMeshFilter& meshFilter = meshFilterPtr();
@ -160,11 +175,13 @@ int main(int argc, char *argv[])
// Filter faces. Pass in the number of bad faces that are present
// from the previous edge filtering to use as a stopping criterion.
meshFilter.filter(fZone);
meshFilter.filter(fSet);
{
polyTopoChange meshMod(newMesh);
meshMod.changeMesh(mesh, false);
polyMeshFilter::copySets(newMesh(), mesh);
}
pointPriority = meshFilter.pointPriority();
@ -184,6 +201,8 @@ int main(int argc, char *argv[])
polyTopoChange meshMod(newMesh);
meshMod.changeMesh(mesh, false);
polyMeshFilter::copySets(newMesh(), mesh);
}
pointPriority = meshFilter.pointPriority();

2
applications/utilities/mesh/generation/foamyHexMesh/Make/options
View File

@ -10,8 +10,10 @@ include $(GENERAL_RULES)/CGAL
EXE_INC = \
${EXE_FROUNDING_MATH} \
${EXE_NDEBUG} \
${CGAL_EXACT} \
${CGAL_INEXACT} \
${CGAL_INC} \
${c++CGALWARN} \
-IconformalVoronoiMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \

29
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C
View File

@ -190,20 +190,25 @@ void Foam::DelaunayMesh<Triangulation>::reset()
resetVertexCount();
resetCellCount();
insertPoints(vertices);
insertPoints(vertices, false);
Info<< "Inserted " << vertexCount() << " fixed points" << endl;
}
template<class Triangulation>
void Foam::DelaunayMesh<Triangulation>::insertPoints(const List<Vb>& vertices)
Foam::Map<Foam::label> Foam::DelaunayMesh<Triangulation>::insertPoints
(
const List<Vb>& vertices,
const bool reIndex
)
{
rangeInsertWithInfo
return rangeInsertWithInfo
(
vertices.begin(),
vertices.end(),
false
false,
reIndex
);
}
@ -268,7 +273,7 @@ const
template<class Triangulation>
template<class PointIterator>
void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
Foam::Map<Foam::label> Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
(
PointIterator begin,
PointIterator end,
@ -307,6 +312,10 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
Vertex_handle hint;
Map<label> oldToNewIndex(points.size());
label maxIndex = -1;
for
(
typename vectorPairPointIndex::const_iterator p = points.begin();
@ -335,11 +344,14 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
{
if (reIndex)
{
const label oldIndex = vert.index();
hint->index() = getNewVertexIndex();
oldToNewIndex.insert(oldIndex, hint->index());
}
else
{
hint->index() = vert.index();
maxIndex = max(maxIndex, vert.index());
}
hint->type() = vert.type();
hint->procIndex() = vert.procIndex();
@ -347,6 +359,13 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
hint->alignment() = vert.alignment();
}
}
if (!reIndex)
{
vertexCount_ = maxIndex + 1;
}
return oldToNewIndex;
}

8
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H
View File

@ -228,7 +228,11 @@ public:
void reset();
//- Insert the list of vertices (calls rangeInsertWithInfo)
void insertPoints(const List<Vb>& vertices);
Map<label> insertPoints
(
const List<Vb>& vertices,
const bool reIndex
);
//- Function inserting points into a triangulation and setting the
// index and type data of the point in the correct order. This is
@ -237,7 +241,7 @@ public:
// Adapted from a post on the CGAL lists: 2010-01/msg00004.html by
// Sebastien Loriot (Geometry Factory).
template<class PointIterator>
void rangeInsertWithInfo
Map<label> rangeInsertWithInfo
(
PointIterator begin,
PointIterator end,

3
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
View File

@ -536,7 +536,8 @@ Foam::label Foam::DistributedDelaunayMesh<Triangulation>::referVertices
labelPairHashSet pointsNotInserted = rangeInsertReferredWithInfo
(
referredVertices.begin(),
referredVertices.end()
referredVertices.end(),
true
);
if (!pointsNotInserted.empty())

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

@ -76,6 +76,7 @@ faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
searchableSurfaceFeatures/searchableSurfaceFeatures.C
searchableSurfaceFeatures/searchableBoxFeatures.C
searchableSurfaceFeatures/searchablePlateFeatures.C
searchableSurfaceFeatures/triSurfaceMeshFeatures.C
LIB = $(FOAM_LIBBIN)/libconformalVoronoiMesh

2
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options
View File

@ -11,8 +11,10 @@ FFLAGS = -DCGAL_FILES='"${CGAL_ARCH_PATH}/share/files"'
EXE_INC = \
${EXE_FROUNDING_MATH} \
${EXE_NDEBUG} \
${CGAL_EXACT} \
${CGAL_INEXACT} \
${CGAL_INC} \
${c++CGALWARN} \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \

117
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
View File

@ -628,49 +628,6 @@ Foam::tensorField Foam::cellShapeControlMesh::dumpAlignments() const
}
void Foam::cellShapeControlMesh::insertBoundingPoints
(
const boundBox& bb,
const cellSizeAndAlignmentControls& sizeControls
)
{
// Loop over bound box points and get cell size and alignment
const pointField bbPoints(bb.points());
forAll(bbPoints, pI)
{
const Foam::point& pt = bbPoints[pI];
// Cell size here will return default cell size
const scalar cellSize = sizeControls.cellSize(pt);
if (debug)
{
Info<< "Insert Bounding Point: " << pt << " " << cellSize << endl;
}
// Get the cell size of the nearest surface.
// geometryToConformTo_.findSurfaceNearest
// (
// pt,
// GREAT,
// surfHit,
// hitSurface
// );
const tensor alignment = tensor::I;
insert
(
pt,
cellSize,
alignment,
Vb::vtInternalNearBoundary
);
}
}
void Foam::cellShapeControlMesh::write() const
{
Info<< "Writing " << meshSubDir << endl;
@ -781,4 +738,78 @@ void Foam::cellShapeControlMesh::write() const
}
Foam::label Foam::cellShapeControlMesh::estimateCellCount
(
const autoPtr<backgroundMeshDecomposition>& decomposition
) const
{
// Loop over all the tets and estimate the cell count in each one
scalar cellCount = 0;
for
(
Finite_cells_iterator cit = finite_cells_begin();
cit != finite_cells_end();
++cit
)
{
if (!cit->hasFarPoint() && !is_infinite(cit))
{
// @todo Check if tet centre is on the processor..
CGAL::Tetrahedron_3<baseK> tet
(
cit->vertex(0)->point(),
cit->vertex(1)->point(),
cit->vertex(2)->point(),
cit->vertex(3)->point()
);
pointFromPoint centre = topoint(CGAL::centroid(tet));
if
(
Pstream::parRun()
&& !decomposition().positionOnThisProcessor(centre)
)
{
continue;
}
scalar volume = CGAL::to_double(tet.volume());
scalar averagedPointCellSize = 0;
//scalar averagedPointCellSize = 1;
// Get an average volume by averaging the cell size of the vertices
for (label vI = 0; vI < 4; ++vI)
{
averagedPointCellSize += cit->vertex(vI)->targetCellSize();
//averagedPointCellSize *= cit->vertex(vI)->targetCellSize();
}
averagedPointCellSize /= 4;
//averagedPointCellSize = ::sqrt(averagedPointCellSize);
// if (averagedPointCellSize < SMALL)
// {
// Pout<< "Volume = " << volume << endl;
//
// for (label vI = 0; vI < 4; ++vI)
// {
// Pout<< "Point " << vI
// << ", point = " << topoint(cit->vertex(vI)->point())
// << ", size = " << cit->vertex(vI)->targetCellSize()
// << endl;
// }
// }
cellCount += volume/pow(averagedPointCellSize, 3);
}
}
return cellCount;
}
// ************************************************************************* //

11
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.H
View File

@ -153,15 +153,14 @@ public:
tensorField dumpAlignments() const;
void insertBoundingPoints
(
const boundBox& bb,
const cellSizeAndAlignmentControls& sizeControls
);
void writeTriangulation();
void write() const;
label estimateCellCount
(
const autoPtr<backgroundMeshDecomposition>& decomposition
) const;
};

5
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
View File

@ -429,6 +429,11 @@ void Foam::searchableSurfaceControl::initialVertices
vectorField normals(1);
searchableSurface_.getNormal(infoList, normals);
if (mag(normals[0]) < SMALL)
{
normals[0] = vector(1, 1, 1);
}
pointAlignment.set(new triad(normals[0]));
if (infoList[0].hit())

2
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C
View File

@ -768,7 +768,7 @@ Foam::label Foam::controlMeshRefinement::refineMesh
}
}
mesh_.insertPoints(verts);
mesh_.insertPoints(verts, false);
return verts.size();
}

197
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
View File

@ -179,15 +179,23 @@ void Foam::conformalVoronoiMesh::insertInternalPoints
}
void Foam::conformalVoronoiMesh::insertPoints
Foam::Map<Foam::label> Foam::conformalVoronoiMesh::insertPointPairs
(
List<Vb>& vertices,
bool distribute
bool distribute,
bool reIndex
)
{
if (Pstream::parRun() && distribute)
{
decomposition_().distributePoints(vertices);
autoPtr<mapDistribute> mapDist =
decomposition_().distributePoints(vertices);
// Re-index the point pairs if one or both have been distributed.
// If both, remove
// If added a point, then need to know its point pair
// If one moved, then update procIndex locally
forAll(vertices, vI)
{
@ -197,7 +205,8 @@ void Foam::conformalVoronoiMesh::insertPoints
label preReinsertionSize(number_of_vertices());
this->DelaunayMesh<Delaunay>::insertPoints(vertices);
Map<label> oldToNewIndices =
this->DelaunayMesh<Delaunay>::insertPoints(vertices, reIndex);
const label nReinserted = returnReduce
(
@ -208,17 +217,18 @@ void Foam::conformalVoronoiMesh::insertPoints
Info<< " Reinserted " << nReinserted << " vertices out of "
<< returnReduce(vertices.size(), sumOp<label>())
<< endl;
return oldToNewIndices;
}
void Foam::conformalVoronoiMesh::insertSurfacePointPairs
(
const pointIndexHitAndFeatureList& surfaceHits,
const fileName fName
const fileName fName,
DynamicList<Vb>& pts
)
{
DynamicList<Vb> pts(2.0*surfaceHits.size());
forAll(surfaceHits, i)
{
vectorField norm(1);
@ -246,6 +256,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
pointPairDistance(surfacePt),
surfacePt,
normal,
true,
pts
);
}
@ -256,6 +267,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
pointPairDistance(surfacePt),
surfacePt,
normal,
true,
pts
);
}
@ -266,6 +278,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
pointPairDistance(surfacePt),
surfacePt,
-normal,
true,
pts
);
}
@ -280,8 +293,6 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
}
}
insertPoints(pts, true);
if (foamyHexMeshControls().objOutput() && fName != fileName::null)
{
DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
@ -292,11 +303,10 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
void Foam::conformalVoronoiMesh::insertEdgePointGroups
(
const pointIndexHitAndFeatureList& edgeHits,
const fileName fName
const fileName fName,
DynamicList<Vb>& pts
)
{
DynamicList<Vb> pts(3.0*edgeHits.size());
forAll(edgeHits, i)
{
if (edgeHits[i].first().hit())
@ -318,10 +328,6 @@ void Foam::conformalVoronoiMesh::insertEdgePointGroups
}
}
pts.shrink();
insertPoints(pts, true);
if (foamyHexMeshControls().objOutput() && fName != fileName::null)
{
DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
@ -348,6 +354,28 @@ bool Foam::conformalVoronoiMesh::nearFeaturePt(const Foam::point& pt) const
}
bool Foam::conformalVoronoiMesh::surfacePtNearFeatureEdge
(
const Foam::point& pt
) const
{
scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
pointIndexHit info;
label featureHit;
geometryToConformTo_.findEdgeNearest
(
pt,
exclusionRangeSqr,
info,
featureHit
);
return info.hit();
}
void Foam::conformalVoronoiMesh::insertInitialPoints()
{
Info<< nl << "Inserting initial points" << endl;
@ -859,6 +887,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
geometryToConformTo_
),
limitBounds_(),
ptPairs_(*this),
ftPtConformer_(*this),
edgeLocationTreePtr_(),
surfacePtLocationTreePtr_(),
@ -1108,7 +1137,7 @@ void Foam::conformalVoronoiMesh::move()
vB->alignment().T() & cartesianDirections
);
Field<vector> alignmentDirs(3);
Field<vector> alignmentDirs(alignmentDirsA);
forAll(alignmentDirsA, aA)
{
@ -1266,7 +1295,7 @@ void Foam::conformalVoronoiMesh::move()
if
(
(
(vA->internalPoint() || vB->internalPoint())
(vA->internalPoint() && vB->internalPoint())
&& (!vA->referred() || !vB->referred())
// ||
// (
@ -1333,12 +1362,17 @@ void Foam::conformalVoronoiMesh::move()
// Point removal
if
(
vA->internalPoint()
&& !vA->referred()
&& !vA->fixed()
&& vB->internalPoint()
&& !vB->referred()
&& !vB->fixed()
(
vA->internalPoint()
&& !vA->referred()
&& !vA->fixed()
)
&&
(
vB->internalPoint()
&& !vB->referred()
&& !vB->fixed()
)
)
{
// Only insert a point at the midpoint of
@ -1507,9 +1541,7 @@ void Foam::conformalVoronoiMesh::move()
Info<< "Writing point displacement vectors to file." << endl;
OFstream str
(
time().path()
+ "displacements_" + runTime_.timeName()
+ ".obj"
time().path()/"displacements_" + runTime_.timeName() + ".obj"
);
for
@ -1607,6 +1639,74 @@ void Foam::conformalVoronoiMesh::move()
time().path()/"boundaryPoints_" + time().timeName() + ".obj",
*this
);
DelaunayMeshTools::writeOBJ
(
time().path()/"internalBoundaryPoints_" + time().timeName()
+ ".obj",
*this,
Foam::indexedVertexEnum::vtInternalSurface,
Foam::indexedVertexEnum::vtInternalFeaturePoint
);
DelaunayMeshTools::writeOBJ
(
time().path()/"externalBoundaryPoints_" + time().timeName()
+ ".obj",
*this,
Foam::indexedVertexEnum::vtExternalSurface,
Foam::indexedVertexEnum::vtExternalFeaturePoint
);
OBJstream multipleIntersections
(
"multipleIntersections_"
+ time().timeName()
+ ".obj"
);
for
(
Delaunay::Finite_edges_iterator eit = finite_edges_begin();
eit != finite_edges_end();
++eit
)
{
Cell_handle c = eit->first;
Vertex_handle vA = c->vertex(eit->second);
Vertex_handle vB = c->vertex(eit->third);
Foam::point ptA(topoint(vA->point()));
Foam::point ptB(topoint(vB->point()));
List<pointIndexHit> surfHits;
labelList hitSurfaces;
geometryToConformTo().findSurfaceAllIntersections
(
ptA,
ptB,
surfHits,
hitSurfaces
);
if
(
surfHits.size() >= 2
|| (
surfHits.size() == 0
&& (
(vA->externalBoundaryPoint()
&& vB->internalBoundaryPoint())
|| (vB->externalBoundaryPoint()
&& vA->internalBoundaryPoint())
)
)
)
{
multipleIntersections.write(linePointRef(ptA, ptB));
}
}
}
}
@ -1629,49 +1729,6 @@ void Foam::conformalVoronoiMesh::move()
}
//Foam::labelListList Foam::conformalVoronoiMesh::overlapsProc
//(
// const List<Foam::point>& centres,
// const List<scalar>& radiusSqrs
//) const
//{
// if (!Pstream::parRun())
// {
// return labelListList(centres.size(), labelList(0));
// }
//
//// DynamicList<Foam::point> pts(number_of_vertices());
//
//// for
//// (
//// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
//// vit != finite_vertices_end();
//// vit++
//// )
//// {
//// pts.append(topoint(vit->point()));
//// }
////
//// dynamicIndexedOctree<dynamicTreeDataPoint> vertexOctree
//// (
//// dynamicTreeDataPoint(pts),
//// treeBoundBox(min(pts), max(pts)),
//// 10, // maxLevel
//// 10, // leafSize
//// 3.0 // duplicity
//// );
//
// return decomposition_().overlapsProcessors
// (
// centres,
// radiusSqrs,
// *this,
// false//,
//// vertexOctree
// );
//}
void Foam::conformalVoronoiMesh::checkCoPlanarCells() const
{
typedef CGAL::Exact_predicates_exact_constructions_kernel Kexact;

53
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
View File

@ -77,7 +77,7 @@ SourceFiles
#include "Pair.H"
#include "DistributedDelaunayMesh.H"
#include "featurePointConformer.H"
#include "pointPairs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -86,12 +86,10 @@ namespace Foam
// Forward declaration of classes
class initialPointsMethod;
class relaxationModel;
class faceAreaWeightModel;
class backgroundMeshDecomposition;
class OBJstream;
/*---------------------------------------------------------------------------*\
Class conformalVoronoiMesh Declaration
@ -165,6 +163,9 @@ private:
//- Limiting bound box before infinity begins
treeBoundBox limitBounds_;
//-
mutable pointPairs<Delaunay> ptPairs_;
//-
featurePointConformer ftPtConformer_;
@ -229,10 +230,11 @@ private:
const bool distribute = false
);
void insertPoints
Map<label> insertPointPairs
(
List<Vb>& vertices,
bool distribute
bool distribute,
bool reIndex
);
//- Create a point-pair at a ppDist distance either side of
@ -242,6 +244,7 @@ private:
const scalar ppDist,
const Foam::point& surfPt,
const vector& n,
const bool ptPair,
DynamicList<Vb>& pts
) const;
@ -254,24 +257,20 @@ private:
const scalar ppDist,
const Foam::point& surfPt,
const vector& n,
const bool ptPair,
DynamicList<Vb>& pts
) const;
//- Check internal point is completely inside the meshable region
inline bool internalPointIsInside(const Foam::point& pt) const;
inline bool isPointPair
(
const Vertex_handle& vA,
const Vertex_handle& vB
) const;
//- Insert pairs of points on the surface with the given normals, at the
// specified spacing
void insertSurfacePointPairs
(
const pointIndexHitAndFeatureList& surfaceHits,
const fileName fName = fileName::null
const fileName fName,
DynamicList<Vb>& pts
);
//- Insert groups of points to conform to an edge given a list of
@ -280,7 +279,8 @@ private:
void insertEdgePointGroups
(
const pointIndexHitAndFeatureList& edgeHits,
const fileName fName = fileName::null
const fileName fName,
DynamicList<Vb>& pts
);
void createEdgePointGroupByCirculating
@ -351,6 +351,9 @@ private:
//- Check if a location is in exclusion range around a feature point
bool nearFeaturePt(const Foam::point& pt) const;
//- Check if a surface point is in exclusion range around a feature edge
bool surfacePtNearFeatureEdge(const Foam::point& pt) const;
//- Insert the initial points into the triangulation, based on the
// initialPointsMethod
void insertInitialPoints();
@ -474,14 +477,11 @@ private:
label& hitSurface
) const;
//- Find the "worst" incursion of the dual cell of a non-internal or
// boundary point through the surface, subject to the
// maxSurfaceProtrusion tolerance
void dualCellLargestSurfaceIncursion
(
const Delaunay::Finite_vertices_iterator& vit,
pointIndexHit& surfHitLargest,
label& hitSurfaceLargest
pointIndexHit& surfHit,
label& hitSurface
) const;
//- Write out vertex-processor occupancy information for debugging
@ -695,7 +695,8 @@ private:
//- Merge vertices that are identical
void mergeIdenticalDualVertices
(
const pointField& pts
const pointField& pts,
labelList& boundaryPts
);
label mergeIdenticalDualVertices
@ -727,6 +728,8 @@ private:
// elements damage the mesh quality, allowing backtracking.
labelHashSet checkPolyMeshQuality(const pointField& pts) const;
label classifyBoundaryPoint(Cell_handle cit) const;
//- Index all of the the Delaunay cells and calculate their
//- dual points
void indexDualVertices
@ -736,7 +739,11 @@ private:
);
//- Re-index all of the Delaunay cells
void reindexDualVertices(const Map<label>& dualPtIndexMap);
void reindexDualVertices
(
const Map<label>& dualPtIndexMap,
labelList& boundaryPts
);
label createPatchInfo
(
@ -846,6 +853,8 @@ private:
const PtrList<dictionary>& patchDicts
) const;
void writePointPairs(const fileName& fName) const;
//- Disallow default bitwise copy construct
conformalVoronoiMesh(const conformalVoronoiMesh&);
@ -992,7 +1001,7 @@ public:
const wordList& patchNames,
const PtrList<dictionary>& patchDicts,
const pointField& cellCentres,
const PackedBoolList& boundaryFacesToRemove
PackedBoolList& boundaryFacesToRemove
) const;
//- Calculate and write a field of the target cell size,

663
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
View File

@ -30,487 +30,11 @@ License
#include "indexedCellChecks.H"
#include "OBJstream.H"
#include "indexedCellOps.H"
#include "ListOps.H"
#include "DelaunayMeshTools.H"
#include "CGAL/Exact_predicates_exact_constructions_kernel.h"
#include "CGAL/Gmpq.h"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void Foam::conformalVoronoiMesh::checkCells()
{
List<List<FixedList<Foam::point, 4> > > cellListList(Pstream::nProcs());
List<FixedList<Foam::point, 4> > cells(number_of_finite_cells());
globalIndex gIndex(number_of_vertices());
label count = 0;
for
(
Delaunay::Finite_cells_iterator cit = finite_cells_begin();
cit != finite_cells_end();
++cit
)
{
if (tetrahedron(cit).volume() == 0)
{
Pout<< "ZERO VOLUME TET" << endl;
Pout<< cit->info();
Pout<< cit->dual();
}
if (cit->hasFarPoint())
{
continue;
}
List<labelPair> cellVerticesPair(4);
List<Foam::point> cellVertices(4);
for (label vI = 0; vI < 4; ++vI)
{
cellVerticesPair[vI] = labelPair
(
cit->vertex(vI)->procIndex(),
cit->vertex(vI)->index()
);
cellVertices[vI] = topoint(cit->vertex(vI)->point());
}
List<Foam::point> cellVerticesOld(cellVertices);
labelList oldToNew;
sortedOrder(cellVerticesPair, oldToNew);
oldToNew = invert(oldToNew.size(), oldToNew);
inplaceReorder(oldToNew, cellVerticesPair);
inplaceReorder(oldToNew, cellVertices);
// FixedList<label, 4> globalTetCell
// (
// cit->globallyOrderedCellVertices(gIndex)
// );
//
// FixedList<Point, 4> cellVertices(Point(0,0,0));
//
// forAll(globalTetCell, gvI)
// {
// label gI = globalTetCell[gvI];
//
// cellVertices[gvI] = cit->vertex(gI)->point();
// }
// if (cit->hasFarPoint())
// {
// continue;
// }
for (label i = 0; i < 4; ++i)
{
//cells[count][i] = topoint(cit->vertex(i)->point());
cells[count][i] = cellVertices[i];
}
count++;
}
cells.setSize(count);
cellListList[Pstream::myProcNo()] = cells;
Pstream::gatherList(cellListList);
if (Pstream::master())
{
Info<< "Checking on master processor the cells of each " << nl
<< "processor point list against the master cell list." << nl
<< "There are " << cellListList.size() << " processors" << nl
<< "The size of each processor's cell list is:" << endl;
forAll(cellListList, cfI)
{
Info<< " Proc " << cfI << " has " << cellListList[cfI].size()
<< " cells" << endl;
}
label nMatches = 0, nMatchFoundDiffOrder = 0;
forAll(cellListList[0], cmI)
{
const FixedList<Foam::point, 4>& masterCell = cellListList[0][cmI];
bool matchFound = false;
bool matchFoundDiffOrder = false;
forAll(cellListList, cpI)
{
if (cpI == 0)
{
continue;
}
forAll(cellListList[cpI], csI)
{
const FixedList<Foam::point, 4>& slaveCell
= cellListList[cpI][csI];
if (masterCell == slaveCell)
{
matchFound = true;
break;
}
else
{
label samePt = 0;
forAll(masterCell, mI)
{
const Foam::point& mPt = masterCell[mI];
forAll(slaveCell, sI)
{
const Foam::point& sPt = slaveCell[sI];
if (mPt == sPt)
{
samePt++;
}
}
}
if (samePt == 4)
{
matchFoundDiffOrder = true;
Pout<< masterCell << nl << slaveCell << endl;
break;
}
}
}
}
if (matchFound)
{
nMatches++;
}
if (matchFoundDiffOrder)
{
nMatchFoundDiffOrder++;
}
}
Info<< "Found " << nMatches << " matching cells and "
<< nMatchFoundDiffOrder << " matching cells with different "
<< "vertex ordering"<< endl;
}
}
//void Foam::conformalVoronoiMesh::checkDuals()
//{
// List<List<Point> > pointFieldList(Pstream::nProcs());
//
// List<Point> duals(number_of_finite_cells());
//
// typedef CGAL::Exact_predicates_exact_constructions_kernel EK2;
// typedef CGAL::Regular_triangulation_euclidean_traits_3<EK2> EK;
// typedef CGAL::Cartesian_converter<baseK::Kernel, EK2> To_exact;
// typedef CGAL::Cartesian_converter<EK2, baseK::Kernel> Back_from_exact;
//
//// PackedBoolList bPoints(number_of_finite_cells());
//
//// indexDualVertices(duals, bPoints);
//
// label count = 0;//duals.size();
//
// duals.setSize(number_of_finite_cells());
//
// globalIndex gIndex(number_of_vertices());
//
// for
// (
// Delaunay::Finite_cells_iterator cit = finite_cells_begin();
// cit != finite_cells_end();
// ++cit
// )
// {
// if (cit->hasFarPoint())
// {
// continue;
// }
//
// duals[count++] = cit->circumcenter();
//
//// List<labelPair> cellVerticesPair(4);
//// List<Point> cellVertices(4);
////
//// for (label vI = 0; vI < 4; ++vI)
//// {
//// cellVerticesPair[vI] = labelPair
//// (
//// cit->vertex(vI)->procIndex(),
//// cit->vertex(vI)->index()
//// );
//// cellVertices[vI] = cit->vertex(vI)->point();
//// }
////
//// labelList oldToNew;
//// sortedOrder(cellVerticesPair, oldToNew);
//// oldToNew = invert(oldToNew.size(), oldToNew);
//// inplaceReorder(oldToNew, cellVerticesPair);
//// inplaceReorder(oldToNew, cellVertices);
////
//// duals[count++] = CGAL::circumcenter
//// (
//// cellVertices[0],
//// cellVertices[1],
//// cellVertices[2],
//// cellVertices[3]
//// );
//
//// To_exact to_exact;
//// Back_from_exact back_from_exact;
//// EK::Construct_circumcenter_3 exact_circumcenter =
//// EK().construct_circumcenter_3_object();
////
//// duals[count++] = topoint
//// (
//// back_from_exact
//// (
//// exact_circumcenter
//// (
//// to_exact(cit->vertex(0)->point()),
//// to_exact(cit->vertex(1)->point()),
//// to_exact(cit->vertex(2)->point()),
//// to_exact(cit->vertex(3)->point())
//// )
//// )
//// );
// }
//
// Pout<< "Duals Calculated " << count << endl;
//
// duals.setSize(count);
//
// pointFieldList[Pstream::myProcNo()] = duals;
//
// Pstream::gatherList(pointFieldList);
//
// if (Pstream::master())
// {
// Info<< "Checking on master processor the dual locations of each" << nl
// << " processor point list against the master dual list." << nl
// << "There are " << pointFieldList.size() << " processors" << nl
// << "The size of each processor's dual list is:" << endl;
//
// forAll(pointFieldList, pfI)
// {
// Info<< " Proc " << pfI << " has " << pointFieldList[pfI].size()
// << " duals" << endl;
// }
//
// label nNonMatches = 0;
// label nNearMatches = 0;
// label nExactMatches = 0;
//
// forAll(pointFieldList[0], pI)
// {
// const Point& masterPoint = pointFieldList[0][pI];
//
// bool foundMatch = false;
// bool foundNearMatch = false;
//
// scalar minCloseness = GREAT;
// Point closestPoint(0, 0, 0);
//
// forAll(pointFieldList, pfI)
// {
// if (pfI == 0)
// {
// continue;
// }
//
//// label pfI = 1;
//
// forAll(pointFieldList[pfI], pISlave)
// {
// const Point& slavePoint
// = pointFieldList[pfI][pISlave];
//
// if (masterPoint == slavePoint)
// {
// foundMatch = true;
// break;
// }
//
// const scalar closeness = mag
// (
// topoint(masterPoint) - topoint(slavePoint)
// );
//
// if (closeness < 1e-12)
// {
// foundNearMatch = true;
// }
// else
// {
// if (closeness < minCloseness)
// {
// minCloseness = closeness;
// closestPoint = slavePoint;
// }
// }
// }
//
// if (!foundMatch)
// {
// if (foundNearMatch)
// {
// CGAL::Gmpq x(CGAL::to_double(masterPoint.x()));
// CGAL::Gmpq y(CGAL::to_double(masterPoint.y()));
// CGAL::Gmpq z(CGAL::to_double(masterPoint.z()));
//
// std::cout<< "master = " << x << " " << y << " " << z
// << std::endl;
//
// CGAL::Gmpq xs(CGAL::to_double(closestPoint.x()));
// CGAL::Gmpq ys(CGAL::to_double(closestPoint.y()));
// CGAL::Gmpq zs(CGAL::to_double(closestPoint.z()));
// std::cout<< "slave = " << xs << " " << ys << " "
// << zs
// << std::endl;
//
// nNearMatches++;
// }
// else
// {
// nNonMatches++;
// Info<< "Closest point to " << masterPoint << " is "
// << closestPoint << nl
// << " Separation is " << minCloseness << endl;
//
// CGAL::Gmpq x(CGAL::to_double(masterPoint.x()));
// CGAL::Gmpq y(CGAL::to_double(masterPoint.y()));
// CGAL::Gmpq z(CGAL::to_double(masterPoint.z()));
//
// std::cout<< "master = " << x << " " << y << " " << z
// << std::endl;
//
// CGAL::Gmpq xs(CGAL::to_double(closestPoint.x()));
// CGAL::Gmpq ys(CGAL::to_double(closestPoint.y()));
// CGAL::Gmpq zs(CGAL::to_double(closestPoint.z()));
// std::cout<< "slave = " << xs << " " << ys << " "
// << zs
// << std::endl;
// }
// }
// else
// {
// nExactMatches++;
// }
// }
// }
//
// Info<< "Found " << nNonMatches << " non-matching duals" << nl
// << " and " << nNearMatches << " near matches"
// << " and " << nExactMatches << " exact matches" << endl;
// }
//}
void Foam::conformalVoronoiMesh::checkVertices()
{
List<pointField> pointFieldList(Pstream::nProcs());
pointField points(number_of_vertices());
labelPairHashSet duplicateVertices;
label count = 0;
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
++vit
)
{
if (duplicateVertices.found(labelPair(vit->procIndex(), vit->index())))
{
Pout<< "DUPLICATE " << vit->procIndex() << vit->index() << endl;
}
else
{
duplicateVertices.insert(labelPair(vit->procIndex(), vit->index()));
}
points[count++] = topoint(vit->point());
}
pointFieldList[Pstream::myProcNo()] = points;
Pstream::gatherList(pointFieldList);
OFstream str("missingPoints.obj");
if (Pstream::master())
{
Info<< "Checking on master processor the point locations of each " << nl
<< "processor point list against the master point list." << nl
<< "There are " << pointFieldList.size() << " processors" << nl
<< "The size of each processor's point list is:" << endl;
forAll(pointFieldList, pfI)
{
Info<< " Proc " << pfI << " has " << pointFieldList[pfI].size()
<< " points" << endl;
}
label nNonMatches = 0;
forAll(pointFieldList[0], pI)
{
const Foam::point& masterPoint = pointFieldList[0][pI];
forAll(pointFieldList, pfI)
{
if (pI == 0)
{
continue;
}
bool foundMatch = false;
forAll(pointFieldList[pfI], pISlave)
{
const Foam::point& slavePoint
= pointFieldList[pfI][pISlave];
if (masterPoint == slavePoint)
{
foundMatch = true;
break;
}
}
if (!foundMatch)
{
Info<< " Proc " << pfI << " Master != Slave -> "
<< masterPoint << endl;
meshTools::writeOBJ(str, masterPoint);
nNonMatches++;
}
}
}
Info<< "Found a total of " << nNonMatches << " non-matching points"
<< endl;
}
}
void Foam::conformalVoronoiMesh::calcDualMesh
(
pointField& points,
@ -528,53 +52,6 @@ void Foam::conformalVoronoiMesh::calcDualMesh
{
timeCheck("Start calcDualMesh");
// if (debug)
// {
// Pout<< nl << "Perfoming some checks . . ." << nl << nl
// << "Total number of vertices = " << number_of_vertices() << nl
// << "Total number of cells = " << number_of_finite_cells()
// << endl;
//
// checkVertices();
// checkCells();
// checkDuals();
//
// Info<< nl << "Finished checks" << nl << endl;
// }
// OFstream str("attachedToFeature.obj");
// label offset = 0;
//
// for
// (
// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
// vit != finite_vertices_end();
// ++vit
// )
// {
// if (vit->featurePoint())
// {
// std::list<Cell_handle> adjacentCells;
//
// finite_incident_cells(vit, std::back_inserter(adjacentCells));
//
// for
// (
// std::list<Cell_handle>::iterator acit = adjacentCells.begin();
// acit != adjacentCells.end();
// ++acit
// )
// {
// if ((*acit)->real())
// {
// drawDelaunayCell(str, (*acit), offset);
// offset++;
//// meshTools::writeOBJ(str, topoint((*acit)->dual()));
// }
// }
// }
// }
setVertexSizeAndAlignment();
timeCheck("After setVertexSizeAndAlignment");
@ -585,7 +62,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
Info<< nl << "Merging identical points" << endl;
// There is no guarantee that a merge of close points is no-risk
mergeIdenticalDualVertices(points);
mergeIdenticalDualVertices(points, boundaryPts);
}
// Final dual face and owner neighbour construction
@ -820,7 +297,8 @@ void Foam::conformalVoronoiMesh::calcTetMesh
void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
(
const pointField& pts
const pointField& pts,
labelList& boundaryPts
)
{
// Assess close points to be merged
@ -838,7 +316,7 @@ void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
dualPtIndexMap
);
reindexDualVertices(dualPtIndexMap);
reindexDualVertices(dualPtIndexMap, boundaryPts);
reduce(nPtsMerged, sumOp<label>());
@ -1245,7 +723,6 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
false
);
//createCellCentres(cellCentres);
cellCentres = DelaunayMeshTools::allPoints(*this);
labelList cellToDelaunayVertex(removeUnusedCells(owner, neighbour));
@ -1327,24 +804,6 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
pMesh.addPatches(patches);
// Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
// forAll(patches, p)
// {
// patches[p] = new polyPatch
// (
// patchNames[p],
// patchSizes[p],
// patchStarts[p],
// p,
// pMesh.boundaryMesh()
// );
// }
// pMesh.addPatches(patches, false);
// pMesh.overrideCellCentres(cellCentres);
return meshPtr;
}
@ -1361,7 +820,7 @@ void Foam::conformalVoronoiMesh::checkCellSizing()
indexDualVertices(ptsField, boundaryPts);
// Merge close dual vertices.
mergeIdenticalDualVertices(ptsField);
mergeIdenticalDualVertices(ptsField, boundaryPts);
autoPtr<polyMesh> meshPtr = createPolyMeshFromPoints(ptsField);
const polyMesh& pMesh = meshPtr();
@ -1740,6 +1199,41 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
}
Foam::label Foam::conformalVoronoiMesh::classifyBoundaryPoint
(
Cell_handle cit
) const
{
if (cit->boundaryDualVertex())
{
if (cit->featurePointDualVertex())
{
return featurePoint;
}
else if (cit->featureEdgeDualVertex())
{
return featureEdge;
}
else
{
return surface;
}
}
else if (cit->baffleSurfaceDualVertex())
{
return surface;
}
else if (cit->baffleEdgeDualVertex())
{
return featureEdge;
}
else
{
return internal;
}
}
void Foam::conformalVoronoiMesh::indexDualVertices
(
pointField& pts,
@ -1989,42 +1483,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
// }
// }
if (cit->boundaryDualVertex())
{
if (cit->featurePointDualVertex())
{
boundaryPts[cit->cellIndex()] = featurePoint;
}
else if (cit->featureEdgeDualVertex())
{
boundaryPts[cit->cellIndex()] = featureEdge;
}
else
{
boundaryPts[cit->cellIndex()] = surface;
}
}
else if
(
cit->baffleBoundaryDualVertex()
)
{
boundaryPts[cit->cellIndex()] = surface;
}
else if
(
cit->vertex(0)->featureEdgePoint()
&& cit->vertex(1)->featureEdgePoint()
&& cit->vertex(2)->featureEdgePoint()
&& cit->vertex(3)->featureEdgePoint()
)
{
boundaryPts[cit->cellIndex()] = featureEdge;
}
else
{
boundaryPts[cit->cellIndex()] = internal;
}
boundaryPts[cit->cellIndex()] = classifyBoundaryPoint(cit);
}
else
{
@ -2040,7 +1499,8 @@ void Foam::conformalVoronoiMesh::indexDualVertices
void Foam::conformalVoronoiMesh::reindexDualVertices
(
const Map<label>& dualPtIndexMap
const Map<label>& dualPtIndexMap,
labelList& boundaryPts
)
{
for
@ -2053,6 +1513,12 @@ void Foam::conformalVoronoiMesh::reindexDualVertices
if (dualPtIndexMap.found(cit->cellIndex()))
{
cit->cellIndex() = dualPtIndexMap[cit->cellIndex()];
boundaryPts[cit->cellIndex()] =
max
(
boundaryPts[cit->cellIndex()],
boundaryPts[dualPtIndexMap[cit->cellIndex()]]
);
}
}
}
@ -2267,11 +1733,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
bool includeEmptyPatches
) const
{
const label defaultPatchIndex = createPatchInfo
(
patchNames,
patchDicts
);
const label defaultPatchIndex = createPatchInfo(patchNames, patchDicts);
const label nPatches = patchNames.size();
@ -2296,6 +1758,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
List<DynamicList<bool> > indirectPatchFace(nPatches, DynamicList<bool>(0));
faces.setSize(number_of_finite_edges());
owner.setSize(number_of_finite_edges());
neighbour.setSize(number_of_finite_edges());
@ -2765,7 +2228,12 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
// If the two vertices are a pair, then the patch face is
// a desired one.
if (!isPointPair(vA, vB))
if
(
vA->boundaryPoint() && vB->boundaryPoint()
&& !ptPairs_.isPointPair(vA, vB)
&& !ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
)
{
indirectPatchFace[patchIndex].append(true);
}
@ -2797,6 +2265,18 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
if (patchIndex != -1)
{
// if
// (
// vA->boundaryPoint() && vB->boundaryPoint()
// && !ptPairs_.isPointPair(vA, vB)
// )
// {
// indirectPatchFace[patchIndex].append(true);
// }
// else
// {
// indirectPatchFace[patchIndex].append(false);
// }
// patchFaces[patchIndex].append(newDualFace);
// patchOwners[patchIndex].append(own);
// indirectPatchFace[patchIndex].append(false);
@ -2813,8 +2293,6 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
// {
// patchPPSlaves[patchIndex].append(vA->index());
// }
// baffleFaces[dualFaceI] = patchIndex;
}
// else
{
@ -2989,7 +2467,6 @@ void Foam::conformalVoronoiMesh::sortProcPatches
faceList& faces = patchFaces[patchI];
labelList& owner = patchOwners[patchI];
DynamicList<label>& slaves = patchPointPairSlaves[patchI];
DynamicList<Pair<labelPair> >& sortingIndices
= patchSortingIndices[patchI];

456
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
View File

@ -62,11 +62,13 @@ void Foam::conformalVoronoiMesh::conformToSurface()
if (Pstream::parRun())
{
sync(decomposition_().procBounds());
sync(decomposition().procBounds());
}
}
else
{
ptPairs_.clear();
// Rebuild, insert and store new surface conformation
buildSurfaceConformation();
@ -74,7 +76,7 @@ void Foam::conformalVoronoiMesh::conformToSurface()
{
if (Pstream::parRun())
{
sync(decomposition_().procBounds());
sync(decomposition().procBounds());
}
}
@ -358,10 +360,13 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
}
DynamicList<Vb> pts(2*surfaceHits.size() + 3*featureEdgeHits.size());
insertSurfacePointPairs
(
surfaceHits,
"surfaceConformationLocations_initial.obj"
"surfaceConformationLocations_initial.obj",
pts
);
// In parallel, synchronise the edge trees
@ -373,9 +378,21 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
insertEdgePointGroups
(
featureEdgeHits,
"edgeConformationLocations_initial.obj"
"edgeConformationLocations_initial.obj",
pts
);
pts.shrink();
Map<label> oldToNewIndices = insertPointPairs(pts, true, true);
// Re-index the point pairs
ptPairs_.reIndex(oldToNewIndices);
//writePointPairs("pointPairs_initial.obj");
// Remove location from surface/edge tree
timeCheck("After initial conformation");
initialTotalHits = nSurfHits + nFeatEdHits;
@ -550,215 +567,6 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
}
}
// for
// (
// Delaunay::Finite_edges_iterator eit = finite_edges_begin();
// eit != finite_edges_end();
// ++eit
// )
// {
// Cell_handle c = eit->first;
// Vertex_handle vA = c->vertex(eit->second);
// Vertex_handle vB = c->vertex(eit->third);
//
// if
// (
// vA->referred()
// || vB->referred()
// )
// {
// continue;
// }
//
// if
// (
// (vA->internalPoint() && vA->referred())
// || (vB->internalPoint() && vB->referred())
// )
// {
// continue;
// }
//
// if
// (
// (vA->internalPoint() && vB->externalBoundaryPoint())
// || (vB->internalPoint() && vA->externalBoundaryPoint())
// || (vA->internalBoundaryPoint() && vB->internalBoundaryPoint())
// )
// {
// pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
// pointIndexHit surfHit;
// label hitSurface;
//
// geometryToConformTo_.findSurfaceNearest
// (
// (
// vA->internalPoint()
// ? topoint(vA->point())
// : topoint(vB->point())
// ),
// magSqr(topoint(vA->point()) - topoint(vB->point())),
// surfHit,
// hitSurface
// );
//
// if (surfHit.hit())
// {
// surfaceIntersections.append
// (
// pointIndexHitAndFeature(surfHit, hitSurface)
// );
//
// addSurfaceAndEdgeHits
// (
// (
// vA->internalPoint()
// ? topoint(vA->point())
// : topoint(vB->point())
// ),
// surfaceIntersections,
// surfacePtReplaceDistCoeffSqr,
// edgeSearchDistCoeffSqr,
// surfaceHits,
// featureEdgeHits,
// surfaceToTreeShape,
// edgeToTreeShape,
// surfacePtToEdgePtDist,
// false
// );
// }
// }
// }
for
(
Delaunay::Finite_cells_iterator cit = finite_cells_begin();
cit != finite_cells_end();
++cit
)
{
label nInternal = 0;
for (label vI = 0; vI < 4; vI++)
{
if (cit->vertex(vI)->internalPoint())
{
nInternal++;
}
}
if (nInternal == 1 && cit->hasBoundaryPoint())
//if (cit->boundaryDualVertex() && !cit->hasReferredPoint())
{
const Foam::point& pt = cit->dual();
const Foam::point cellCentre =
topoint
(
CGAL::centroid
(
CGAL::Tetrahedron_3<baseK>
(
cit->vertex(0)->point(),
cit->vertex(1)->point(),
cit->vertex(2)->point(),
cit->vertex(3)->point()
)
)
);
pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
pointIndexHit surfHit;
label hitSurface;
geometryToConformTo_.findSurfaceNearestIntersection
(
cellCentre,
pt,
surfHit,
hitSurface
);
if (surfHit.hit())
{
surfaceIntersections.append
(
pointIndexHitAndFeature(surfHit, hitSurface)
);
addSurfaceAndEdgeHits
(
pt,
surfaceIntersections,
surfacePtReplaceDistCoeffSqr,
edgeSearchDistCoeffSqr,
surfaceHits,
featureEdgeHits,
surfaceToTreeShape,
edgeToTreeShape,
surfacePtToEdgePtDist,
false
);
}
}
}
// for
// (
// Delaunay::Finite_cells_iterator cit = finite_cells_begin();
// cit != finite_cells_end();
// ++cit
// )
// {
// if (cit->boundaryDualVertex() && !cit->hasReferredPoint())
// {
// const Foam::point& pt = cit->dual();
//
// pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
// pointIndexHit surfHit;
// label hitSurface;
//
// geometryToConformTo_.findSurfaceNearest
// (
// pt,
// sqr(0.1*targetCellSize(pt)),
// surfHit,
// hitSurface
// );
//
// if (!surfHit.hit())
// {
// geometryToConformTo_.findSurfaceNearest
// (
// pt,
// GREAT,
// surfHit,
// hitSurface
// );
//
// if (surfHit.hit())
// {
// surfaceIntersections.append
// (
// pointIndexHitAndFeature(surfHit, hitSurface)
// );
//
// addSurfaceAndEdgeHits
// (
// pt,
// surfaceIntersections,
// surfacePtReplaceDistCoeffSqr,
// edgeSearchDistCoeffSqr,
// surfaceHits,
// featureEdgeHits,
// surfaceToTreeShape,
// edgeToTreeShape,
// false
// );
// }
// }
// }
// }
label nVerts = number_of_vertices();
label nSurfHits = surfaceHits.size();
label nFeatEdHits = featureEdgeHits.size();
@ -789,10 +597,16 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
}
DynamicList<Vb> pts
(
2*surfaceHits.size() + 3*featureEdgeHits.size()
);
insertSurfacePointPairs
(
surfaceHits,
"surfaceConformationLocations_" + name(iterationNo) + ".obj"
"surfaceConformationLocations_" + name(iterationNo) + ".obj",
pts
);
// In parallel, synchronise the edge trees
@ -805,9 +619,19 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
insertEdgePointGroups
(
featureEdgeHits,
"edgeConformationLocations_" + name(iterationNo) + ".obj"
"edgeConformationLocations_" + name(iterationNo) + ".obj",
pts
);
pts.shrink();
Map<label> oldToNewIndices = insertPointPairs(pts, true, true);
// Reindex the point pairs
ptPairs_.reIndex(oldToNewIndices);
//writePointPairs("pointPairs_" + name(iterationNo) + ".obj");
if (Pstream::parRun())
{
sync
@ -1155,8 +979,6 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
|| is_infinite(fit->first->neighbor(fit->second))
|| !fit->first->hasInternalPoint()
|| !fit->first->neighbor(fit->second)->hasInternalPoint()
// || fit->first->hasFarPoint()
// || fit->first->neighbor(fit->second)->hasFarPoint()
)
{
continue;
@ -1209,14 +1031,13 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
const scalar d = p.normalIntersect(r);
const Foam::point newPoint = vertex + d*n;
Foam::point newPoint = vertex + d*n;
pointIndexHitAndFeature info;
geometryToConformTo_.findSurfaceNearest
(
newPoint,
surfaceSearchDistanceSqr(newPoint),
4.0*magSqr(newPoint - vertex),
info.first(),
info.second()
);
@ -1359,7 +1180,10 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
if
(
is_infinite(c1) || is_infinite(c2)
|| (!c1->hasInternalPoint() && !c2->hasInternalPoint())
|| (
!c1->internalOrBoundaryDualVertex()
|| !c2->internalOrBoundaryDualVertex()
)
|| !c1->real() || !c2->real()
)
{
@ -1369,19 +1193,24 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
// Foam::point endPt = 0.5*(c1->dual() + c2->dual());
Foam::point endPt = c1->dual();
if
(
magSqr(vert - c1->dual())
< magSqr(vert - c2->dual())
)
if (magSqr(vert - c1->dual()) < magSqr(vert - c2->dual()))
{
endPt = c2->dual();
}
if
(
magSqr(vert - endPt)
> magSqr(geometryToConformTo().globalBounds().mag())
)
{
continue;
}
pointIndexHit surfHit;
label hitSurface;
geometryToConformTo_.findSurfaceAnyIntersection
geometryToConformTo_.findSurfaceNearestIntersection
(
vert,
endPt,
@ -1401,15 +1230,46 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
const vector& n = norm[0];
const scalar normalProtrusionDistance =
(endPt - surfHit.hitPoint()) & n;
const scalar normalProtrusionDistance
(
(endPt - surfHit.hitPoint()) & n
);
if (normalProtrusionDistance > maxProtrusionDistance)
{
surfHitLargest = surfHit;
hitSurfaceLargest = hitSurface;
const plane p(surfHit.hitPoint(), n);
maxProtrusionDistance = normalProtrusionDistance;
const plane::ray r(endPt, -n);
const scalar d = p.normalIntersect(r);
Foam::point newPoint = endPt - d*n;
pointIndexHitAndFeature info;
geometryToConformTo_.findSurfaceNearest
(
newPoint,
4.0*magSqr(newPoint - endPt),
info.first(),
info.second()
);
if (info.first().hit())
{
if
(
surfaceLocationConformsToInside
(
pointIndexHitAndFeature(info.first(), info.second())
)
)
{
surfHitLargest = info.first();
hitSurfaceLargest = info.second();
maxProtrusionDistance = normalProtrusionDistance;
}
}
}
}
}
@ -1465,7 +1325,10 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
if
(
is_infinite(c1) || is_infinite(c2)
|| (!c1->hasInternalPoint() && !c2->hasInternalPoint())
|| (
!c1->internalOrBoundaryDualVertex()
|| !c2->internalOrBoundaryDualVertex()
)
|| !c1->real() || !c2->real()
)
{
@ -1475,19 +1338,24 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
// Foam::point endPt = 0.5*(c1->dual() + c2->dual());
Foam::point endPt = c1->dual();
if
(
magSqr(vert - c1->dual())
< magSqr(vert - c2->dual())
)
if (magSqr(vert - c1->dual()) < magSqr(vert - c2->dual()))
{
endPt = c2->dual();
}
if
(
magSqr(vert - endPt)
> magSqr(geometryToConformTo().globalBounds().mag())
)
{
continue;
}
pointIndexHit surfHit;
label hitSurface;
geometryToConformTo_.findSurfaceAnyIntersection
geometryToConformTo_.findSurfaceNearestIntersection
(
vert,
endPt,
@ -1507,19 +1375,46 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
const vector& n = norm[0];
scalar normalIncursionDistance = (endPt - surfHit.hitPoint()) & n;
scalar normalIncursionDistance
(
(endPt - surfHit.hitPoint()) & n
);
if (normalIncursionDistance < minIncursionDistance)
{
surfHitLargest = surfHit;
hitSurfaceLargest = hitSurface;
const plane p(surfHit.hitPoint(), n);
minIncursionDistance = normalIncursionDistance;
const plane::ray r(endPt, n);
// Info<< nl << "# Incursion: " << endl;
// meshTools::writeOBJ(Info, vert);
// meshTools::writeOBJ(Info, edgeMid);
// Info<< "l Na Nb" << endl;
const scalar d = p.normalIntersect(r);
Foam::point newPoint = endPt + d*n;
pointIndexHitAndFeature info;
geometryToConformTo_.findSurfaceNearest
(
newPoint,
4.0*magSqr(newPoint - endPt),
info.first(),
info.second()
);
if (info.first().hit())
{
if
(
surfaceLocationConformsToInside
(
pointIndexHitAndFeature(info.first(), info.second())
)
)
{
surfHitLargest = info.first();
hitSurfaceLargest = info.second();
minIncursionDistance = normalIncursionDistance;
}
}
}
}
}
@ -2162,9 +2057,11 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
bool isNearFeaturePt = nearFeaturePt(surfPt);
bool isNearFeatureEdge = surfacePtNearFeatureEdge(surfPt);
bool isNearSurfacePt = nearSurfacePoint(surfHitI);
if (isNearFeaturePt || isNearSurfacePt)
if (isNearFeaturePt || isNearSurfacePt || isNearFeatureEdge)
{
keepSurfacePoint = false;
}
@ -2219,6 +2116,9 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
// feature edge, give control to edge control points
// instead, this will prevent "pits" forming.
// Allow if different surfaces
keepSurfacePoint = false;
}
@ -2298,11 +2198,13 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
surfaceHits.append(surfHitI);
appendToSurfacePtTree(surfPt);
surfaceToTreeShape.append(existingSurfacePtLocations_.size() - 1);
// addedPoints.write(surfPt);
}
else
{
// removedPoints.write(surfPt);
}
// else
// {
// surfaceToTreeShape.remove();
// }
}
}
@ -2338,7 +2240,9 @@ void Foam::conformalVoronoiMesh::storeSurfaceConformation()
Vb
(
vit->point(),
vit->type()
vit->index(),
vit->type(),
Pstream::myProcNo()
)
);
}
@ -2362,24 +2266,40 @@ void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
{
Info<< nl << "Reinserting stored surface conformation" << endl;
const label preReinsertionSize(number_of_vertices());
Map<label> oldToNewIndices =
insertPointPairs(surfaceConformationVertices_, true, true);
// It is assumed that the stored surface conformation is on the correct
// processor and does not need distributed
rangeInsertWithInfo
ptPairs_.reIndex(oldToNewIndices);
PackedBoolList selectedElems(surfaceConformationVertices_.size(), true);
forAll(surfaceConformationVertices_, vI)
{
Vb& v = surfaceConformationVertices_[vI];
label& vIndex = v.index();
Map<label>::const_iterator iter = oldToNewIndices.find(vIndex);
if (iter != oldToNewIndices.end())
{
const label newIndex = iter();
if (newIndex != -1)
{
vIndex = newIndex;
}
else
{
selectedElems[vI] = false;
}
}
}
inplaceSubset<PackedBoolList, List<Vb> >
(
surfaceConformationVertices_.begin(),
surfaceConformationVertices_.end(),
true
selectedElems,
surfaceConformationVertices_
);
const label nInserted = label(number_of_vertices()) - preReinsertionSize;
const label nFailed = surfaceConformationVertices_.size() - nInserted;
Info<< " " << returnReduce(nInserted, sumOp<label>())
<< " points reinserted, failed to insert "
<< returnReduce(nFailed, sumOp<label>())
<< endl;
}

429
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
View File

@ -29,6 +29,7 @@ License
#include "tetrahedron.H"
#include "const_circulator.H"
#include "DelaunayMeshTools.H"
#include "OBJstream.H"
using namespace Foam::vectorTools;
@ -225,21 +226,48 @@ void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
vector masterPtVec(normalDir + nextNormalDir);
masterPtVec /= mag(masterPtVec) + SMALL;
if (((normalDir ^ nextNormalDir) & edDir) < SMALL)
if
(
((normalDir ^ nextNormalDir) & edDir) < SMALL
|| mag(masterPtVec) < SMALL
)
{
// Info<< " IGNORE REGION" << endl;
addedMasterPreviously = false;
if (circ.size() == 2 && mag((normal & nextNormal) - 1) < SMALL)
if
(
circ.size() == 2
&& mag((normal & nextNormal) - 1) < SMALL
)
{
// Add an extra point
const vector n = 0.5*(normal + nextNormal);
vector s = ppDist*(edDir ^ n);
const vector s = ppDist*(edDir ^ n);
createPointPair(ppDist, edgePt + s, n, pts);
createPointPair(ppDist, edgePt - s, n, pts);
if (volType == extendedFeatureEdgeMesh::BOTH)
{
createBafflePointPair(ppDist, edgePt + s, n, true, pts);
createBafflePointPair(ppDist, edgePt - s, n, true, pts);
}
else
{
WarningIn
(
"Foam::conformalVoronoiMesh::"
"createEdgePointGroupByCirculating"
"("
" const extendedFeatureEdgeMesh&,"
" const pointIndexHit&,"
" DynamicList<Vb>&"
")"
) << "Failed to insert flat/open edge as volType is "
<< extendedFeatureEdgeMesh::sideVolumeTypeNames_
[
volType
]
<< endl;
}
break;
}
@ -247,34 +275,6 @@ void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
continue;
}
if (mag(masterPtVec) < SMALL)
{
if (circ.size() == 2)
{
// Add an extra point
// Average normal to remove any bias to one side, although as
// it is a flat edge, the normals should be essentially the same
const vector n = 0.5*(normal + nextNormal);
// Direction along the surface to the control point, sense of
// direction not important, as +s and -s can be used because it
// is a flat edge
vector s = ppDist*(edDir ^ n);
createPointPair(ppDist, edgePt + s, n, pts);
createPointPair(ppDist, edgePt - s, n, pts);
break;
}
else
{
// Info<< " IGNORE REGION" << endl;
addedMasterPreviously = false;
continue;
}
}
const Foam::point masterPt = edgePt + ppDist*masterPtVec;
// Check that region is inside or outside
@ -490,11 +490,19 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
feMesh.normalVolumeTypes();
// As this is an external edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
const extendedFeatureEdgeMesh::sideVolumeType& volTypeA =
normalVolumeTypes[edNormalIs[0]];
const extendedFeatureEdgeMesh::sideVolumeType& volTypeB =
normalVolumeTypes[edNormalIs[1]];
if (areParallel(nA, nB))
{
// The normals are nearly parallel, so this is too sharp a feature to
@ -524,13 +532,6 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
// 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
if (!geometryToConformTo_.inside(refPt))
@ -559,7 +560,11 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
(
reflectedA,
vertexCount() + pts.size(),
Vb::vtExternalFeatureEdge,
(
volTypeA == extendedFeatureEdgeMesh::BOTH
? Vb::vtInternalFeatureEdge
: Vb::vtExternalFeatureEdge
),
Pstream::myProcNo()
)
);
@ -571,10 +576,26 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
(
reflectedB,
vertexCount() + pts.size(),
Vb::vtExternalFeatureEdge,
(
volTypeB == extendedFeatureEdgeMesh::BOTH
? Vb::vtInternalFeatureEdge
: Vb::vtExternalFeatureEdge
),
Pstream::myProcNo()
)
);
ptPairs_.addPointPair
(
pts[pts.size() - 3].index(),
pts[pts.size() - 1].index()
);
ptPairs_.addPointPair
(
pts[pts.size() - 3].index(),
pts[pts.size() - 2].index()
);
}
@ -591,11 +612,19 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
feMesh.normalVolumeTypes();
// As this is an external edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
const extendedFeatureEdgeMesh::sideVolumeType& volTypeA =
normalVolumeTypes[edNormalIs[0]];
// const extendedFeatureEdgeMesh::sideVolumeType& volTypeB =
// normalVolumeTypes[edNormalIs[1]];
if (areParallel(nA, nB))
{
// The normals are nearly parallel, so this is too sharp a feature to
@ -705,14 +734,30 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
(
reflMasterPt,
vertexCount() + pts.size(),
Vb::vtExternalFeatureEdge,
(
volTypeA == extendedFeatureEdgeMesh::BOTH
? Vb::vtInternalFeatureEdge
: Vb::vtExternalFeatureEdge
),
Pstream::myProcNo()
)
);
ptPairs_.addPointPair
(
pts[pts.size() - 2].index(),
pts[pts.size() - 1].index()
);
ptPairs_.addPointPair
(
pts[pts.size() - 3].index(),
pts[pts.size() - 1].index()
);
if (nAddPoints == 1)
{
// One additinal point is the reflection of the slave point,
// One additional point is the reflection of the slave point,
// i.e. the original reference point
pts.append
(
@ -767,6 +812,8 @@ void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
feMesh.normalVolumeTypes();
// As this is a flat edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
@ -781,8 +828,21 @@ void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
// is a flat edge
vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
createPointPair(ppDist, edgePt + s, n, pts);
createPointPair(ppDist, edgePt - s, n, pts);
if (normalVolumeTypes[edNormalIs[0]] == extendedFeatureEdgeMesh::OUTSIDE)
{
createPointPair(ppDist, edgePt + s, -n, true, pts);
createPointPair(ppDist, edgePt - s, -n, true, pts);
}
else if (normalVolumeTypes[edNormalIs[0]] == extendedFeatureEdgeMesh::BOTH)
{
createBafflePointPair(ppDist, edgePt + s, n, true, pts);
createBafflePointPair(ppDist, edgePt - s, n, true, pts);
}
else
{
createPointPair(ppDist, edgePt + s, n, true, pts);
createPointPair(ppDist, edgePt - s, n, true, pts);
}
}
@ -814,16 +874,23 @@ void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
plane facePlane(edgePt, n);
Foam::point pt1 = edgePt + s + ppDist*n;
Foam::point pt2 = edgePt - s + ppDist*n;
const label initialPtsSize = pts.size();
Foam::point pt3 = facePlane.mirror(pt1);
Foam::point pt4 = facePlane.mirror(pt2);
if
(
!geometryToConformTo_.inside(edgePt)
)
{
return;
}
pts.append(Vb(pt1, Vb::vtInternalSurface));
pts.append(Vb(pt2, Vb::vtInternalSurface));
pts.append(Vb(pt3, Vb::vtInternalSurface));
pts.append(Vb(pt4, Vb::vtInternalSurface));
createBafflePointPair(ppDist, edgePt - s, n, true, pts);
createBafflePointPair(ppDist, edgePt + s, n, false, pts);
for (label ptI = initialPtsSize; ptI < pts.size(); ++ptI)
{
pts[ptI].type() = Vb::vtInternalFeatureEdge;
}
}
else
{
@ -846,26 +913,245 @@ void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
const scalar ppDist = pointPairDistance(edgePt);
const vector edDir = feMesh.edgeDirections()[edHit.index()];
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
const labelList& normalDirs = feMesh.normalDirections()[edHit.index()];
Info<< edNormalIs.size() << endl;
const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
feMesh.normalVolumeTypes();
// As this is a flat edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
labelList nNormalTypes(4, 0);
// 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);
forAll(edNormalIs, edgeNormalI)
{
const extendedFeatureEdgeMesh::sideVolumeType sType =
normalVolumeTypes[edNormalIs[edgeNormalI]];
// 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);
nNormalTypes[sType]++;
}
createPointPair(ppDist, edgePt + s, n, pts);
createPointPair(ppDist, edgePt - s, n, pts);
if (nNormalTypes[extendedFeatureEdgeMesh::BOTH] == 4)
{
label masterEdgeNormalIndex = -1;
forAll(edNormalIs, edgeNormalI)
{
const extendedFeatureEdgeMesh::sideVolumeType sType =
normalVolumeTypes[edNormalIs[edgeNormalI]];
if (sType == extendedFeatureEdgeMesh::BOTH)
{
masterEdgeNormalIndex = edgeNormalI;
break;
}
}
const vector& n = feNormals[edNormalIs[masterEdgeNormalIndex]];
label nDir = normalDirs[masterEdgeNormalIndex];
vector normalDir =
(feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
normalDir *= nDir/mag(normalDir);
Foam::point pt1 = edgePt + ppDist*normalDir + ppDist*n;
Foam::point pt2 = edgePt + ppDist*normalDir - ppDist*n;
plane plane3(edgePt, normalDir);
Foam::point pt3 = plane3.mirror(pt1);
Foam::point pt4 = plane3.mirror(pt2);
pts.append
(
Vb
(
pt1,
vertexCount() + pts.size(),
Vb::vtInternalSurface,
Pstream::myProcNo()
)
);
pts.append
(
Vb
(
pt2,
vertexCount() + pts.size(),
Vb::vtInternalSurface,
Pstream::myProcNo()
)
);
ptPairs_.addPointPair
(
pts[pts.size() - 2].index(), // external 0 -> slave
pts[pts.size() - 1].index()
);
pts.append
(
Vb
(
pt3,
vertexCount() + pts.size(),
Vb::vtInternalSurface,
Pstream::myProcNo()
)
);
ptPairs_.addPointPair
(
pts[pts.size() - 3].index(), // external 0 -> slave
pts[pts.size() - 1].index()
);
pts.append
(
Vb
(
pt4,
vertexCount() + pts.size(),
Vb::vtInternalSurface,
Pstream::myProcNo()
)
);
ptPairs_.addPointPair
(
pts[pts.size() - 3].index(), // external 0 -> slave
pts[pts.size() - 1].index()
);
ptPairs_.addPointPair
(
pts[pts.size() - 2].index(), // external 0 -> slave
pts[pts.size() - 1].index()
);
}
else if
(
nNormalTypes[extendedFeatureEdgeMesh::BOTH] == 1
&& nNormalTypes[extendedFeatureEdgeMesh::INSIDE] == 2
)
{
label masterEdgeNormalIndex = -1;
forAll(edNormalIs, edgeNormalI)
{
const extendedFeatureEdgeMesh::sideVolumeType sType =
normalVolumeTypes[edNormalIs[edgeNormalI]];
if (sType == extendedFeatureEdgeMesh::BOTH)
{
masterEdgeNormalIndex = edgeNormalI;
break;
}
}
const vector& n = feNormals[edNormalIs[masterEdgeNormalIndex]];
label nDir = normalDirs[masterEdgeNormalIndex];
vector normalDir =
(feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
normalDir *= nDir/mag(normalDir);
const label nextNormalI =
(masterEdgeNormalIndex + 1) % edNormalIs.size();
if ((normalDir & feNormals[edNormalIs[nextNormalI]]) > 0)
{
normalDir *= -1;
}
Foam::point pt1 = edgePt + ppDist*normalDir + ppDist*n;
Foam::point pt2 = edgePt + ppDist*normalDir - ppDist*n;
plane plane3(edgePt, normalDir);
Foam::point pt3 = plane3.mirror(pt1);
Foam::point pt4 = plane3.mirror(pt2);
pts.append
(
Vb
(
pt1,
vertexCount() + pts.size(),
Vb::vtInternalSurface,
Pstream::myProcNo()
)
);
pts.append
(
Vb
(
pt2,
vertexCount() + pts.size(),
Vb::vtInternalSurface,
Pstream::myProcNo()
)
);
ptPairs_.addPointPair
(
pts[pts.size() - 2].index(), // external 0 -> slave
pts[pts.size() - 1].index()
);
pts.append
(
Vb
(
pt3,
vertexCount() + pts.size(),
Vb::vtExternalSurface,
Pstream::myProcNo()
)
);
ptPairs_.addPointPair
(
pts[pts.size() - 3].index(), // external 0 -> slave
pts[pts.size() - 1].index()
);
pts.append
(
Vb
(
pt4,
vertexCount() + pts.size(),
Vb::vtExternalSurface,
Pstream::myProcNo()
)
);
ptPairs_.addPointPair
(
pts[pts.size() - 3].index(), // external 0 -> slave
pts[pts.size() - 1].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);
//
// createBafflePointPair(ppDist, edgePt + s, n, true, pts);
// createBafflePointPair(ppDist, edgePt - s, n, true, pts);
}
@ -884,7 +1170,10 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints(bool distribute)
const List<Vb>& ftPtVertices = ftPtConformer_.featurePointVertices();
// Insert the created points directly as already distributed.
this->DelaunayMesh<Delaunay>::insertPoints(ftPtVertices);
Map<label> oldToNewIndices =
this->DelaunayMesh<Delaunay>::insertPoints(ftPtVertices, true);
ftPtConformer_.reIndexPointPairs(oldToNewIndices);
label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
reduce(nFeatureVertices, sumOp<label>());

109
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
View File

@ -225,6 +225,7 @@ inline void Foam::conformalVoronoiMesh::createPointPair
const scalar ppDist,
const Foam::point& surfPt,
const vector& n,
const bool ptPair,
DynamicList<Vb>& pts
) const
{
@ -260,14 +261,14 @@ inline void Foam::conformalVoronoiMesh::createPointPair
)
);
// if (ptPair)
// {
// ptPairs_.addPointPair
// (
// pts[pts.size() - 2].index(),
// pts[pts.size() - 1].index() // external 0 -> slave
// );
// }
if (ptPair)
{
ptPairs_.addPointPair
(
pts[pts.size() - 2].index(),
pts[pts.size() - 1].index() // external 0 -> slave
);
}
}
// else
// {
@ -305,6 +306,7 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
const scalar ppDist,
const Foam::point& surfPt,
const vector& n,
const bool ptPair,
DynamicList<Vb>& pts
) const
{
@ -315,8 +317,8 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
Vb
(
surfPt - ppDistn,
vertexCount() + 1 + pts.size(),
Vb::vtInternalSurface,
vertexCount() + pts.size(),
Vb::vtInternalSurfaceBaffle,
Pstream::myProcNo()
)
);
@ -326,11 +328,20 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
Vb
(
surfPt + ppDistn,
vertexCount() + 1 + pts.size(),
Vb::vtInternalSurface,
vertexCount() + pts.size(),
Vb::vtExternalSurfaceBaffle,
Pstream::myProcNo()
)
);
if (ptPair)
{
ptPairs_.addPointPair
(
pts[pts.size() - 2].index(), // external 0 -> slave
pts[pts.size() - 1].index()
);
}
}
@ -341,8 +352,8 @@ inline bool Foam::conformalVoronoiMesh::internalPointIsInside
{
if
(
!geometryToConformTo_.inside(pt)
|| !geometryToConformTo_.globalBounds().contains(pt)
!geometryToConformTo_.globalBounds().contains(pt)
|| !geometryToConformTo_.inside(pt)
)
{
return false;
@ -352,76 +363,6 @@ inline bool Foam::conformalVoronoiMesh::internalPointIsInside
}
inline bool Foam::conformalVoronoiMesh::isPointPair
(
const Vertex_handle& vA,
const Vertex_handle& vB
) const
{
// Want to do this topologically, but problem if vertices are redistributed
// so that one of the point pair is one processor and the other is on
// another.
const Foam::point& ptA = topoint(vA->point());
const Foam::point& ptB = topoint(vB->point());
if
(
(
vA->type() == Vb::vtInternalSurface
&& vB->type() == Vb::vtExternalSurface
)
||
(
vB->type() == Vb::vtInternalSurface
&& vA->type() == Vb::vtExternalSurface
)
||
(
vA->type() == Vb::vtInternalFeatureEdge
&& vB->type() == Vb::vtExternalFeatureEdge
)
||
(
vB->type() == Vb::vtInternalFeatureEdge
&& vA->type() == Vb::vtExternalFeatureEdge
)
||
(
vA->type() == Vb::vtInternalSurface
&& vB->type() == Vb::vtExternalFeatureEdge
)
||
(
vB->type() == Vb::vtInternalSurface
&& vA->type() == Vb::vtExternalFeatureEdge
)
||
(
vA->type() == Vb::vtExternalSurface
&& vB->type() == Vb::vtInternalFeatureEdge
)
||
(
vB->type() == Vb::vtExternalSurface
&& vA->type() == Vb::vtInternalFeatureEdge
)
)
{
const scalar distSqr = magSqr(ptA - ptB);
const scalar ppDistSqr = sqr(2*pointPairDistance(0.5*(ptA + ptB)));
if (distSqr > 1.001*ppDistSqr)
{
return false;
}
}
return true;
}
inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
(
const Delaunay::Finite_edges_iterator& eit

142
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
View File

@ -36,6 +36,8 @@ License
#include "indexedVertexOps.H"
#include "DelaunayMeshTools.H"
#include "syncTools.H"
#include "faceSet.H"
#include "OBJstream.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -114,7 +116,6 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
wordList patchNames;
PtrList<dictionary> patchDicts;
pointField cellCentres;
PackedBoolList boundaryFacesToRemove;
calcDualMesh
@ -789,7 +790,7 @@ void Foam::conformalVoronoiMesh::writeMesh
const wordList& patchNames,
const PtrList<dictionary>& patchDicts,
const pointField& cellCentres,
const PackedBoolList& boundaryFacesToRemove
PackedBoolList& boundaryFacesToRemove
) const
{
if (foamyHexMeshControls().objOutput())
@ -912,42 +913,84 @@ void Foam::conformalVoronoiMesh::writeMesh
// Add indirectPatchFaces to a face zone
Info<< indent << "Add pointZones" << endl;
{
labelList addr(boundaryFacesToRemove.count());
label count = 0;
label sz = mesh.pointZones().size();
forAll(boundaryFacesToRemove, faceI)
DynamicList<label> bPts(boundaryPts.size());
forAll(dualMeshPointTypeNames_, typeI)
{
if
(
boundaryFacesToRemove[faceI]
&& mesh.faceZones().whichZone(faceI) == -1
)
forAll(boundaryPts, ptI)
{
addr[count++] = faceI;
const label& bPtType = boundaryPts[ptI];
if (bPtType == typeI)
{
bPts.append(ptI);
}
}
}
addr.setSize(count);
// syncTools::syncPointList(mesh, bPts, maxEqOp<label>(), -1);
label sz = mesh.faceZones().size();
boolList flip(addr.size(), false);
mesh.faceZones().setSize(sz + 1);
mesh.faceZones().set
(
sz,
new faceZone
Info<< incrIndent << indent
<< "Adding " << bPts.size()
<< " points of type " << dualMeshPointTypeNames_.words()[typeI]
<< decrIndent << endl;
mesh.pointZones().append
(
"indirectPatchFaces",
addr,
flip,
sz,
mesh.faceZones()
)
);
new pointZone
(
dualMeshPointTypeNames_.words()[typeI],
bPts,
sz + typeI,
mesh.pointZones()
)
);
bPts.clear();
}
}
// Add indirectPatchFaces to a face zone
Info<< indent << "Adding indirect patch faces set" << endl;
syncTools::syncFaceList
(
mesh,
boundaryFacesToRemove,
orOp<bool>()
);
labelList addr(boundaryFacesToRemove.count());
label count = 0;
forAll(boundaryFacesToRemove, faceI)
{
if (boundaryFacesToRemove[faceI])
{
addr[count++] = faceI;
}
}
addr.setSize(count);
faceSet indirectPatchFaces
(
mesh,
"indirectPatchFaces",
addr,
IOobject::AUTO_WRITE
);
indirectPatchFaces.sync(mesh);
Info<< decrIndent;
timeCheck("Before fvMesh filtering");
autoPtr<polyMeshFilter> meshFilter;
@ -966,9 +1009,11 @@ void Foam::conformalVoronoiMesh::writeMesh
{
const autoPtr<fvMesh>& newMesh = meshFilter().filteredMesh();
polyTopoChange meshMod(newMesh);
polyTopoChange meshMod(newMesh());
meshMod.changeMesh(mesh, false);
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
polyMeshFilter::copySets(newMesh(), mesh);
}
}
@ -1001,9 +1046,11 @@ void Foam::conformalVoronoiMesh::writeMesh
{
const autoPtr<fvMesh>& newMesh = meshFilter().filteredMesh();
polyTopoChange meshMod(newMesh);
polyTopoChange meshMod(newMesh());
meshMod.changeMesh(mesh, false);
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
polyMeshFilter::copySets(newMesh(), mesh);
}
}
@ -1068,7 +1115,7 @@ void Foam::conformalVoronoiMesh::writeMesh
// writeCellCentres(mesh);
// findRemainingProtrusionSet(mesh);
findRemainingProtrusionSet(mesh);
}
@ -1375,4 +1422,33 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findRemainingProtrusionSet
}
void Foam::conformalVoronoiMesh::writePointPairs
(
const fileName& fName
) const
{
OBJstream os(fName);
for
(
Delaunay::Finite_edges_iterator eit = finite_edges_begin();
eit != finite_edges_end();
++eit
)
{
Cell_handle c = eit->first;
Vertex_handle vA = c->vertex(eit->second);
Vertex_handle vB = c->vertex(eit->third);
if (ptPairs_.isPointPair(vA, vB))
{
os.write
(
linePointRef(topoint(vA->point()), topoint(vB->point()))
);
}
}
}
// ************************************************************************* //

34
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.C
View File

@ -146,7 +146,7 @@ void Foam::featurePointConformer::addMasterAndSlavePoints
)
);
// const label masterIndex = pts[pts.size() - 1].index();
const label masterIndex = pts.last().index();
//meshTools::writeOBJ(strMasters, masterPt);
@ -180,6 +180,11 @@ void Foam::featurePointConformer::addMasterAndSlavePoints
)
);
ftPtPairs_.addPointPair
(
masterIndex,
pts.last().index()
);
//meshTools::writeOBJ(strSlaves, slavePt);
}
@ -525,7 +530,8 @@ Foam::featurePointConformer::featurePointConformer
foamyHexMesh_(foamyHexMesh),
foamyHexMeshControls_(foamyHexMesh.foamyHexMeshControls()),
geometryToConformTo_(foamyHexMesh.geometryToConformTo()),
featurePointVertices_()
featurePointVertices_(),
ftPtPairs_(foamyHexMesh)
{
Info<< nl
<< "Conforming to feature points" << endl;
@ -598,6 +604,30 @@ void Foam::featurePointConformer::distribute
{
featurePointVertices_[vI].procIndex() = Pstream::myProcNo();
}
// Also update the feature point pairs
}
void Foam::featurePointConformer::reIndexPointPairs
(
const Map<label>& oldToNewIndices
)
{
forAll(featurePointVertices_, vI)
{
const label currentIndex = featurePointVertices_[vI].index();
Map<label>::const_iterator newIndexIter =
oldToNewIndices.find(currentIndex);
if (newIndexIter != oldToNewIndices.end())
{
featurePointVertices_[vI].index() = newIndexIter();
}
}
ftPtPairs_.reIndex(oldToNewIndices);
}

10
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.H
View File

@ -44,6 +44,7 @@ SourceFiles
#include "DynamicList.H"
#include "List.H"
#include "extendedFeatureEdgeMesh.H"
#include "pointPairs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -83,6 +84,9 @@ class featurePointConformer
// triangulation clear.
List<Vb> featurePointVertices_;
//-
mutable pointPairs<Delaunay> ftPtPairs_;
// Private Member Functions
@ -163,12 +167,18 @@ public:
// triangulation.
inline const List<Vb>& featurePointVertices() const;
//- Return the feature point pair table
inline const pointPairs<Delaunay>& featurePointPairs() const;
// Edit
//- Distribute the feature point vertices according to the
// supplied background mesh
void distribute(const backgroundMeshDecomposition& decomposition);
//- Reindex the feature point pairs using the map.
void reIndexPointPairs(const Map<label>& oldToNewIndices);
};

6
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerI.H
View File

@ -28,5 +28,11 @@ const Foam::List<Vb>& Foam::featurePointConformer::featurePointVertices() const
return featurePointVertices_;
}
const Foam::pointPairs<Delaunay>&
Foam::featurePointConformer::featurePointPairs() const
{
return ftPtPairs_;
}
// ************************************************************************* //

196
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerSpecialisations.C
View File

@ -168,9 +168,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(internalPtA, Vb::vtInternalFeaturePoint)
Vb
(
internalPtA,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtInternalFeaturePoint,
Pstream::myProcNo()
)
);
const label internalPtAIndex(pts.last().index());
const Foam::point internalPtB =
concaveEdgeExternalPt
- 2.0*planeB.distance(concaveEdgeExternalPt)
@ -178,9 +186,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(internalPtB, Vb::vtInternalFeaturePoint)
Vb
(
internalPtB,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtInternalFeaturePoint,
Pstream::myProcNo()
)
);
const label internalPtBIndex(pts.last().index());
// Add the external points
Foam::point externalPtD;
@ -288,7 +304,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(externalPtD, Vb::vtExternalFeaturePoint)
Vb
(
externalPtD,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtExternalFeaturePoint,
Pstream::myProcNo()
)
);
ftPtPairs_.addPointPair
(
internalPtAIndex,
pts.last().index()
);
}
}
@ -319,7 +347,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(externalPtE, Vb::vtExternalFeaturePoint)
Vb
(
externalPtE,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtExternalFeaturePoint,
Pstream::myProcNo()
)
);
ftPtPairs_.addPointPair
(
internalPtBIndex,
pts.last().index()
);
}
}
@ -328,9 +368,29 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(concaveEdgeExternalPt, Vb::vtExternalFeaturePoint)
Vb
(
concaveEdgeExternalPt,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtExternalFeaturePoint,
Pstream::myProcNo()
)
);
ftPtPairs_.addPointPair
(
internalPtBIndex,
pts.last().index()
);
ftPtPairs_.addPointPair
(
internalPtAIndex,
pts.last().index()
);
const label concaveEdgeExternalPtIndex(pts.last().index());
const scalar totalAngle = radToDeg
(
constant::mathematical::pi
@ -377,7 +437,21 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(internalPtF, Vb::vtInternalFeaturePoint)
Vb
(
internalPtF,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtInternalFeaturePoint,
Pstream::myProcNo()
)
);
const label internalPtFIndex(pts.last().index());
ftPtPairs_.addPointPair
(
concaveEdgeExternalPtIndex,
pts.last().index()
);
const Foam::point externalPtG =
@ -386,7 +460,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(externalPtG, Vb::vtExternalFeaturePoint)
Vb
(
externalPtG,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtExternalFeaturePoint,
Pstream::myProcNo()
)
);
ftPtPairs_.addPointPair
(
internalPtFIndex,
pts.last().index()
);
}
@ -520,9 +606,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(internalPtA, Vb::vtExternalFeaturePoint)
Vb
(
internalPtA,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtExternalFeaturePoint,
Pstream::myProcNo()
)
);
const label internalPtAIndex(pts.last().index());
const Foam::point internalPtB =
convexEdgeExternalPt
+ 2.0*planeB.distance(convexEdgeExternalPt)
@ -530,9 +624,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(internalPtB, Vb::vtExternalFeaturePoint)
Vb
(
internalPtB,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtExternalFeaturePoint,
Pstream::myProcNo()
)
);
const label internalPtBIndex(pts.last().index());
// Add the internal points
Foam::point externalPtD;
@ -643,7 +745,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(externalPtD, Vb::vtInternalFeaturePoint)
Vb
(
externalPtD,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtInternalFeaturePoint,
Pstream::myProcNo()
)
);
ftPtPairs_.addPointPair
(
internalPtAIndex,
pts.last().index()
);
}
}
@ -674,7 +788,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(externalPtE, Vb::vtInternalFeaturePoint)
Vb
(
externalPtE,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtInternalFeaturePoint,
Pstream::myProcNo()
)
);
ftPtPairs_.addPointPair
(
internalPtBIndex,
pts.last().index()
);
}
}
@ -683,7 +809,25 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(convexEdgeExternalPt, Vb::vtInternalFeaturePoint)
Vb
(
convexEdgeExternalPt,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtInternalFeaturePoint,
Pstream::myProcNo()
)
);
ftPtPairs_.addPointPair
(
internalPtBIndex,
pts.last().index()
);
ftPtPairs_.addPointPair
(
internalPtAIndex,
pts.last().index()
);
const scalar totalAngle = radToDeg
@ -732,7 +876,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(internalPtF, Vb::vtExternalFeaturePoint)
Vb
(
internalPtF,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtExternalFeaturePoint,
Pstream::myProcNo()
)
);
ftPtPairs_.addPointPair
(
pts[pts.size() - 2].index(),
pts.last().index()
);
const Foam::point externalPtG =
@ -741,7 +897,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
pts.append
(
Vb(externalPtG, Vb::vtInternalFeaturePoint)
Vb
(
externalPtG,
foamyHexMesh_.vertexCount() + pts.size(),
Vb::vtInternalFeaturePoint,
Pstream::myProcNo()
)
);
ftPtPairs_.addPointPair
(
pts[pts.size() - 2].index(),
pts.last().index()
);
}

6
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H
View File

@ -193,6 +193,8 @@ public:
//- Does the Dual vertex form part of a processor patch
inline bool parallelDualVertex() const;
inline Foam::label vertexLowestProc() 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
@ -219,7 +221,9 @@ public:
// least one Delaunay vertex outside and at least one inside
inline bool boundaryDualVertex() const;
inline bool baffleBoundaryDualVertex() const;
inline bool baffleSurfaceDualVertex() const;
inline bool baffleEdgeDualVertex() const;
//- A dual vertex on a feature edge will result from this Delaunay cell
inline bool featureEdgeDualVertex() const;

56
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H
View File

@ -291,6 +291,23 @@ inline bool CGAL::indexedCell<Gt, Cb>::parallelDualVertex() const
}
template<class Gt, class Cb>
inline Foam::label CGAL::indexedCell<Gt, Cb>::vertexLowestProc() const
{
Foam::label lowestProc = -1;
for (int i = 0; i < 4; ++i)
{
if (this->vertex(i)->referred())
{
lowestProc = min(lowestProc, this->vertex(i)->procIndex());
}
}
return lowestProc;
}
template<class Gt, class Cb>
inline Foam::tetCell CGAL::indexedCell<Gt, Cb>::vertexGlobalIndices
(
@ -426,21 +443,42 @@ inline bool CGAL::indexedCell<Gt, Cb>::boundaryDualVertex() const
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::baffleBoundaryDualVertex() const
inline bool CGAL::indexedCell<Gt, Cb>::baffleSurfaceDualVertex() const
{
return
(
(
this->vertex(0)->internalBafflePoint()
|| this->vertex(1)->internalBafflePoint()
|| this->vertex(2)->internalBafflePoint()
|| this->vertex(3)->internalBafflePoint()
this->vertex(0)->internalBaffleSurfacePoint()
|| this->vertex(1)->internalBaffleSurfacePoint()
|| this->vertex(2)->internalBaffleSurfacePoint()
|| this->vertex(3)->internalBaffleSurfacePoint()
)
&& (
this->vertex(0)->externalBafflePoint()
|| this->vertex(1)->externalBafflePoint()
|| this->vertex(2)->externalBafflePoint()
|| this->vertex(3)->externalBafflePoint()
this->vertex(0)->externalBaffleSurfacePoint()
|| this->vertex(1)->externalBaffleSurfacePoint()
|| this->vertex(2)->externalBaffleSurfacePoint()
|| this->vertex(3)->externalBaffleSurfacePoint()
)
);
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::baffleEdgeDualVertex() const
{
return
(
(
this->vertex(0)->internalBaffleEdgePoint()
|| this->vertex(1)->internalBaffleEdgePoint()
|| this->vertex(2)->internalBaffleEdgePoint()
|| this->vertex(3)->internalBaffleEdgePoint()
)
&& (
this->vertex(0)->externalBaffleEdgePoint()
|| this->vertex(1)->externalBaffleEdgePoint()
|| this->vertex(2)->externalBaffleEdgePoint()
|| this->vertex(3)->externalBaffleEdgePoint()
)
);
}

6
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertex.H
View File

@ -243,10 +243,12 @@ public:
inline bool surfacePoint() const;
inline bool internalBoundaryPoint() const;
inline bool internalBafflePoint() const;
inline bool internalBaffleSurfacePoint() const;
inline bool internalBaffleEdgePoint() const;
inline bool externalBoundaryPoint() const;
inline bool externalBafflePoint() const;
inline bool externalBaffleSurfacePoint() const;
inline bool externalBaffleEdgePoint() const;
inline bool constrained() const;

27
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexI.H
View File

@ -308,15 +308,16 @@ inline bool CGAL::indexedVertex<Gt, Vb>::internalBoundaryPoint() const
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::internalBafflePoint() const
inline bool CGAL::indexedVertex<Gt, Vb>::internalBaffleSurfacePoint() const
{
return
(
type_ == vtInternalSurfaceBaffle
|| type_ == vtInternalFeatureEdgeBaffle
);
return (type_ == vtInternalSurfaceBaffle);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::internalBaffleEdgePoint() const
{
return (type_ == vtInternalFeatureEdgeBaffle);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::externalBoundaryPoint() const
@ -325,13 +326,15 @@ inline bool CGAL::indexedVertex<Gt, Vb>::externalBoundaryPoint() const
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::externalBafflePoint() const
inline bool CGAL::indexedVertex<Gt, Vb>::externalBaffleSurfacePoint() const
{
return
(
type_ == vtExternalSurfaceBaffle
|| type_ == vtExternalFeatureEdgeBaffle
);
return (type_ == vtExternalSurfaceBaffle);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::externalBaffleEdgePoint() const
{
return (type_ == vtExternalFeatureEdgeBaffle);
}

29
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
View File

@ -83,6 +83,13 @@ void Foam::conformationSurfaces::hasBoundedVolume
const pointField& surfPts = triSurf.points();
Info<< " Checking " << surface.name() << endl;
label nBaffles = 0;
Info<< " Index = " << surfaces_[s] << endl;
Info<< " Offset = " << regionOffset_[s] << endl;
forAll(triSurf, sI)
{
const label patchID =
@ -98,7 +105,13 @@ void Foam::conformationSurfaces::hasBoundedVolume
{
sum += triSurf[sI].normal(surfPts);
}
else
{
nBaffles++;
}
}
Info<< " has " << nBaffles << " baffles out of "
<< triSurf.size() << " triangles" << endl;
totalTriangles += triSurf.size();
}
@ -153,8 +166,9 @@ void Foam::conformationSurfaces::readFeatures
{
const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
Info<< " features: " << surface.name() << " of type "
<< surface.type() << endl;
Info<< " features: " << surface.name()
<< " of type " << surface.type()
<< ", id: " << featureIndex << endl;
autoPtr<searchableSurfaceFeatures> ssFeatures
(
@ -210,7 +224,8 @@ void Foam::conformationSurfaces::readFeatures
{
fileName feMeshName(featureDict.lookup("extendedFeatureEdgeMesh"));
Info<< " features: " << feMeshName << endl;
Info<< " features: " << feMeshName << ", id: " << featureIndex
<< endl;
features_.set
(
@ -716,7 +731,11 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
const searchableSurface& surface(allGeometry_[surfaces_[s]]);
if (!surface.hasVolumeType())
if
(
!surface.hasVolumeType()
//&& surfaceVolumeTests[s][i] == volumeType::UNKNOWN
)
{
pointField sample(1, samplePts[i]);
scalarField nearestDistSqr(1, GREAT);
@ -733,7 +752,7 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
findSurfaceNearestIntersection
(
samplePts[i],
info[0].rawPoint(),
info[0].rawPoint() - 1e-3*mag(hitDir)*hitDir,
surfHit,
hitSurface
);

20
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
View File

@ -187,7 +187,7 @@ Foam::label Foam::autoDensity::recurseAndFill
word recursionName
) const
{
label maxDepth = 0;
label treeDepth = 0;
for (direction i = 0; i < 8; i++)
{
@ -201,10 +201,10 @@ Foam::label Foam::autoDensity::recurseAndFill
{
if (levelLimit > 0)
{
maxDepth =
treeDepth =
max
(
maxDepth,
treeDepth,
recurseAndFill
(
initialPoints,
@ -229,10 +229,10 @@ Foam::label Foam::autoDensity::recurseAndFill
if (!fillBox(initialPoints, subBB, true))
{
maxDepth =
treeDepth =
max
(
maxDepth,
treeDepth,
recurseAndFill
(
initialPoints,
@ -259,10 +259,10 @@ Foam::label Foam::autoDensity::recurseAndFill
if (!fillBox(initialPoints, subBB, false))
{
maxDepth =
treeDepth =
max
(
maxDepth,
treeDepth,
recurseAndFill
(
initialPoints,
@ -286,7 +286,7 @@ Foam::label Foam::autoDensity::recurseAndFill
}
}
return maxDepth + 1;
return treeDepth + 1;
}
@ -941,7 +941,7 @@ List<Vb::Point> autoDensity::initialPoints() const
Pout<< " Filling box " << hierBB << endl;
}
label maxDepth = recurseAndFill
label treeDepth = recurseAndFill
(
initialPoints,
hierBB,
@ -966,7 +966,7 @@ List<Vb::Point> autoDensity::initialPoints() const
<< scalar(nInitialPoints)/scalar(max(globalTrialPoints_, 1))
<< " success rate" << nl
<< indent
<< returnReduce(maxDepth, maxOp<label>())
<< returnReduce(treeDepth, maxOp<label>())
<< " levels of recursion (maximum)"
<< decrIndent << decrIndent
<< endl;

244
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.C Normal file
View File

@ -0,0 +1,244 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "pointPairs.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class Triangulation>
inline Foam::Pair<Foam::labelPair>
Foam::pointPairs<Triangulation>::orderPointPair
(
const labelPair& vA,
const labelPair& vB
) const
{
return
(
(vA < vB)
? Pair<labelPair>(vA, vB)
: Pair<labelPair>(vB, vA)
);
}
template<class Triangulation>
inline bool Foam::pointPairs<Triangulation>::findPointPair
(
const labelPair& vA,
const labelPair& vB
) const
{
if (vA == vB)
{
return false;
}
else if (find(orderPointPair(vA, vB)) == end())
{
return false;
}
return true;
}
template<class Triangulation>
inline bool Foam::pointPairs<Triangulation>::insertPointPair
(
const labelPair& vA,
const labelPair& vB
)
{
if (vA == vB)
{
return false;
}
return insert(orderPointPair(vA, vB));
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Triangulation>
Foam::pointPairs<Triangulation>::pointPairs(const Triangulation& triangulation)
:
ptPairTable(),
triangulation_(triangulation)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Triangulation>
Foam::pointPairs<Triangulation>::~pointPairs()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class Triangulation>
inline bool Foam::pointPairs<Triangulation>::addPointPair
(
const labelPair& vA,
const labelPair& vB
)
{
return insertPointPair(vA, vB);
}
template<class Triangulation>
inline bool Foam::pointPairs<Triangulation>::addPointPair
(
const labelPair& master,
const DynamicList<labelPair>& slaves
)
{
forAll(slaves, sI)
{
const labelPair& slave = slaves[sI];
addPointPair(master, slave);
}
return true;
}
template<class Triangulation>
inline bool Foam::pointPairs<Triangulation>::addPointPair
(
const label vA,
const label vB
)
{
const label procNo = Pstream::myProcNo();
labelPair a(vA, procNo);
labelPair b(vB, procNo);
return addPointPair(a, b);
}
template<class Triangulation>
inline bool Foam::pointPairs<Triangulation>::isPointPair
(
const Vertex_handle& vA,
const Vertex_handle& vB
) const
{
if (vA->boundaryPoint() && vB->boundaryPoint())
{
labelPair a(vA->index(), vA->procIndex());
labelPair b(vB->index(), vB->procIndex());
return findPointPair(a, b);
}
return false;
}
template<class Triangulation>
inline bool Foam::pointPairs<Triangulation>::isPointPair
(
const labelPair& vA,
const labelPair& vB
) const
{
return findPointPair(vA, vB);
}
template<class Triangulation>
void Foam::pointPairs<Triangulation>::reIndex(const Map<label>& oldToNewIndices)
{
pointPairs<Triangulation> newTable(triangulation_);
forAllConstIter(pointPairs, *this, iter)
{
Pair<labelPair> e = iter.key();
labelPair& start = e.first();
labelPair& end = e.second();
bool insert = true;
if (start.second() == Pstream::myProcNo())
{
Map<label>::const_iterator iter2 =
oldToNewIndices.find(start.first());
if (iter2 != oldToNewIndices.end())
{
if (iter2() != -1)
{
start.first() = iter2();
}
else
{
insert = false;
}
}
}
if (end.second() == Pstream::myProcNo())
{
Map<label>::const_iterator iter2 =
oldToNewIndices.find(end.first());
if (iter2 != oldToNewIndices.end())
{
if (iter2() != -1)
{
end.first() = iter2();
}
else
{
insert = false;
}
}
}
if (insert)
{
if (e.first() < e.second())
{
newTable.insert(e);
}
else if (e.first() > e.second())
{
newTable.insert(reverse(e));
}
}
}
this->transfer(newTable);
}
// ************************************************************************* //

163
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.H Normal file
View File

@ -0,0 +1,163 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
pointPairs
Description
HashSet of unique edges. The edges are stored as a pair of pairs:
( (local index, processor index) (local index, processor index) )
e.g.,
( (0 1) (3 1) )
( (0 2) (5 1) )
\*---------------------------------------------------------------------------*/
#ifndef pointPairs_H
#define pointPairs_H
#include "labelPair.H"
#include "HashSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
typedef HashSet
<
Pair<labelPair>,
FixedList<labelPair, 2>::Hash<>
> ptPairTable;
/*---------------------------------------------------------------------------*\
Class pointPairs Declaration
\*---------------------------------------------------------------------------*/
template<class Triangulation>
class pointPairs
:
public ptPairTable
{
// Private typedefs
typedef typename Triangulation::Vertex_handle Vertex_handle;
// Private data
const Triangulation& triangulation_;
// Private Member Functions
inline Pair<labelPair> orderPointPair
(
const labelPair& vA,
const labelPair& vB
) const;
inline bool insertPointPair
(
const labelPair& vA,
const labelPair& vB
);
inline bool findPointPair
(
const labelPair& vA,
const labelPair& vB
) const;
public:
// Constructors
//- Construct from triangulation
pointPairs(const Triangulation& triangulation);
//- Destructor
~pointPairs();
// Member Functions
// Access
inline bool isPointPair
(
const Vertex_handle& vA,
const Vertex_handle& vB
) const;
inline bool isPointPair
(
const labelPair& vA,
const labelPair& vB
) const;
// Edit
inline bool addPointPair
(
const labelPair& vA,
const labelPair& vB
);
inline bool addPointPair
(
const labelPair& master,
const DynamicList<labelPair>& slaves
);
inline bool addPointPair
(
const label vA,
const label vB
);
void reIndex(const Map<label>& oldToNewIndices);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "pointPairs.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

11
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchableBoxFeatures.C
View File

@ -104,11 +104,14 @@ Foam::searchableBoxFeatures::features() const
edgeNormals[10][0] = 1; edgeNormals[10][1] = 3;
edgeNormals[11][0] = 3; edgeNormals[11][1] = 0;
tmp<pointField> surfacePointsTmp(surface().points());
pointField& surfacePoints = surfacePointsTmp();
forAll(edgeDirections, eI)
{
edgeDirections[eI] =
surface().points()()[treeBoundBox::edges[eI].end()]
- surface().points()()[treeBoundBox::edges[eI].start()];
surfacePoints[treeBoundBox::edges[eI].end()]
- surfacePoints[treeBoundBox::edges[eI].start()];
normalDirections[eI] = labelList(2, 0);
for (label j = 0; j < 2; ++j)
@ -116,8 +119,8 @@ Foam::searchableBoxFeatures::features() const
const vector cross =
(faceNormals[edgeNormals[eI][j]] ^ edgeDirections[eI]);
const vector fC0tofE0 =
0.5*(max(surface().points()() + min(surface().points()())))
- surface().points()()[treeBoundBox::edges[eI].start()];
0.5*(max(surfacePoints + min(surfacePoints)))
- surfacePoints[treeBoundBox::edges[eI].start()];
normalDirections[eI][j] =
(

218
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchablePlateFeatures.C Normal file
View File

@ -0,0 +1,218 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "searchablePlateFeatures.H"
#include "addToRunTimeSelectionTable.H"
#include "treeBoundBox.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchablePlateFeatures, 0);
addToRunTimeSelectionTable
(
searchableSurfaceFeatures,
searchablePlateFeatures,
dict
);
//! \cond - skip documentation : local scope only
const Foam::label edgesArray[4][2] =
{
{0, 1}, // 0
{0, 3},
{2, 1}, // 2
{2, 3}
};
//! \endcond
const edgeList searchablePlateFeatures::edges(calcEdges(edgesArray));
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::edgeList Foam::searchablePlateFeatures::calcEdges
(
const label edgesArray[4][2]
)
{
edgeList edges(4);
forAll(edges, edgeI)
{
edges[edgeI][0] = edgesArray[edgeI][0];
edges[edgeI][1] = edgesArray[edgeI][1];
}
return edges;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchablePlateFeatures::searchablePlateFeatures
(
const searchableSurface& surface,
const dictionary& dict
)
:
searchableSurfaceFeatures(surface, dict),
mode_
(
extendedFeatureEdgeMesh::sideVolumeTypeNames_
[
dict.lookupOrDefault<word>("meshableSide", "inside")
]
)
{
Info<< indent
<< " Meshable region = "
<< extendedFeatureEdgeMesh::sideVolumeTypeNames_[mode_]
<< endl;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::searchablePlateFeatures::~searchablePlateFeatures()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::autoPtr<Foam::extendedFeatureEdgeMesh>
Foam::searchablePlateFeatures::features() const
{
autoPtr<extendedFeatureEdgeMesh> features;
vectorField faceNormals(1);
surface().getNormal(List<pointIndexHit>(1), faceNormals);
vectorField edgeDirections(4);
labelListList normalDirections(4);
labelListList edgeNormals(4);
forAll(edgeNormals, eI)
{
edgeNormals[eI].setSize(2, 0);
}
edgeNormals[0][0] = 0; edgeNormals[0][1] = 0;
edgeNormals[1][0] = 0; edgeNormals[1][1] = 0;
edgeNormals[2][0] = 0; edgeNormals[2][1] = 0;
edgeNormals[3][0] = 0; edgeNormals[3][1] = 0;
forAll(edgeDirections, eI)
{
edgeDirections[eI] =
surface().points()()[edges[eI].end()]
- surface().points()()[edges[eI].start()];
normalDirections[eI] = labelList(2, 0);
for (label j = 0; j < 2; ++j)
{
const vector cross =
(faceNormals[edgeNormals[eI][j]] ^ edgeDirections[eI]);
const vector fC0tofE0 =
0.5*(max(surface().points()() + min(surface().points()())))
- surface().points()()[edges[eI].start()];
normalDirections[eI][j] =
(
(
(cross/(mag(cross) + VSMALL))
& (fC0tofE0/(mag(fC0tofE0)+ VSMALL))
)
> 0.0
? 1
: -1
);
}
}
labelListList featurePointNormals(4);
labelListList featurePointEdges(4);
forAll(featurePointNormals, pI)
{
labelList& ftPtEdges = featurePointEdges[pI];
ftPtEdges.setSize(2, 0);
label edgeI = 0;
forAll(edges, eI)
{
const edge& e = edges[eI];
if (e.start() == pI)
{
ftPtEdges[edgeI++] = eI;
}
else if (e.end() == pI)
{
ftPtEdges[edgeI++] = eI;
}
}
labelList& ftPtNormals = featurePointNormals[pI];
ftPtNormals.setSize(1, 0);
ftPtNormals[0] = edgeNormals[ftPtEdges[0]][0];
}
labelList regionEdges;
features.set
(
new extendedFeatureEdgeMesh
(
IOobject
(
surface().name(),
surface().instance(),
"extendedFeatureEdgeMesh",
surface().db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
surface().points(),
edges,
4, 4, 4,
0, 0, 4, 4, // 4 flat edges
faceNormals,
List<extendedFeatureEdgeMesh::sideVolumeType>(4, mode_),
edgeDirections,
normalDirections,
edgeNormals,
featurePointNormals,
featurePointEdges,
regionEdges
)
);
return features;
}
// ************************************************************************* //

120
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchablePlateFeatures.H Normal file
View File

@ -0,0 +1,120 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::searchablePlateFeatures
Description
SourceFiles
searchablePlateFeatures.C
\*---------------------------------------------------------------------------*/
#ifndef searchablePlateFeatures_H
#define searchablePlateFeatures_H
#include "searchableSurfaceFeatures.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class searchablePlateFeatures Declaration
\*---------------------------------------------------------------------------*/
class searchablePlateFeatures
:
public searchableSurfaceFeatures
{
private:
//- To initialise edges.
static edgeList calcEdges(const label[4][2]);
// Private Member Data
//- Which side of the plate to mesh
const extendedFeatureEdgeMesh::sideVolumeType mode_;
// Private Member Functions
//- Disallow default bitwise copy construct
searchablePlateFeatures(const searchablePlateFeatures&);
//- Disallow default bitwise assignment
void operator=(const searchablePlateFeatures&);
public:
//- Runtime type information
TypeName("searchablePlateFeatures");
// Static data members
//- Edge to point addressing
static const edgeList edges;
// Constructors
//- Construct from searchable surface and dictionary
searchablePlateFeatures
(
const searchableSurface& surface,
const dictionary& dict
);
//- Destructor
virtual ~searchablePlateFeatures();
// Member Functions
//- Return true for a searchable plate having features
virtual bool hasFeatures() const
{
return true;
}
//- Return an extendedFeatureEdgeMesh containing the features
virtual autoPtr<extendedFeatureEdgeMesh> features() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

24
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/triSurfaceMeshFeatures.C
View File

@ -57,10 +57,19 @@ Foam::triSurfaceMeshFeatures::triSurfaceMeshFeatures
)
:
searchableSurfaceFeatures(surface, dict),
includedAngle_(readScalar(dict.lookup("includedAngle")))
includedAngle_(readScalar(dict.lookup("includedAngle"))),
mode_
(
extendedFeatureEdgeMesh::sideVolumeTypeNames_
[
dict.lookupOrDefault<word>("meshableSide", "inside")
]
)
{
Info<< indent
<< " Included angle = " << includedAngle_
<< " Included angle = " << includedAngle_ << nl
<< " Meshable region = "
<< extendedFeatureEdgeMesh::sideVolumeTypeNames_[mode_]
<< endl;
}
@ -82,7 +91,12 @@ Foam::triSurfaceMeshFeatures::features() const
surfaceFeatures sFeat(surfMesh, includedAngle_);
boolList surfBaffleRegions(surfMesh.patches().size(), false);
// @todo Need to read on a per region basis
boolList surfBaffleRegions
(
surfMesh.patches().size(),
(mode_ == extendedFeatureEdgeMesh::BOTH ? true : false)
);
features.set
(
@ -90,8 +104,8 @@ Foam::triSurfaceMeshFeatures::features() const
(
sFeat,
surface().db(),
surface().name() + ".extendedFeatureEdgeMesh"//,
//surfBaffleRegions
surface().name() + ".extendedFeatureEdgeMesh",
surfBaffleRegions
)
);

3
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/triSurfaceMeshFeatures.H
View File

@ -55,6 +55,9 @@ private:
const scalar includedAngle_;
//- Which side of the surface to mesh
const extendedFeatureEdgeMesh::sideVolumeType mode_;
// Private Member Functions

12
applications/utilities/mesh/generation/foamyHexMesh/foamyHexMesh.C
View File

@ -54,8 +54,6 @@ int main(int argc, char *argv[])
"conform to the initial points without any point motion"
);
#include "addOverwriteOption.H"
#include "setRootCase.H"
#include "createTime.H"
@ -63,7 +61,6 @@ int main(int argc, char *argv[])
const bool checkGeometry = args.optionFound("checkGeometry");
const bool conformationOnly = args.optionFound("conformationOnly");
const bool overwrite = args.optionFound("overwrite");
IOdictionary foamyHexMeshDict
(
@ -123,10 +120,7 @@ int main(int argc, char *argv[])
{
mesh.initialiseForConformation();
if (!overwrite)
{
runTime++;
}
runTime++;
mesh.writeMesh(runTime.timeName());
}
@ -134,8 +128,10 @@ int main(int argc, char *argv[])
{
mesh.initialiseForMotion();
while (runTime.loop())
while (runTime.run())
{
runTime++;
Info<< nl << "Time = " << runTime.timeName() << endl;
mesh.move();

1
applications/utilities/mesh/generation/foamyQuadMesh/Make/options
View File

@ -9,6 +9,7 @@ EXE_INC = \
${EXE_FROUNDING_MATH} \
${EXE_NDEBUG} \
${CGAL_INC} \
${c++CGALWARN} \
-I$(FOAM_APP)/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/lnInclude \
-I../cvMesh/vectorTools \
-IconformalVoronoi2DMesh/lnInclude \

1
applications/utilities/mesh/manipulation/orientFaceZone/Make/files
View File

@ -1,4 +1,3 @@
patchFaceOrientation.C
orientFaceZone.C
EXE = $(FOAM_APPBIN)/orientFaceZone

2
applications/utilities/mesh/manipulation/orientFaceZone/Make/options
View File

@ -1,7 +1,9 @@
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/mesh/autoMesh/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude
EXE_LIBS = \
-lmeshTools \
-lautoMesh \
-ltriSurface

285
applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
View File

@ -78,11 +78,79 @@ Description
#include "booleanSurface.H"
#include "edgeIntersections.H"
#include "meshTools.H"
#include "labelPair.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
bool intersectSurfaces
(
triSurface& surf,
edgeIntersections& edgeCuts
)
{
bool hasMoved = false;
for (label iter = 0; iter < 10; iter++)
{
Info<< "Determining intersections of surface edges with itself" << endl;
// Determine surface edge intersections. Allow surface to be moved.
// Number of iterations needed to resolve degenerates
label nIters = 0;
{
triSurfaceSearch querySurf(surf);
scalarField surfPointTol
(
1e-3*edgeIntersections::minEdgeLength(surf)
);
// Determine raw intersections
edgeCuts = edgeIntersections
(
surf,
querySurf,
surfPointTol
);
// Shuffle a bit to resolve degenerate edge-face hits
{
pointField points(surf.points());
nIters =
edgeCuts.removeDegenerates
(
5, // max iterations
surf,
querySurf,
surfPointTol,
points // work array
);
if (nIters != 0)
{
// Update geometric quantities
surf.movePoints(points);
hasMoved = true;
}
}
}
}
if (hasMoved)
{
fileName newFile("surf.obj");
Info<< "Surface has been moved. Writing to " << newFile << endl;
surf.write(newFile);
}
return hasMoved;
}
// Keep on shuffling surface points until no more degenerate intersections.
// Moves both surfaces and updates set of edge cuts.
bool intersectSurfaces
@ -238,6 +306,106 @@ label calcNormalDirection
}
void calcEdgeCuts
(
triSurface& surf1,
triSurface& surf2,
const bool perturb,
edgeIntersections& edge1Cuts,
edgeIntersections& edge2Cuts
)
{
if (perturb)
{
intersectSurfaces
(
surf1,
edge1Cuts,
surf2,
edge2Cuts
);
}
else
{
triSurfaceSearch querySurf2(surf2);
Info<< "Determining intersections of surf1 edges with surf2 faces"
<< endl;
edge1Cuts = edgeIntersections
(
surf1,
querySurf2,
1e-3*edgeIntersections::minEdgeLength(surf1)
);
triSurfaceSearch querySurf1(surf1);
Info<< "Determining intersections of surf2 edges with surf1 faces"
<< endl;
edge2Cuts = edgeIntersections
(
surf2,
querySurf1,
1e-3*edgeIntersections::minEdgeLength(surf2)
);
}
}
void calcFeaturePoints(const pointField& points, const edgeList& edges)
{
edgeMesh eMesh(points, edges);
const labelListList& pointEdges = eMesh.pointEdges();
// Get total number of feature points
label nFeaturePoints = 0;
forAll(pointEdges, pI)
{
const labelList& pEdges = pointEdges[pI];
if (pEdges.size() == 1)
{
nFeaturePoints++;
}
}
// Calculate addressing from feature point to cut point and cut edge
labelList featurePointToCutPoint(nFeaturePoints);
labelList featurePointToCutEdge(nFeaturePoints);
label nFeatPts = 0;
forAll(pointEdges, pI)
{
const labelList& pEdges = pointEdges[pI];
if (pEdges.size() == 1)
{
featurePointToCutPoint[nFeatPts] = pI;
featurePointToCutEdge[nFeatPts] = pEdges[0];
nFeatPts++;
}
}
label concaveStart = 0;
label mixedStart = 0;
label nonFeatureStart = nFeaturePoints;
labelListList featurePointNormals(nFeaturePoints);
labelListList featurePointEdges(nFeaturePoints);
labelList regionEdges;
}
int main(int argc, char *argv[])
{
argList::noParallel();
@ -319,52 +487,33 @@ int main(int argc, char *argv[])
<< exit(FatalError);
}
if (args.optionFound("perturb"))
{
intersectSurfaces
(
surf1,
edge1Cuts,
surf2,
edge2Cuts
);
}
else
{
triSurfaceSearch querySurf2(surf2);
// Calculate the points where the edges are cut by the other surface
calcEdgeCuts
(
surf1,
surf2,
args.optionFound("perturb"),
edge1Cuts,
edge2Cuts
);
Info<< "Determining intersections of surf1 edges with surf2 faces"
<< endl;
// Determine intersection edges from the edge cuts
surfaceIntersection inter
(
surf1,
edge1Cuts,
surf2,
edge2Cuts
);
edge1Cuts = edgeIntersections
(
surf1,
querySurf2,
1e-3*edgeIntersections::minEdgeLength(surf1)
);
triSurfaceSearch querySurf1(surf1);
Info<< "Determining intersections of surf2 edges with surf1 faces"
<< endl;
edge2Cuts = edgeIntersections
(
surf2,
querySurf1,
1e-3*edgeIntersections::minEdgeLength(surf2)
);
}
// Determine intersection edges
surfaceIntersection inter(surf1, edge1Cuts, surf2, edge2Cuts);
fileName sFeatFileName =
fileName sFeatFileName
(
surf1Name.lessExt().name()
+ "_"
+ surf2Name.lessExt().name()
+ "_"
+ action;
+ action
);
label nFeatEds = inter.cutEdges().size();
@ -393,6 +542,7 @@ int main(int argc, char *argv[])
edgeDirections[cutEdgeI] = fE.vec(inter.cutPoints());
normals.append(norm1);
eNormals.append(normals.size() - 1);
if (surf1Baffle)
{
@ -416,9 +566,8 @@ int main(int argc, char *argv[])
);
}
eNormals.append(normals.size() - 1);
normals.append(norm2);
eNormals.append(normals.size() - 1);
if (surf2Baffle)
{
@ -443,8 +592,6 @@ int main(int argc, char *argv[])
);
}
eNormals.append(normals.size() - 1);
if (surf1Baffle)
{
@ -509,6 +656,38 @@ int main(int argc, char *argv[])
label internalStart = -1;
label nIntOrExt = 0;
label nFlat = 0;
label nOpen = 0;
label nMultiple = 0;
forAll(edgeNormals, eI)
{
label nEdNorms = edgeNormals[eI].size();
if (nEdNorms == 1)
{
nOpen++;
}
else if (nEdNorms == 2)
{
const vector& n0(normals[edgeNormals[eI][0]]);
const vector& n1(normals[edgeNormals[eI][1]]);
if ((n0 & n1) > extendedFeatureEdgeMesh::cosNormalAngleTol_)
{
nFlat++;
}
else
{
nIntOrExt++;
}
}
else if (nEdNorms > 2)
{
nMultiple++;
}
}
if (validActions[action] == booleanSurface::UNION)
{
@ -520,7 +699,7 @@ int main(int argc, char *argv[])
else
{
// All edges are external
internalStart = nFeatEds;
internalStart = nIntOrExt;
}
}
else if (validActions[action] == booleanSurface::INTERSECTION)
@ -528,7 +707,7 @@ int main(int argc, char *argv[])
if (!invertedSpace)
{
// All edges are external
internalStart = nFeatEds;
internalStart = nIntOrExt;
}
else
{
@ -539,7 +718,7 @@ int main(int argc, char *argv[])
else if (validActions[action] == booleanSurface::DIFFERENCE)
{
// All edges are external
internalStart = nFeatEds;
internalStart = nIntOrExt;
}
else
{
@ -569,6 +748,8 @@ int main(int argc, char *argv[])
normalDirectionsTmp[i] = normalDirections[i];
}
calcFeaturePoints(inter.cutPoints(), inter.cutEdges());
extendedFeatureEdgeMesh feMesh
(
IOobject
@ -582,18 +763,22 @@ int main(int argc, char *argv[])
),
inter.cutPoints(),
inter.cutEdges(),
0, // concaveStart,
0, // mixedStart,
0, // nonFeatureStart,
internalStart, // internalStart,
nFeatEds, // flatStart,
nFeatEds, // openStart,
nFeatEds, // multipleStart,
nIntOrExt, // flatStart,
nIntOrExt + nFlat, // openStart,
nIntOrExt + nFlat + nOpen, // multipleStart,
normalsTmp,
normalVolumeTypesTmp,
edgeDirections,
normalDirectionsTmp,
edgeNormalsTmp,
labelListList(0), // featurePointNormals,
labelListList(0), // featurePointEdges,
labelList(0) // regionEdges

29
applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
View File

@ -925,7 +925,7 @@ void writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os)
<< " open edges : "
<< (fem.multipleStart()- fem.openStart()) << nl
<< " multiply connected : "
<< (fem.edges().size()- fem.multipleStart()) << nl;
<< (fem.edges().size()- fem.multipleStart()) << endl;
}
@ -1223,17 +1223,34 @@ int main(int argc, char *argv[])
<< " internal edges : " << newSet.nInternalEdges() << nl
<< endl;
//if (writeObj)
//{
// newSet.writeObj("final");
//}
if (writeObj)
{
newSet.writeObj("final");
}
boolList surfBaffleRegions(surf.patches().size(), false);
wordList surfBaffleNames;
surfaceDict.readIfPresent("baffles", surfBaffleNames);
forAll(surf.patches(), pI)
{
const word& name = surf.patches()[pI].name();
if (findIndex(surfBaffleNames, name) != -1)
{
Info<< "Adding baffle region " << name << endl;
surfBaffleRegions[pI] = true;
}
}
// Extracting and writing a extendedFeatureEdgeMesh
extendedFeatureEdgeMesh feMesh
(
newSet,
runTime,
sFeatFileName + ".extendedFeatureEdgeMesh"
sFeatFileName + ".extendedFeatureEdgeMesh",
surfBaffleRegions
);

6
src/OpenFOAM/db/error/error.H
View File

@ -28,11 +28,11 @@ Description
Class to handle errors and exceptions in a simple, consistent stream-based
manner.
The error class is globaly instantiated with a title string. Errors,
The error class is globally instantiated with a title string. Errors,
messages and other data are piped to the messageStream class in the
standard manner. Manipulators are supplied for exit and abort which may
terminate the program or throw an exception depending of if the exception
handling has beed switched on (off by default).
terminate the program or throw an exception depending on whether the
exception handling has been switched on (off by default).
Usage
\code

6
src/OpenFOAM/memory/autoPtr/autoPtrI.H
View File

@ -52,10 +52,14 @@ inline Foam::autoPtr<T>::autoPtr(const autoPtr<T>& ap, const bool reUse)
ptr_ = ap.ptr_;
ap.ptr_ = 0;
}
else
else if (ap.valid())
{
ptr_ = ap().clone().ptr();
}
else
{
ptr_ = NULL;
}
}

32
src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
View File

@ -32,7 +32,9 @@ License
#include "polyTopoChange.H"
#include "globalIndex.H"
#include "PackedBoolList.H"
#include "faceZone.H"
#include "pointSet.H"
#include "faceSet.H"
#include "cellSet.H"
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
@ -42,6 +44,26 @@ defineTypeNameAndDebug(polyMeshFilter, 0);
}
void Foam::polyMeshFilter::updateSets(const mapPolyMesh& map)
{
updateSets<pointSet>(map);
updateSets<faceSet>(map);
updateSets<cellSet>(map);
}
void Foam::polyMeshFilter::copySets
(
const polyMesh& oldMesh,
const polyMesh& newMesh
)
{
copySets<pointSet>(oldMesh, newMesh);
copySets<faceSet>(oldMesh, newMesh);
copySets<cellSet>(oldMesh, newMesh);
}
Foam::autoPtr<Foam::fvMesh> Foam::polyMeshFilter::copyMesh(const fvMesh& mesh)
{
polyTopoChange originalMeshToNewMesh(mesh);
@ -72,6 +94,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::polyMeshFilter::copyMesh(const fvMesh& mesh)
meshCopy().movePoints(map.preMotionPoints());
}
copySets(mesh, meshCopy());
return meshCopy;
}
@ -359,6 +383,7 @@ Foam::label Foam::polyMeshFilter::filterFaces
{
newMesh.movePoints(newMap.preMotionPoints());
}
updateSets(newMap);
updatePointPriorities(newMesh, newMap.pointMap());
@ -475,6 +500,7 @@ Foam::label Foam::polyMeshFilter::filterEdges
{
newMesh.movePoints(newMap.preMotionPoints());
}
updateSets(newMap);
// Synchronise the factors
mapOldMeshEdgeFieldToNewMesh
@ -926,14 +952,14 @@ Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
}
Foam::label Foam::polyMeshFilter::filter(const faceZone& fZone)
Foam::label Foam::polyMeshFilter::filter(const faceSet& fSet)
{
minEdgeLen_.resize(mesh_.nEdges(), minLen());
faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor());
forAll(faceFilterFactor_, fI)
{
if (fZone.whichFace(fI) == -1)
if (!fSet.found(fI))
{
faceFilterFactor_[fI] = -1;
}

38
src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
View File

@ -25,13 +25,14 @@ Class
Foam::polyMeshFilter
Description
Filter the edges and faces of a polyMesh whilst satisfying the given mesh
Remove the edges and faces of a polyMesh whilst satisfying the given mesh
quality criteria.
Works on a copy of the mesh.
SourceFiles
polyMeshFilter.C
polyMeshFilterTemplates.C
\*---------------------------------------------------------------------------*/
@ -53,7 +54,7 @@ namespace Foam
class polyMesh;
class fvMesh;
class PackedBoolList;
class faceZone;
class faceSet;
/*---------------------------------------------------------------------------*\
Class polyMeshFilter Declaration
@ -71,8 +72,12 @@ class polyMeshFilter
//- Copy of the original mesh to perform the filtering on
autoPtr<fvMesh> newMeshPtr_;
//-
//- Original point priorities. If a point has a higher priority than
// another point then the edge between them collapses towards the
// point with the higher priority. (e.g. 2 is a higher priority than 1)
labelList originalPointPriority_;
//- Point priority associated with the new mesh
autoPtr<labelList> pointPriority_;
//- The minimum edge length for each edge
@ -84,6 +89,12 @@ class polyMeshFilter
// Private Member Functions
template<typename T>
static void updateSets(const mapPolyMesh& map);
template<typename T>
static void copySets(const polyMesh& oldMesh, const polyMesh& newMesh);
label filterFacesLoop(const label nOriginalBadFaces);
label filterFaces
@ -209,6 +220,7 @@ public:
// mesh has actually been filtered.
const autoPtr<fvMesh>& filteredMesh() const;
//- Return the new pointPriority list.
const autoPtr<labelList>& pointPriority() const;
@ -217,11 +229,21 @@ public:
//- Return a copy of an fvMesh
static autoPtr<fvMesh> copyMesh(const fvMesh& mesh);
//- Copy loaded topoSets from the old mesh to the new mesh
static void copySets
(
const polyMesh& oldMesh,
const polyMesh& newMesh
);
//- Update the loaded topoSets
static void updateSets(const mapPolyMesh& map);
//- Filter edges and faces
label filter(const label nOriginalBadFaces);
//- Filter all faces that are in the face zone indirectPatchFaces
label filter(const faceZone& fZone);
//- Filter all faces that are in the face set
label filter(const faceSet& fSet);
//- Filter edges only.
label filterEdges(const label nOriginalBadFaces);
@ -234,6 +256,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "polyMeshFilterTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

1
src/dynamicMesh/polyMeshFilter/polyMeshFilterSettings.H
View File

@ -29,6 +29,7 @@ Description
SourceFiles
polyMeshFilterSettings.C
polyMeshFilterSettingsI.H
\*---------------------------------------------------------------------------*/

105
src/dynamicMesh/polyMeshFilter/polyMeshFilterTemplates.C Normal file
View File

@ -0,0 +1,105 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "polyMeshFilter.H"
#include "polyMesh.H"
#include "mapPolyMesh.H"
#include "IOobjectList.H"
// * * * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * //
template<typename SetType>
void Foam::polyMeshFilter::updateSets(const mapPolyMesh& map)
{
HashTable<const SetType*> sets =
map.mesh().objectRegistry::lookupClass<const SetType>();
forAllIter(typename HashTable<const SetType*>, sets, iter)
{
SetType& set = const_cast<SetType&>(*iter());
set.updateMesh(map);
set.sync(map.mesh());
}
IOobjectList Objects
(
map.mesh().time(),
map.mesh().facesInstance(),
"polyMesh/sets"
);
IOobjectList fileSets(Objects.lookupClass(SetType::typeName));
forAllConstIter(IOobjectList, fileSets, iter)
{
if (!sets.found(iter.key()))
{
// Not in memory. Load it.
SetType set(*iter());
set.updateMesh(map);
set.write();
}
}
}
template<typename SetType>
void Foam::polyMeshFilter::copySets
(
const polyMesh& oldMesh,
const polyMesh& newMesh
)
{
HashTable<const SetType*> sets =
oldMesh.objectRegistry::lookupClass<const SetType>();
forAllConstIter(typename HashTable<const SetType*>, sets, iter)
{
const SetType& set = *iter();
if (newMesh.objectRegistry::foundObject<SetType>(set.name()))
{
const SetType& origSet =
newMesh.objectRegistry::lookupObject<SetType>(set.name());
const_cast<SetType&>(origSet) = set;
const_cast<SetType&>(origSet).sync(newMesh);
}
else
{
SetType* newSet
(
new SetType(newMesh, set.name(), set, set.writeOpt())
);
newSet->store();
newSet->sync(newMesh);
}
}
}
// ************************************************************************* //

32
src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
View File

@ -263,7 +263,8 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
(
const surfaceFeatures& sFeat,
const objectRegistry& obr,
const fileName& sFeatFileName
const fileName& sFeatFileName,
const boolList& surfBaffleRegions
)
:
regIOobject
@ -337,7 +338,12 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
// Check to see if the points have been already used
if (faceMap[eFI] == -1)
{
normalVolumeTypes_[nAdded++] = INSIDE;
normalVolumeTypes_[nAdded++] =
(
surfBaffleRegions[surf[eFI].region()]
? BOTH
: INSIDE
);
faceMap[eFI] = nAdded - 1;
}
@ -492,7 +498,7 @@ Foam::extendedFeatureEdgeMesh::classifyEdge
const List<vector>& norms,
const labelList& edNorms,
const vector& fC0tofC1
) const
)
{
label nEdNorms = edNorms.size();
@ -502,8 +508,8 @@ Foam::extendedFeatureEdgeMesh::classifyEdge
}
else if (nEdNorms == 2)
{
const vector n0(norms[edNorms[0]]);
const vector n1(norms[edNorms[1]]);
const vector& n0(norms[edNorms[0]]);
const vector& n1(norms[edNorms[1]]);
if ((n0 & n1) > cosNormalAngleTol_)
{
@ -1342,6 +1348,22 @@ void Foam::extendedFeatureEdgeMesh::writeObj
meshTools::writeOBJ(regionStr, points()[e[1]]); verti++;
regionStr << "l " << verti-1 << ' ' << verti << endl;
}
OFstream edgeDirsStr(prefix + "_edgeDirections.obj");
Info<< "Writing edge directions to " << edgeDirsStr.name() << endl;
forAll(edgeDirections_, i)
{
const vector& eVec = edgeDirections_[i];
const edge& e = edges()[i];
meshTools::writeOBJ
(
edgeDirsStr,
points()[e.start()],
eVec + points()[e.start()]
);
}
}

31
src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
View File

@ -119,14 +119,14 @@ public:
static const Foam::NamedEnum<sideVolumeType, 4> sideVolumeTypeNames_;
//- Angular closeness tolerance for treating normals as the same
static scalar cosNormalAngleTol_;
private:
// Static data
//- Angular closeness tolerance for treating normals as the same
static scalar cosNormalAngleTol_;
//- Index of the start of the convex feature points - static as 0
static label convexStart_;
@ -202,15 +202,6 @@ private:
// data for edges and normals.
pointStatus classifyFeaturePoint(label ptI) const;
//- Classify the type of feature edge. Requires face centre 0 to face
// centre 1 vector to distinguish internal from external
edgeStatus classifyEdge
(
const List<vector>& norms,
const labelList& edNorms,
const vector& fC0tofC1
) const;
template<class Patch>
void sortPointsAndEdges
(
@ -259,7 +250,8 @@ public:
(
const surfaceFeatures& sFeat,
const objectRegistry& obr,
const fileName& sFeatFileName
const fileName& sFeatFileName,
const boolList& surfBaffleRegions
);
//- Construct from PrimitivePatch
@ -434,6 +426,9 @@ public:
//- Return the edgeStatus of a specified edge
inline edgeStatus getEdgeStatus(label edgeI) const;
//- Return the baffle faces of a specified edge
inline PackedList<2> edgeBaffles(label edgeI) const;
//- Demand driven construction of octree for feature points
const indexedOctree<treeDataPoint>& pointTree() const;
@ -468,6 +463,16 @@ public:
friend Istream& operator>>(Istream& is, sideVolumeType& vt);
friend Ostream& operator<<(Ostream& os, const sideVolumeType& vt);
//- Classify the type of feature edge. Requires face centre 0 to face
// centre 1 vector to distinguish internal from external
static edgeStatus classifyEdge
(
const List<vector>& norms,
const labelList& edNorms,
const vector& fC0tofC1
);
};

23
src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
View File

@ -264,4 +264,27 @@ Foam::extendedFeatureEdgeMesh::getEdgeStatus(label edgeI) const
}
inline Foam::PackedList<2> Foam::extendedFeatureEdgeMesh::edgeBaffles
(
label edgeI
) const
{
const labelList& eNormals = edgeNormals_[edgeI];
DynamicList<label> edgeBaffles(eNormals.size());
forAll(eNormals, enI)
{
const label normI = eNormals[enI];
if (normalVolumeTypes_[normI])
{
edgeBaffles.append(normI);
}
}
return PackedList<2>(edgeBaffles);
}
// ************************************************************************* //

18
src/fvAgglomerationMethods/Allwmake
View File

@ -12,4 +12,22 @@ fi
wmake $makeType pairPatchAgglomeration
## get SCOTCH_VERSION, SCOTCH_ARCH_PATH
#if settings=`$WM_PROJECT_DIR/bin/foamEtcFile config/scotch.sh`
#then
# . $settings
# echo "using SCOTCH_ARCH_PATH=$SCOTCH_ARCH_PATH"
#else
# echo
# echo "Error: no config/scotch.sh settings"
# echo
#fi
#
#if [ -n "$SCOTCH_ARCH_PATH" ]
#then
# wmake $makeType scotchGamgAgglomeration
#fi
# ----------------------------------------------------------------- end-of-file

2
src/mesh/autoMesh/Make/files
View File

@ -17,6 +17,8 @@ $(autoHexMesh)/meshRefinement/meshRefinement.C
$(autoHexMesh)/meshRefinement/meshRefinementMerge.C
$(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C
$(autoHexMesh)/meshRefinement/meshRefinementRefine.C
$(autoHexMesh)/meshRefinement/patchFaceOrientation.C
$(autoHexMesh)/refinementFeatures/refinementFeatures.C
$(autoHexMesh)/refinementSurfaces/surfaceZonesInfo.C
$(autoHexMesh)/refinementSurfaces/refinementSurfaces.C

10
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
View File

@ -112,8 +112,6 @@ const
if (localCellI != -1)
{
Pout<< "Found point " << keepPoint << " in cell " << localCellI
<< " on processor " << Pstream::myProcNo() << endl;
globalCellI = globalCells.toGlobal(localCellI);
}
@ -130,6 +128,14 @@ const
<< exit(FatalError);
}
label procI = globalCells.whichProcID(globalCellI);
label procCellI = globalCells.toLocal(procI, globalCellI);
Info<< "Found point " << keepPoint << " in cell " << procCellI
<< " on processor " << procI << endl;
if (globalCells.isLocal(globalCellI))
{
cellLabels[i] = localCellI;

16
src/meshTools/regionSplit/regionSplit.C
View File

@ -379,6 +379,7 @@ Foam::label Foam::regionSplit::calcLocalRegionSplit
Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
(
const bool doGlobalRegions,
const boolList& blockedFace,
const List<labelPair>& explicitConnections,
@ -395,7 +396,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
cellRegion
);
if (!Pstream::parRun())
if (!doGlobalRegions)
{
return autoPtr<globalIndex>(new globalIndex(nLocalRegions));
}
@ -422,7 +423,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
// (this will create gaps in the global region list so they will get
// merged later on)
while (Pstream::parRun())
while (true)
{
if (debug)
{
@ -690,13 +691,14 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::regionSplit::regionSplit(const polyMesh& mesh)
Foam::regionSplit::regionSplit(const polyMesh& mesh, const bool doGlobalRegions)
:
MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
labelList(mesh.nCells(), -1)
{
globalNumberingPtr_ = calcRegionSplit
(
doGlobalRegions, //do global regions
boolList(0, false), //blockedFaces
List<labelPair>(0), //explicitConnections,
*this
@ -707,7 +709,8 @@ Foam::regionSplit::regionSplit(const polyMesh& mesh)
Foam::regionSplit::regionSplit
(
const polyMesh& mesh,
const boolList& blockedFace
const boolList& blockedFace,
const bool doGlobalRegions
)
:
MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
@ -715,6 +718,7 @@ Foam::regionSplit::regionSplit
{
globalNumberingPtr_ = calcRegionSplit
(
doGlobalRegions,
blockedFace, //blockedFaces
List<labelPair>(0), //explicitConnections,
*this
@ -726,7 +730,8 @@ Foam::regionSplit::regionSplit
(
const polyMesh& mesh,
const boolList& blockedFace,
const List<labelPair>& explicitConnections
const List<labelPair>& explicitConnections,
const bool doGlobalRegions
)
:
MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
@ -734,6 +739,7 @@ Foam::regionSplit::regionSplit
{
globalNumberingPtr_ = calcRegionSplit
(
doGlobalRegions,
blockedFace, //blockedFaces
explicitConnections, //explicitConnections,
*this

26
src/meshTools/regionSplit/regionSplit.H
View File

@ -87,6 +87,9 @@ Description
proc0 | proc1 | proc2
Can optionally keep all regions local to the processor.
SourceFiles
regionSplit.C
@ -155,6 +158,7 @@ class regionSplit
//- Calculate global region split. Return globalIndex.
autoPtr<globalIndex> calcRegionSplit
(
const bool doGlobalRegions,
const boolList& blockedFace,
const List<labelPair>& explicitConnections,
labelList& cellRegion
@ -170,11 +174,20 @@ public:
// Constructors
//- Construct from mesh
regionSplit(const polyMesh&);
regionSplit
(
const polyMesh&,
const bool doGlobalRegions = Pstream::parRun()
);
//- Construct from mesh and whether face is blocked
// NOTE: blockedFace has to be consistent across coupled faces!
regionSplit(const polyMesh&, const boolList& blockedFace);
regionSplit
(
const polyMesh&,
const boolList& blockedFace,
const bool doGlobalRegions = Pstream::parRun()
);
//- Construct from mesh and whether face is blocked. Additional explicit
// connections between normal boundary faces.
@ -183,7 +196,8 @@ public:
(
const polyMesh&,
const boolList& blockedFace,
const List<labelPair>&
const List<labelPair>&,
const bool doGlobalRegions = Pstream::parRun()
);
@ -195,6 +209,12 @@ public:
return globalNumberingPtr_();
}
//- Return local number of regions
label nLocalRegions() const
{
return globalNumbering().localSize(Pstream::myProcNo());
}
//- Return total number of regions
label nRegions() const
{

136
src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
View File

@ -117,19 +117,16 @@ void Foam::edgeIntersections::intersectEdges
<< " triangles ..." << endl;
}
// Construct octree.
const indexedOctree<treeDataTriSurface>& tree = querySurf2.tree();
label nHits = 0;
pointField start(edgeLabels.size());
pointField end(edgeLabels.size());
vectorField edgeDirs(edgeLabels.size());
// Go through all edges, calculate intersections
forAll(edgeLabels, i)
{
label edgeI = edgeLabels[i];
if (debug && (i % 1000 == 0))
if (debug)// && (i % 1000 == 0))
{
Pout<< "Intersecting edge " << edgeI << " with surface" << endl;
}
@ -140,85 +137,81 @@ void Foam::edgeIntersections::intersectEdges
const point& pEnd = points1[meshPoints[e.end()]];
const vector eVec(pEnd - pStart);
const scalar eMag = mag(eVec);
const vector n(eVec/(eMag + VSMALL));
const vector n(eVec/(mag(eVec) + VSMALL));
// Smallish length for intersection calculation errors.
const point tolVec = 1e-6*eVec;
// Start tracking somewhat before pStart and upto somewhat after p1.
// Start tracking somewhat before pStart and up to somewhat after p1.
// Note that tolerances here are smaller than those used to classify
// hit below.
// This will cause this hit to be marked as degenerate and resolved
// later on.
point p0 = pStart - 0.5*surf1PointTol[e[0]]*n;
const point p1 = pEnd + 0.5*surf1PointTol[e[1]]*n;
const scalar maxS = mag(p1 - pStart);
start[i] = pStart - 0.5*surf1PointTol[e[0]]*n;
end[i] = pEnd + 0.5*surf1PointTol[e[1]]*n;
// Get all intersections of the edge with the surface
edgeDirs[i] = n;
}
DynamicList<pointIndexHit> currentIntersections(100);
DynamicList<label> currentIntersectionTypes(100);
List<List<pointIndexHit> > edgeIntersections;
querySurf2.findLineAll
(
start,
end,
edgeIntersections
);
while (true)
label nHits = 0;
// Classify the hits
forAll(edgeLabels, i)
{
const label edgeI = edgeLabels[i];
labelList& intersectionTypes = classification_[edgeI];
intersectionTypes.setSize(edgeIntersections[i].size(), -1);
this->operator[](edgeI).transfer(edgeIntersections[i]);
forAll(intersectionTypes, hitI)
{
pointIndexHit pHit = tree.findLine(p0, p1);
const pointIndexHit& pHit = this->operator[](edgeI)[hitI];
label& hitType = intersectionTypes[hitI];
if (pHit.hit())
if (!pHit.hit())
{
nHits++;
currentIntersections.append(pHit);
// Classify point on surface1 edge.
label edgeEnd = -1;
if (mag(pHit.hitPoint() - pStart) < surf1PointTol[e[0]])
{
edgeEnd = 0;
}
else if (mag(pHit.hitPoint() - pEnd) < surf1PointTol[e[1]])
{
edgeEnd = 1;
}
else if (mag(n & normals2[pHit.index()]) < alignedCos_)
{
Pout<< "Flat angle edge:" << edgeI
<< " face:" << pHit.index()
<< " cos:" << mag(n & normals2[pHit.index()])
<< endl;
edgeEnd = 2;
}
currentIntersectionTypes.append(edgeEnd);
if (edgeEnd == 1)
{
// Close to end
break;
}
else
{
// Continue tracking. Shift by a small amount.
p0 = pHit.hitPoint() + tolVec;
if (((p0-pStart) & n) >= maxS)
{
break;
}
}
continue;
}
else
const edge& e = surf1.edges()[edgeI];
// Classify point on surface1 edge.
if (mag(pHit.hitPoint() - start[i]) < surf1PointTol[e[0]])
{
// No hit.
break;
// Intersection is close to edge start
hitType = 0;
}
else if (mag(pHit.hitPoint() - end[i]) < surf1PointTol[e[1]])
{
// Intersection is close to edge end
hitType = 1;
}
else if (mag(edgeDirs[i] & normals2[pHit.index()]) < alignedCos_)
{
// Edge is almost coplanar with the face
Pout<< "Flat angle edge:" << edgeI
<< " face:" << pHit.index()
<< " cos:" << mag(edgeDirs[i] & normals2[pHit.index()])
<< endl;
hitType = 2;
}
if (debug)
{
Info<< " hit " << pHit << " classify = " << hitType << endl;
}
nHits++;
}
// Done current edge. Transfer all data into *this
operator[](edgeI).transfer(currentIntersections);
classification_[edgeI].transfer(currentIntersectionTypes);
}
if (debug)
@ -226,7 +219,6 @@ void Foam::edgeIntersections::intersectEdges
Pout<< "Found " << nHits << " intersections of edges with surface ..."
<< endl;
}
}

11
src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
View File

@ -202,7 +202,6 @@ void Foam::surfaceIntersection::storeIntersection
DynamicList<point>& allCutPoints
)
{
forAll(facesA, facesAI)
{
label faceA = facesA[facesAI];
@ -969,10 +968,10 @@ Foam::surfaceIntersection::surfaceIntersection
if (!usedPoints.found(pointI))
{
FatalErrorIn("surfaceIntersection::surfaceIntersection")
WarningIn("surfaceIntersection::surfaceIntersection")
<< "Problem: cut point:" << pointI
<< " coord:" << cutPoints_[pointI]
<< " not used by any edge" << abort(FatalError);
<< " not used by any edge" << endl;
}
}
}
@ -1147,6 +1146,12 @@ const Foam::edgeList& Foam::surfaceIntersection::cutEdges() const
}
const Foam::labelPairLookup& Foam::surfaceIntersection::facePairToVertex() const
{
return facePairToVertex_;
}
const Foam::labelPairLookup& Foam::surfaceIntersection::facePairToEdge() const
{
return facePairToEdge_;

2
src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H
View File

@ -273,7 +273,7 @@ public:
const edgeList& cutEdges() const;
//const labelPairLookup& facePairToVertex() const;
const labelPairLookup& facePairToVertex() const;
const labelPairLookup& facePairToEdge() const;

60
src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
View File

@ -167,7 +167,7 @@ void Foam::decompositionMethod::calcCellCells
(
const polyMesh& mesh,
const labelList& agglom,
const label nCoarse,
const label nLocalCoarse,
const bool parallel,
CompactListList<label>& cellCells
)
@ -182,7 +182,7 @@ void Foam::decompositionMethod::calcCellCells
globalIndex globalAgglom
(
nCoarse,
nLocalCoarse,
Pstream::msgType(),
Pstream::worldComm,
parallel
@ -224,7 +224,7 @@ void Foam::decompositionMethod::calcCellCells
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Number of faces per coarse cell
labelList nFacesPerCell(nCoarse, 0);
labelList nFacesPerCell(nLocalCoarse, 0);
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
@ -374,7 +374,7 @@ void Foam::decompositionMethod::calcCellCells
// const boolList& blockedFace,
// const List<labelPair>& explicitConnections,
// const labelList& agglom,
// const label nCoarse,
// const label nLocalCoarse,
// const bool parallel,
// CompactListList<label>& cellCells
//)
@ -389,7 +389,7 @@ void Foam::decompositionMethod::calcCellCells
//
// globalIndex globalAgglom
// (
// nCoarse,
// nLocalCoarse,
// Pstream::msgType(),
// Pstream::worldComm,
// parallel
@ -431,7 +431,7 @@ void Foam::decompositionMethod::calcCellCells
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// // Number of faces per coarse cell
// labelList nFacesPerCell(nCoarse, 0);
// labelList nFacesPerCell(nLocalCoarse, 0);
//
// // 1. Internal faces
// for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
@ -764,6 +764,24 @@ Foam::labelList Foam::decompositionMethod::decompose
// Any weights specified?
label nWeights = returnReduce(cellWeights.size(), sumOp<label>());
if (nWeights > 0 && cellWeights.size() != mesh.nCells())
{
FatalErrorIn
(
"decompositionMethod::decompose\n"
"(\n"
" const polyMesh&,\n"
" const scalarField&,\n"
" const boolList&,\n"
" const PtrList<labelList>&,\n"
" const labelList&,\n"
" const List<labelPair>&\n"
) << "Number of weights " << cellWeights.size()
<< " differs from number of cells " << mesh.nCells()
<< exit(FatalError);
}
// Any processor sets?
label nProcSets = 0;
forAll(specifiedProcessorFaces, setI)
@ -828,10 +846,18 @@ Foam::labelList Foam::decompositionMethod::decompose
<< nProcSets << endl << endl;
}
// Determine global regions, separated by blockedFaces
regionSplit globalRegion(mesh, blockedFace, explicitConnections);
// Determine local regions, separated by blockedFaces
regionSplit localRegion(mesh, blockedFace, explicitConnections, false);
if (debug)
{
Info<< "Constrained decomposition:" << endl
<< " split into " << localRegion.nLocalRegions()
<< " regions."
<< endl;
}
// Determine region cell centres
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -840,11 +866,11 @@ Foam::labelList Foam::decompositionMethod::decompose
// somewhere in the middle of the domain which might not be anywhere
// near any of the cells.
pointField regionCentres(globalRegion.nRegions(), point::max);
pointField regionCentres(localRegion.nLocalRegions(), point::max);
forAll(globalRegion, cellI)
forAll(localRegion, cellI)
{
label regionI = globalRegion[cellI];
label regionI = localRegion[cellI];
if (regionCentres[regionI] == point::max)
{
@ -855,22 +881,22 @@ Foam::labelList Foam::decompositionMethod::decompose
// Do decomposition on agglomeration
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalarField regionWeights(globalRegion.nRegions(), 0);
scalarField regionWeights(localRegion.nLocalRegions(), 0);
if (nWeights > 0)
{
forAll(globalRegion, cellI)
forAll(localRegion, cellI)
{
label regionI = globalRegion[cellI];
label regionI = localRegion[cellI];
regionWeights[regionI] += cellWeights[cellI];
}
}
else
{
forAll(globalRegion, cellI)
forAll(localRegion, cellI)
{
label regionI = globalRegion[cellI];
label regionI = localRegion[cellI];
regionWeights[regionI] += 1.0;
}
@ -879,7 +905,7 @@ Foam::labelList Foam::decompositionMethod::decompose
finalDecomp = decompose
(
mesh,
globalRegion,
localRegion,
regionCentres,
regionWeights
);

27
src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
View File

@ -219,7 +219,7 @@ public:
// Other
//- Helper: determine (local or global) cellCells from mesh
// agglomeration.
// agglomeration. Agglomeration is local to the processor.
// local : connections are in local indices. Coupled across
// cyclics but not processor patches.
// global : connections are in global indices. Coupled across
@ -228,34 +228,11 @@ public:
(
const polyMesh& mesh,
const labelList& agglom,
const label nCoarse,
const label nLocalCoarse,
const bool global,
CompactListList<label>& cellCells
);
//- Helper: determine (local or global) cellCells from mesh
// agglomeration and additional specification:
// - any additional connections between non-coupled internal
// or boundary faces.
// - any internal or coupled faces (or additional connections)
// are blocked
//
// local : connections are in local indices. Coupled across
// cyclics but not processor patches.
// global : connections are in global indices. Coupled across
// cyclics and processor patches.
//static void calcCellCells
//(
// const polyMesh& mesh,
// const boolList& blockedFace,
// const List<labelPair>& explicitConnections,
// const labelList& agglom,
// const label nCoarse,
// const bool global,
// CompactListList<label>& cellCells
//);
//- Helper: extract constraints:
// blockedface: existing faces where owner and neighbour on same
// proc

4
src/thermophysicalModels/laminarFlameSpeed/Gulders/Gulders.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -217,7 +217,7 @@ Foam::laminarFlameSpeedModels::Gulders::operator()() const
dimensionedScalar
(
psiuReactionThermo_.lookup("stoichiometricAirFuelMassRatio")
)*ft/((1 + SMALL) - ft)
)*ft/max(1 - ft, SMALL)
);
}
else

6
tutorials/mesh/foamyHexMesh/blob/Allrun
View File

@ -8,7 +8,13 @@ cd ${0%/*} || exit 1 # run from this directory
cp $FOAM_TUTORIALS/resources/geometry/blob.stl.gz constant/triSurface/
runApplication foamyHexMesh
runApplication collapseEdges -latestTime -collapseFaces
mv log.collapseEdges log.collapseFaces
runApplication collapseEdges -latestTime -collapseFaceSet indirectPatchFaces
mv log.collapseEdges log.collapseFaceSet
runApplication checkMesh -latestTime -allGeometry -allTopology

6
tutorials/mesh/foamyHexMesh/blob/Allrun-parallel
View File

@ -14,7 +14,13 @@ runApplication blockMesh -region backgroundMeshDecomposition
runApplication decomposePar -region backgroundMeshDecomposition
runParallel foamyHexMesh $nProc
runParallel collapseEdges $nProc -latestTime -collapseFaces
mv log.collapseEdges log.collapseFaces
runParallel collapseEdges $nProc -latestTime -collapseFaceSet indirectPatchFaces
mv log.collapseEdges log.collapseFaceSet
runParallel checkMesh $nProc -latestTime -allTopology -allGeometry
runApplication reconstructParMesh -latestTime

6
tutorials/mesh/foamyHexMesh/flange/Allrun
View File

@ -8,7 +8,13 @@ cd ${0%/*} || exit 1 # run from this directory
cp $FOAM_TUTORIALS/resources/geometry/flange.stl.gz constant/triSurface/
runApplication foamyHexMesh
runApplication collapseEdges -latestTime -collapseFaces
mv log.collapseEdges log.collapseFaces
runApplication collapseEdges -latestTime -collapseFaceSet indirectPatchFaces
mv log.collapseEdges log.collapseFaceSet
runApplication checkMesh -latestTime -allGeometry -allTopology

6
tutorials/mesh/foamyHexMesh/flange/Allrun-parallel
View File

@ -15,7 +15,13 @@ runApplication blockMesh -region backgroundMeshDecomposition
runApplication decomposePar -region backgroundMeshDecomposition
runParallel foamyHexMesh $nProc
runParallel collapseEdges $nProc -latestTime -collapseFaces
mv log.collapseEdges log.collapseFaces
runParallel collapseEdges $nProc -latestTime -collapseFaceSet indirectPatchFaces
mv log.collapseEdges log.collapseFaceSet
runParallel checkMesh $nProc -latestTime -allTopology -allGeometry
runApplication reconstructParMesh -latestTime

10
tutorials/mesh/foamyHexMesh/mixerVessel/system/foamyHexMeshDict
View File

@ -59,9 +59,9 @@ surfaceConformation
extendedFeatureEdgeMesh "vessel.extendedFeatureEdgeMesh";
regions
{
patch0 {}
patch1 {}
patch2 {}
patch3 {}
}
}
@ -118,7 +118,7 @@ motionControl
minimumCellSizeCoeff 0.1;
maxRefinementIterations 0;
maxRefinementIterations 1;
maxSmoothingIterations 100;
@ -142,7 +142,7 @@ motionControl
regions
{
patch1
patch0
{
surfaceCellSizeFunction uniformValue;
uniformValueCoeffs
@ -154,7 +154,7 @@ motionControl
uniformCoeffs{}
}
patch2
patch1
{
priority 2;
surfaceCellSizeFunction uniformValue;
@ -170,7 +170,7 @@ motionControl
}
}
patch3
patch2
{
priority 2;
surfaceCellSizeFunction uniformValue;

2
tutorials/mesh/foamyHexMesh/simpleShapes/Allclean
View File

@ -15,6 +15,8 @@ rm -r snapToSurface?.obj tetsToSnapTo.obj > /dev/null 2>&1
rm domain coneAndSphere > /dev/null 2>&1
rm -rf 0/
cleanCase
# ----------------------------------------------------------------- end-of-file

2
tutorials/mesh/foamyHexMesh/simpleShapes/Allrun
View File

@ -22,5 +22,7 @@ runApplication surfaceBooleanFeatures intersection \
runApplication foamyHexMesh
runApplication collapseEdges -latestTime -collapseFaceSet indirectPatchFaces
# ----------------------------------------------------------------- end-of-file

1
wmake/rules/General/CGAL
View File

@ -1,5 +1,4 @@
CGAL_INC = \
-Wno-old-style-cast \
-I$(CGAL_ARCH_PATH)/include \
-I$(MPFR_ARCH_PATH)/include \
-I$(GMP_ARCH_PATH)/include \

4
wmake/rules/linux64Clang/c++
View File

@ -3,6 +3,10 @@
# -Woverloaded-virtual may produce spurious warnings, disable for now
c++WARN = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -Wno-overloaded-virtual -Wno-unused-comparison
# Suppress CGAL warnings
c++CGALWARN = -Wno-c++11-extensions -Wno-sometimes-uninitialized -Wno-mismatched-tags
CC = clang++ -m64
include $(RULES)/c++$(WM_COMPILE_OPTION)

4
wmake/rules/linuxClang/c++
View File

@ -3,6 +3,10 @@
# -Woverloaded-virtual may produce spurious warnings, disable for now
c++WARN = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -Wno-overloaded-virtual -Wno-unused-comparison
# Suppress CGAL warnings
c++CGALWARN = -Wno-c++11-extensions -Wno-sometimes-uninitialized -Wno-mismatched-tags
CC = clang++ -m32
include $(RULES)/c++$(WM_COMPILE_OPTION)