diff --git a/applications/utilities/mesh/generation/cvMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/Make/options
index a0a776d977b1f754a01560aa9e671c9b108dd7c7..5efdb125725138a1df13a4eae74625e3a0fc8cc0 100644
--- a/applications/utilities/mesh/generation/cvMesh/Make/options
+++ b/applications/utilities/mesh/generation/cvMesh/Make/options
@@ -20,7 +20,8 @@ EXE_LIBS = \
$(CGAL_LIBS) \
-lconformalVoronoiMesh \
-lmeshTools \
- -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \
+ -ldecompositionMethods \
+ -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
-ledgeMesh \
-ltriSurface \
-ldynamicMesh
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
index 1a1902b559ddda2a3ff6c37acca1316f89b6b7ea..2d2d1a1178abede674023f52331496cd82e291a6 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
@@ -19,7 +19,6 @@ EXE_INC = \
EXE_LIBS = \
-lmeshTools \
- -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \
-ledgeMesh \
-ltriSurface \
-ldynamicMesh
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C
index d40427abf7e5fa9bb6f60945b9b66a57a13f3133..788c1d69c6852f1c51e71143c59f5cb1b07f7303 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C
@@ -63,7 +63,7 @@ Foam::fieldFromFile::fieldFromFile
Foam::triSurfaceScalarField Foam::fieldFromFile::load()
{
- Info<< "Loading: " << fileName_ << endl;
+ Info<< indent << "Loading: " << fileName_ << endl;
triSurfaceScalarField surfaceCellSize
(
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C
index 5e2d032d91e1c66979f1c646bc836845620dd4d7..fc91ed279b04a7979eface6cf063fec1ab5b36c4 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C
@@ -82,9 +82,9 @@ Foam::nonUniformField::nonUniformField
Info<< decrIndent;
Info<< indent << "Cell size field statistics:" << nl
- << indent << " Minimum: " << min(surfaceCellSize_).value() << nl
- << indent << " Average: " << average(surfaceCellSize_).value() << nl
- << indent << " Maximum: " << max(surfaceCellSize_).value() << endl;
+ << indent << " Minimum: " << min(surfaceCellSize_).value() << nl
+ << indent << " Average: " << average(surfaceCellSize_).value() << nl
+ << indent << " Maximum: " << max(surfaceCellSize_).value() << endl;
Info<< decrIndent;
}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/surfaceCellSizeFunction/surfaceCellSizeFunction.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/surfaceCellSizeFunction/surfaceCellSizeFunction.C
index 7022511c4c2dbdf594ecc6666cf4c6fa5712a833..04d9af824f0980840a2e1846f9144efb404b891f 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/surfaceCellSizeFunction/surfaceCellSizeFunction.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/surfaceCellSizeFunction/surfaceCellSizeFunction.C
@@ -46,7 +46,10 @@ Foam::surfaceCellSizeFunction::surfaceCellSizeFunction
dictionary(surfaceCellSizeFunctionDict),
surface_(surface),
coeffsDict_(subDict(type + "Coeffs")),
- refinementFactor_(readScalar(lookup("refinementFactor")))
+ refinementFactor_
+ (
+ lookupOrDefault<scalar>("refinementFactor", 1.0)
+ )
{}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 09d27dfa3620f9e110a7704e2a56d06ae88c6069..de23690ea43b26ca788d49724cb42e2d17a6e192 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -44,6 +44,142 @@ const Foam::scalar Foam::conformalVoronoiMesh::tolParallel = 1e-3;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+Foam::scalar Foam::conformalVoronoiMesh::requiredSize
+(
+ const Foam::point& pt
+) const
+{
+ pointIndexHit surfHit;
+ label hitSurface;
+
+ DynamicList<scalar> cellSizeHits;
+
+ geometryToConformTo_.findSurfaceNearest
+ (
+ pt,
+ sqr(GREAT),
+ surfHit,
+ hitSurface
+ );
+
+ if (!surfHit.hit())
+ {
+ FatalErrorIn
+ (
+ "Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment"
+ )
+ << "findSurfaceNearest did not find a hit across the surfaces."
+ << exit(FatalError) << endl;
+ }
+
+ cellSizeHits.append(cellSizeControl().cellSize(pt));
+
+ // Primary alignment
+
+ vectorField norm(1);
+
+ allGeometry_[hitSurface].getNormal
+ (
+ List<pointIndexHit>(1, surfHit),
+ norm
+ );
+
+ const vector np = norm[0];
+
+ // Generate equally spaced 'spokes' in a circle normal to the
+ // direction from the vertex to the closest point on the surface
+ // and look for a secondary intersection.
+
+ const vector d = surfHit.hitPoint() - pt;
+
+ const tensor Rp = rotationTensor(vector(0,0,1), np);
+
+ const label s = cvMeshControls().alignmentSearchSpokes();
+
+ const scalar spanMag = geometryToConformTo_.globalBounds().mag();
+
+ scalar totalDist = 0;
+
+ for (label i = 0; i < s; i++)
+ {
+ vector spoke
+ (
+ Foam::cos(i*constant::mathematical::twoPi/s),
+ Foam::sin(i*constant::mathematical::twoPi/s),
+ 0
+ );
+
+ spoke *= spanMag;
+
+ spoke = Rp & spoke;
+
+ pointIndexHit spokeHit;
+
+ label spokeSurface = -1;
+
+ // internal spoke
+
+ geometryToConformTo_.findSurfaceNearestIntersection
+ (
+ pt,
+ pt + spoke,
+ spokeHit,
+ spokeSurface
+ );
+
+ if (spokeHit.hit())
+ {
+ const Foam::point& hitPt = spokeHit.hitPoint();
+
+ scalar spokeHitDistance = mag(hitPt - pt);
+
+ cellSizeHits.append
+ (
+ cellSizeControl().cellSize(hitPt)
+ );
+
+ totalDist += spokeHitDistance;
+ }
+
+ //external spoke
+
+ Foam::point mirrorPt = pt + 2*d;
+
+ geometryToConformTo_.findSurfaceNearestIntersection
+ (
+ mirrorPt,
+ mirrorPt + spoke,
+ spokeHit,
+ spokeSurface
+ );
+
+ if (spokeHit.hit())
+ {
+ const Foam::point& hitPt = spokeHit.hitPoint();
+
+ scalar spokeHitDistance = mag(hitPt - mirrorPt);
+
+ cellSizeHits.append
+ (
+ cellSizeControl().cellSize(hitPt)
+ );
+
+ totalDist += spokeHitDistance;
+ }
+ }
+
+ scalar cellSize = 0;
+
+ forAll(cellSizeHits, hitI)
+ {
+ cellSize += cellSizeHits[hitI];
+ }
+
+ return cellSize/cellSizeHits.size();
+ //return cellSizeControl().cellSize(pt);
+}
+
+
Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
(
const Foam::point& pt
@@ -828,10 +964,7 @@ void Foam::conformalVoronoiMesh::storeSizesAndAlignments
{
sizeAndAlignmentLocations_[i] = topoint(*pit);
- storedSizes_[i] = cellSizeControl().cellSize
- (
- sizeAndAlignmentLocations_[i]
- );
+ storedSizes_[i] = requiredSize(sizeAndAlignmentLocations_[i]);
storedAlignments_[i] = requiredAlignment(sizeAndAlignmentLocations_[i]);
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 81ce496c7173d96a4ef4f021d9d4a80d440603db..f9648dc1e7749c6cd78f444c5a872d8aefbc6fe2 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -268,6 +268,9 @@ private:
//- Return the local maximum surface protrusion distance
inline scalar maxSurfaceProtrusion(const Foam::point& pt) const;
+ //- Return the required cell size at the given location
+ scalar requiredSize(const Foam::point& pt) const;
+
//- Return the required alignment directions at the given location
tensor requiredAlignment(const Foam::point& pt) const;
@@ -533,6 +536,11 @@ private:
const Delaunay::Finite_edges_iterator& eit
) const;
+ boolList dualFaceBoundaryPoints
+ (
+ const Delaunay::Finite_edges_iterator& eit
+ ) const;
+
//- Finds the maximum filterCount of the dual vertices
// (Delaunay cells) that form the dual face produced by the
// supplied edge
@@ -839,6 +847,13 @@ private:
const Delaunay::Finite_facets_iterator& fit
) const;
+ //- Merge adjacent edges that are not attached to other faces
+ label mergeNearlyParallelEdges
+ (
+ const pointField& pts,
+ const scalar maxCosAngle
+ );
+
//- Merge vertices that are very close together
void mergeCloseDualVertices
(
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index 5586446cb8a17934c284d07a017d50ff32565908..e0fe697e88e437c90713bb15ba378c9b5d1ca30d 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -212,6 +212,20 @@ void Foam::conformalVoronoiMesh::calcDualMesh
Info<< nl << "Collapsing unnecessary faces" << endl;
collapseFaces(points, boundaryPts, deferredCollapseFaces);
+
+ const scalar maxCosAngle
+ = cos(degToRad(cvMeshControls().edgeMergeAngle()));
+
+ Info<< nl << "Merging adjacent edges which have an angle "
+ << "of greater than "
+ << cvMeshControls().edgeMergeAngle() << ": " << endl;
+
+ label nRemovedEdges =
+ mergeNearlyParallelEdges(points, maxCosAngle);
+
+ reduce(nRemovedEdges, sumOp<label>());
+
+ Info<< " Merged " << nRemovedEdges << " edges" << endl;
}
labelHashSet wrongFaces = checkPolyMeshQuality(points);
@@ -652,6 +666,113 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
}
+Foam::label Foam::conformalVoronoiMesh::mergeNearlyParallelEdges
+(
+ const pointField& pts,
+ const scalar maxCosAngle
+)
+{
+ List<HashSet<label> > pointFaceCount(number_of_cells());
+ labelList pointNeighbour(number_of_cells(), -1);
+
+ Map<label> dualPtIndexMap;
+
+ 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 (isBoundaryDualFace(eit))
+ {
+ const face f = buildDualFace(eit);
+
+ forAll(f, pI)
+ {
+ const label pIndex = f[pI];
+
+ const label prevPointI = f.prevLabel(pI);
+ const label nextPointI = f.nextLabel(pI);
+
+ pointFaceCount[pIndex].insert(prevPointI);
+ pointFaceCount[pIndex].insert(nextPointI);
+ pointNeighbour[pIndex] = nextPointI;
+ }
+ }
+ else if
+ (
+ vA->internalOrBoundaryPoint()
+ || vB->internalOrBoundaryPoint()
+ )
+ {
+ const face f = buildDualFace(eit);
+ const boolList faceBoundaryPoints = dualFaceBoundaryPoints(eit);
+
+ forAll(f, pI)
+ {
+ const label pIndex = f[pI];
+
+ const label prevPointI = f.prevLabel(pI);
+ const label nextPointI = f.nextLabel(pI);
+
+ pointFaceCount[pIndex].insert(prevPointI);
+ pointFaceCount[pIndex].insert(nextPointI);
+
+ if (faceBoundaryPoints[pI] == false)
+ {
+ pointNeighbour[pIndex] = nextPointI;
+ }
+ else
+ {
+ if (faceBoundaryPoints[prevPointI] == true)
+ {
+ pointNeighbour[pIndex] = prevPointI;
+ }
+ else if (faceBoundaryPoints[nextPointI] == true)
+ {
+ pointNeighbour[pIndex] = nextPointI;
+ }
+ else
+ {
+ pointNeighbour[pIndex] = pIndex;
+ }
+ }
+ }
+ }
+ }
+
+ forAll(pointFaceCount, pI)
+ {
+ if (pointFaceCount[pI].size() == 2)
+ {
+ List<vector> edges(2, vector(0, 0, 0));
+
+ label count = 0;
+ forAllConstIter(HashSet<label>, pointFaceCount[pI], iter)
+ {
+ edges[count] = pts[pI] - pts[iter.key()];
+ edges[count] /= mag(edges[count]) + VSMALL;
+ count++;
+ }
+
+ if (mag(edges[0] & edges[1]) > maxCosAngle)
+ {
+ dualPtIndexMap.insert(pI, pointNeighbour[pI]);
+ }
+ }
+ }
+
+ reindexDualVertices(dualPtIndexMap);
+
+ return dualPtIndexMap.size();
+}
+
+
void Foam::conformalVoronoiMesh::smoothSurface
(
pointField& pts,
@@ -696,7 +817,6 @@ void Foam::conformalVoronoiMesh::smoothSurface
} while (nCollapsedFaces > 0);
// Force all points of boundary faces to be on the surface
-
for
(
Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@@ -917,6 +1037,7 @@ void Foam::conformalVoronoiMesh::collapseFaces
mergeCloseDualVertices(pts, boundaryPts);
+
if (nCollapsedFaces > 0)
{
Info<< " Collapsed " << nCollapsedFaces << " faces" << endl;
@@ -932,6 +1053,7 @@ void Foam::conformalVoronoiMesh::collapseFaces
}
} while (nCollapsedFaces > 0);
+
}
@@ -1064,6 +1186,25 @@ Foam::conformalVoronoiMesh::collapseFace
bool allowEarlyCollapseToPoint = true;
+
+// // Quick exit
+// label smallEdges = 0;
+// const edgeList& fEdges = f.edges();
+// forAll(fEdges, eI)
+// {
+// const edge& e = fEdges[eI];
+//
+// if (e.mag(pts) < 0.2*targetFaceSize)
+// {
+// smallEdges++;
+// }
+// }
+// if (smallEdges == 0)
+// {
+// return fcmNone;
+// }
+
+
// if (maxFC > cvMeshControls().filterCountSkipThreshold() - 3)
// {
// limitToQuadsOrTris = true;
@@ -1071,6 +1212,7 @@ Foam::conformalVoronoiMesh::collapseFace
// allowEarlyCollapseToPoint = false;
// }
+
collapseSizeLimitCoeff *= pow
(
cvMeshControls().filterErrorReductionCoeff(),
@@ -1158,6 +1300,91 @@ Foam::conformalVoronoiMesh::collapseFace
}
}
+
+ scalar maxDist = 0;
+ scalar minDist = GREAT;
+//
+// if (f.size() <= 3)
+// {
+// const edgeList& fEdges = f.edges();
+//
+// forAll(fEdges, eI)
+// {
+// const edge& e = fEdges[eI];
+// const scalar d = e.mag(pts);
+//
+// if (d > maxDist)
+// {
+// maxDist = d;
+// collapseAxis = e.vec(pts);
+// }
+// else if (d < minDist && d != 0)
+// {
+// minDist = d;
+// }
+// }
+// }
+// else
+// {
+// forAll(f, pI)
+// {
+// for (label i = pI + 1; i < f.size(); ++i)
+// {
+// if
+// (
+// f[i] != f.nextLabel(pI)
+// && f[i] != f.prevLabel(pI)
+// )
+// {
+// scalar d = mag(pts[f[pI]] - pts[f[i]]);
+//
+// if (d > maxDist)
+// {
+// maxDist = d;
+// collapseAxis = pts[f[pI]] - pts[f[i]];
+// }
+// else if (d < minDist && d != 0)
+// {
+// minDist = d;
+// }
+// }
+// }
+// }
+// }
+//
+// const edgeList& fEdges = f.edges();
+//
+// scalar perimeter = 0;
+//
+// forAll(fEdges, eI)
+// {
+// const edge& e = fEdges[eI];
+// const scalar d = e.mag(pts);
+//
+// perimeter += d;
+//
+//// collapseAxis += e.vec(pts);
+//
+// if (d > maxDist)
+// {
+// collapseAxis = e.vec(pts);
+// maxDist = d;
+// }
+// else if (d < minDist && d != 0)
+// {
+// minDist = d;
+// }
+// }
+//
+// collapseAxis /= mag(collapseAxis);
+//
+//// Info<< f.size() << " " << minDist << " " << maxDist << " | "
+//// << collapseAxis << endl;
+//
+// aspectRatio = maxDist/minDist;
+//
+// aspectRatio = min(aspectRatio, sqr(perimeter)/(16.0*fA));
+
if (magSqr(collapseAxis) < VSMALL)
{
WarningIn
@@ -1170,9 +1397,9 @@ Foam::conformalVoronoiMesh::collapseFace
// Output face and collapse axis for visualisation
Pout<< "# Aspect ratio = " << aspectRatio << nl
- << "# inertia = " << J << nl
- << "# determinant = " << detJ << nl
- << "# eigenvalues = " << eigenValues(J) << nl
+// << "# inertia = " << J << nl
+// << "# determinant = " << detJ << nl
+// << "# eigenvalues = " << eigenValues(J) << nl
<< "# collapseAxis = " << collapseAxis << nl
<< "# facePts = " << facePts << nl
<< endl;
@@ -1280,8 +1507,6 @@ Foam::conformalVoronoiMesh::collapseFace
{
scalar guardFraction = cvMeshControls().edgeCollapseGuardFraction();
- cvMeshControls().maxCollapseFaceToPointSideLengthCoeff();
-
if
(
allowEarlyCollapseToPoint
@@ -1337,6 +1562,8 @@ Foam::conformalVoronoiMesh::collapseFace
Foam::point collapseToPt =
collapseAxis*(sum(dNeg)/dNeg.size() - dShift) + fC;
+// DynamicList<label> faceBoundaryPts(f.size());
+
forAll(facePtsNeg, fPtI)
{
if (boundaryPts[facePtsNeg[fPtI]] == true)
@@ -1351,9 +1578,36 @@ Foam::conformalVoronoiMesh::collapseFace
collapseToPt = pts[collapseToPtI];
break;
+
+// faceBoundaryPts.append(facePtsNeg[fPtI]);
}
}
+// if (!faceBoundaryPts.empty())
+// {
+// if (faceBoundaryPts.size() == 2)
+// {
+// collapseToPtI = faceBoundaryPts[0];
+//
+// collapseToPt =
+// 0.5
+// *(
+// pts[faceBoundaryPts[0]]
+// + pts[faceBoundaryPts[1]]
+// );
+// }
+// else if (faceBoundaryPts.size() < f.size())
+// {
+// face bFace(faceBoundaryPts);
+//
+// collapseToPtI = faceBoundaryPts.first();
+//
+// collapseToPt = bFace.centre(pts);
+// }
+// }
+//
+// faceBoundaryPts.clear();
+
// ...otherwise arbitrarily choosing the most distant
// point as the index to collapse to.
@@ -1384,9 +1638,34 @@ Foam::conformalVoronoiMesh::collapseFace
collapseToPt = pts[collapseToPtI];
break;
+
+// faceBoundaryPts.append(facePtsNeg[fPtI]);
}
}
+// if (!faceBoundaryPts.empty())
+// {
+// if (faceBoundaryPts.size() == 2)
+// {
+// collapseToPtI = faceBoundaryPts[0];
+//
+// collapseToPt =
+// 0.5
+// *(
+// pts[faceBoundaryPts[0]]
+// + pts[faceBoundaryPts[1]]
+// );
+// }
+// else if (faceBoundaryPts.size() < f.size())
+// {
+// face bFace(faceBoundaryPts);
+//
+// collapseToPtI = faceBoundaryPts.first();
+//
+// collapseToPt = bFace.centre(pts);
+// }
+// }
+
// ...otherwise arbitrarily choosing the most distant
// point as the index to collapse to.
@@ -1406,6 +1685,8 @@ Foam::conformalVoronoiMesh::collapseFace
Foam::point collapseToPt = fC;
+ DynamicList<label> faceBoundaryPts(f.size());
+
forAll(facePts, fPtI)
{
if (boundaryPts[facePts[fPtI]] == true)
@@ -1415,11 +1696,32 @@ Foam::conformalVoronoiMesh::collapseFace
// use the first boundary point encountered if
// there are multiple boundary points.
- collapseToPtI = facePts[fPtI];
+// collapseToPtI = facePts[fPtI];
+//
+// collapseToPt = pts[collapseToPtI];
+//
+// break;
- collapseToPt = pts[collapseToPtI];
+ faceBoundaryPts.append(facePts[fPtI]);
+ }
+ }
- break;
+ if (!faceBoundaryPts.empty())
+ {
+ if (faceBoundaryPts.size() == 2)
+ {
+ collapseToPtI = faceBoundaryPts[0];
+
+ collapseToPt =
+ 0.5*(pts[faceBoundaryPts[0]] + pts[faceBoundaryPts[1]]);
+ }
+ else if (faceBoundaryPts.size() < f.size())
+ {
+ face bFace(faceBoundaryPts);
+
+ collapseToPtI = faceBoundaryPts.first();
+
+ collapseToPt = bFace.centre(pts);
}
}
@@ -1791,9 +2093,9 @@ void Foam::conformalVoronoiMesh::checkCellSizing()
}
Info<< " Automatically re-sizing " << cellsToResize.size()
- << " cells that are attached to the bad faces: DISABLED" << endl;
+ << " cells that are attached to the bad faces: " << endl;
- //cellSizeControl_.setCellSizes(cellsToResize);
+ cellSizeControl_.setCellSizes(cellsToResize);
}
timeCheck("End of Cell Sizing");
@@ -1858,7 +2160,10 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findOffsetPatchFaces
}
}
- offsetBoundaryCells.write();
+ if (cvMeshControls().objOutput())
+ {
+ offsetBoundaryCells.write();
+ }
return offsetBoundaryCells;
}
@@ -2092,13 +2397,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
pts[dualVertI] = cit->dual();
- if
- (
- !cit->vertex(0)->internalOrBoundaryPoint()
- || !cit->vertex(1)->internalOrBoundaryPoint()
- || !cit->vertex(2)->internalOrBoundaryPoint()
- || !cit->vertex(3)->internalOrBoundaryPoint()
- )
+ if (cit->boundaryDualVertex())
{
// This is a boundary dual vertex
boundaryPts[dualVertI] = true;
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index bee04b2fbddc8a165fc72f6f46cdee5b9db1185b..989fe2cf95c42f70514e866d3b50bcb4f8d03898 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -719,12 +719,15 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
const Foam::point& p = infoList[hitI].hitPoint();
- const scalar separationDistance
- = mag(p - info.hitPoint());
+ const scalar separationDistance =
+ mag(p - info.hitPoint());
- const scalar minSepDist
- = sqr(cvMeshControls().removalDistCoeff()
- *targetCellSize(p));
+ const scalar minSepDist =
+ sqr
+ (
+ cvMeshControls().removalDistCoeff()
+ *targetCellSize(p)
+ );
// Reject the point if it is too close to another
// surface point.
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index 1e24b7d8ef29cf030425647cfc70fc03f50bf26a..43d8c077fc785724ad83b37ac8e48a7a76f77ba9 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -366,6 +366,49 @@ inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
}
+inline Foam::List<bool> Foam::conformalVoronoiMesh::dualFaceBoundaryPoints
+(
+ const Delaunay::Finite_edges_iterator& eit
+) const
+{
+ Cell_circulator ccStart = incident_cells(*eit);
+ Cell_circulator cc1 = ccStart;
+ Cell_circulator cc2 = cc1;
+
+ // Advance the second circulator so that it always stays on the next
+ // cell around the edge;
+ cc2++;
+
+ DynamicList<bool> tmpFaceBoundaryPoints;
+
+ do
+ {
+ label cc1I = cc1->cellIndex();
+
+ label cc2I = cc2->cellIndex();
+
+ if (cc1I != cc2I)
+ {
+ if (cc1->boundaryDualVertex())
+ {
+ tmpFaceBoundaryPoints.append(true);
+ }
+ else
+ {
+ tmpFaceBoundaryPoints.append(false);
+ }
+ }
+
+ cc1++;
+
+ cc2++;
+
+ } while (cc1 != ccStart);
+
+ return tmpFaceBoundaryPoints;
+}
+
+
inline Foam::List<Foam::label> Foam::conformalVoronoiMesh::processorsAttached
(
const Delaunay::Finite_facets_iterator& fit
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
index c99f9324ff1c02da60ee84c7275f76ba7415941e..996aa71abfcf13e62ba123e3519090121d5ff095 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
@@ -291,6 +291,11 @@ Foam::cvControls::cvControls
filteringDict.lookup("mergeClosenessCoeff")
);
+ edgeMergeAngle_ = readScalar
+ (
+ filteringDict.lookup("edgeMergeAngle")
+ );
+
continueFilteringOnBadInitialPolyMesh_ = Switch
(
filteringDict.lookupOrDefault<Switch>
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
index dae83ed06dcec1bec4b9b5be5b5c0053953e19c3..58dd06fc741ee6b1c88dc14a5d1cd6756de93692 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
@@ -234,6 +234,11 @@ class cvControls
// being merged, fraction of the local target cell size
scalar mergeClosenessCoeff_;
+ //- If the angle between two dual edges that are connected by a single
+ // point is less than this angle, then the edges will be merged into a
+ // single edge.
+ scalar edgeMergeAngle_;
+
//- If the mesh quality criteria cannot be satisfied, continue
// with filtering anyway?
Switch continueFilteringOnBadInitialPolyMesh_;
@@ -405,6 +410,9 @@ public:
//- Return the mergeClosenessCoeff
inline scalar mergeClosenessCoeff() const;
+ //- Return the edgeMergeAngle
+ inline scalar edgeMergeAngle() const;
+
//- Return the continueFilteringOnBadInitialPolyMesh Switch
inline Switch continueFilteringOnBadInitialPolyMesh() const;
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
index ce45e5f52242a1f8b38c6d8285802ca72e18e859..50df7e8f5841b36d61322b0931f6840a4de52cf1 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
@@ -166,6 +166,12 @@ inline Foam::scalar Foam::cvControls::mergeClosenessCoeff() const
}
+inline Foam::scalar Foam::cvControls::edgeMergeAngle() const
+{
+ return edgeMergeAngle_;
+}
+
+
inline Foam::Switch
Foam::cvControls::continueFilteringOnBadInitialPolyMesh() const
{
diff --git a/applications/utilities/surface/surfaceFeatureExtract/Allwmake b/applications/utilities/surface/surfaceFeatureExtract/Allwmake
index e54f15ae2829ee050db4964162262739637b8226..369428e1713ee11a3f49aadc110376bcaf2aaa79 100755
--- a/applications/utilities/surface/surfaceFeatureExtract/Allwmake
+++ b/applications/utilities/surface/surfaceFeatureExtract/Allwmake
@@ -3,30 +3,25 @@ cd ${0%/*} || exit 1 # run from this directory
set -x
-if [ ! -e "Make/files" ] || [ ! -e "Make/options" ]
+if [ -n "$CGAL_ARCH_PATH" ]
then
- mkdir -p Make
-
- if [ -n "$CGAL_ARCH_PATH" ]
- then
- cp -r MakeWithCGAL/* Make
-
- echo
- echo Compiling surfaceFeatureExtract with CGAL support for curvature
- echo
-
- wmake
- else
- cp -r MakeWithoutCGAL/* Make
-
- echo
- echo Compiling surfaceFeatureExtract without CGAL support for curvature
- echo
-
- wmake
- fi
+ echo
+ echo "Compiling surfaceFeatureExtract with CGAL curvature support"
+ echo
+
+ wmake "ENABLE_CURVATURE=-DENABLE_CURVATURE \
+ EXE_FROUNDING_MATH=-frounding-math \
+ USE_F2C=-DCGAL_USE_F2C \
+ CGAL_LIBDIR=-L$CGAL_ARCH_PATH/lib \
+ LAPACK_LIB=-llapack \
+ BLAS_LIB=-lblas \
+ CGAL_LIB=-lCGAL"
else
- echo surfaceFeatureExtract already has a Make folder
+ echo
+ echo "Compiling surfaceFeatureExtract without CGAL curvature support"
+ echo
+
+ wmake
fi
diff --git a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H
index 3001a247d3ee8435b64545e7e63b5e39e95ef8ae..30627c7463ffd0a91f2b3b5cc8bf0f3deaf260ec 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H
+++ b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H
@@ -53,7 +53,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
- Class buildCGALPolyhedron Declaration
+ Class buildCGALPolyhedron Declaration
\*---------------------------------------------------------------------------*/
class buildCGALPolyhedron
@@ -80,18 +80,56 @@ public:
// Constructors
//- Construct with reference to triSurface
- buildCGALPolyhedron(const triSurface& surf);
+ explicit buildCGALPolyhedron(const triSurface& surf)
+ :
+ CGAL::Modifier_base<HalfedgeDS>(),
+ surf_(surf)
+ {}
//- Destructor
- ~buildCGALPolyhedron();
+ ~buildCGALPolyhedron(){}
// Member Operators
//- operator() of this `modifier' called by delegate function of
// Polyhedron
- void operator()(HalfedgeDS& hds);
+ void operator()(HalfedgeDS& hds)
+ {
+ typedef HalfedgeDS::Traits Traits;
+ typedef Traits::Point_3 Point;
+
+ // Postcondition: `hds' is a valid polyhedral surface.
+ CGAL::Polyhedron_incremental_builder_3<HalfedgeDS> B(hds, false);
+
+ B.begin_surface
+ (
+ surf_.points().size(), // n points
+ surf_.size(), // n facets
+ 2*surf_.edges().size() // n halfedges
+ );
+
+ forAll(surf_.points(), pI)
+ {
+ const Foam::point& p = surf_.points()[pI];
+
+ B.add_vertex(Point(p.x(), p.y(), p.z()));
+ }
+
+ forAll(surf_, fI)
+ {
+ B.begin_facet();
+
+ B.add_vertex_to_facet(surf_[fI][0]);
+ B.add_vertex_to_facet(surf_[fI][1]);
+ B.add_vertex_to_facet(surf_[fI][2]);
+
+ B.end_facet();
+ }
+
+ B.end_surface();
+ }
};
diff --git a/applications/utilities/surface/surfaceFeatureExtract/Make/options b/applications/utilities/surface/surfaceFeatureExtract/Make/options
index b733d5fde961f788415c0aa33e9e74bd1caa539d..b2e5dcf5f7abef5f52c4472ffc0030b91944366e 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/Make/options
+++ b/applications/utilities/surface/surfaceFeatureExtract/Make/options
@@ -1,4 +1,11 @@
+include $(GENERAL_RULES)/CGAL
+
EXE_INC = \
+ ${ENABLE_CURVATURE}\
+ ${EXE_FROUNDING_MATH} \
+ ${USE_F2C} \
+ ${CGAL_INC} \
+ -ICGALPolyhedron \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/edgeMesh/lnInclude \
@@ -7,6 +14,11 @@ EXE_INC = \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
+ $(CGAL_LIBS) \
+ ${CGAL_LIBDIR} \
+ ${LAPACK_LIB} \
+ ${BLAS_LIB} \
+ ${CGAL_LIB} \
-lmeshTools \
-ledgeMesh \
-ltriSurface \
diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files b/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files
deleted file mode 100644
index 0c0f6f7966efb238b4a4b605e978696384379fc5..0000000000000000000000000000000000000000
--- a/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files
+++ /dev/null
@@ -1,4 +0,0 @@
-surfaceFeatureExtract.C
-CGALPolyhedron/buildCGALPolyhedron.C
-
-EXE = $(FOAM_APPBIN)/surfaceFeatureExtract
diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options b/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options
deleted file mode 100644
index 60833572b75ba153f30697475b5ba1eb623b209f..0000000000000000000000000000000000000000
--- a/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options
+++ /dev/null
@@ -1,29 +0,0 @@
-EXE_FROUNDING_MATH = -frounding-math
-EXE_NDEBUG = -DNDEBUG
-USE_F2C = -DCGAL_USE_F2C
-include $(GENERAL_RULES)/CGAL
-
-EXE_INC = \
- -DENABLE_CURVATURE \
- ${EXE_FROUNDING_MATH} \
- ${EXE_NDEBUG} \
- ${USE_F2C} \
- ${CGAL_INC} \
- -ICGALPolyhedron \
- -I$(LIB_SRC)/finiteVolume/lnInclude \
- -I$(LIB_SRC)/meshTools/lnInclude \
- -I$(LIB_SRC)/edgeMesh/lnInclude \
- -I$(LIB_SRC)/triSurface/lnInclude \
- -I$(LIB_SRC)/surfMesh/lnInclude \
- -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
- $(CGAL_LIBS) \
- -L$(CGAL_ARCH_PATH)/lib \
- -llapack \
- -lblas \
- -lCGAL \
- -lmeshTools \
- -ledgeMesh \
- -ltriSurface \
- -lsampling
diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files b/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files
deleted file mode 100644
index a57125cd6a7c24a2892c292ee358cce52a6529ce..0000000000000000000000000000000000000000
--- a/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files
+++ /dev/null
@@ -1,3 +0,0 @@
-surfaceFeatureExtract.C
-
-EXE = $(FOAM_APPBIN)/surfaceFeatureExtract
diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options b/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options
deleted file mode 100644
index b733d5fde961f788415c0aa33e9e74bd1caa539d..0000000000000000000000000000000000000000
--- a/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options
+++ /dev/null
@@ -1,13 +0,0 @@
-EXE_INC = \
- -I$(LIB_SRC)/finiteVolume/lnInclude \
- -I$(LIB_SRC)/meshTools/lnInclude \
- -I$(LIB_SRC)/edgeMesh/lnInclude \
- -I$(LIB_SRC)/triSurface/lnInclude \
- -I$(LIB_SRC)/surfMesh/lnInclude \
- -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
- -lmeshTools \
- -ledgeMesh \
- -ltriSurface \
- -lsampling
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index da71cafc90b8e63cc3f0d0b0fbb72d2c071c1e23..b02c4c815d77ab492ed5c57d745ae0416a95dcf7 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -74,30 +74,6 @@ scalarField calcCurvature(const triSurface& surf)
// The rest of this function adapted from
// CGAL-3.7/examples/Jet_fitting_3/Mesh_estimation.cpp
- // Licensed under CGAL-3.7/LICENSE.FREE_USE
-
- // Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
- // Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
- // Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
- // Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
- // Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
- // (Israel). All rights reserved.
-
- // Permission is hereby granted, free of charge, to any person obtaining a
- // copy of this software and associated documentation files (the
- // "Software"), to deal in the Software without restriction, including
- // without limitation the rights to use, copy, modify, merge, publish,
- // distribute, sublicense, and/or sell copies of the Software, and to permit
- // persons to whom the Software is furnished to do so, subject to the
- // following conditions:
-
- // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
- // IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
- // CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT
- // OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
- // THE USE OR OTHER DEALINGS IN THE SOFTWARE.
//Vertex property map, with std::map
typedef std::map<Vertex*, int> Vertex2int_map_type;
@@ -181,12 +157,12 @@ scalarField calcCurvature(const triSurface& surf)
// << monge_fit.condition_number() << nl << std::endl;
// Use the maximum curvature to give smaller cell sizes later.
- k[vertI++]
- = max
- (
- mag(monge_form.principal_curvatures(0)),
- mag(monge_form.principal_curvatures(1))
- );
+ k[vertI++] =
+ max
+ (
+ mag(monge_form.principal_curvatures(0)),
+ mag(monge_form.principal_curvatures(1))
+ );
}
return k;
@@ -198,8 +174,10 @@ bool edgesConnected(const edge& e1, const edge& e2)
{
if
(
- e1.start() == e2.start() || e1.start() == e2.end()
- || e1.end() == e2.start() || e1.end() == e2.end()
+ e1.start() == e2.start()
+ || e1.start() == e2.end()
+ || e1.end() == e2.start()
+ || e1.end() == e2.end()
)
{
return true;
@@ -289,8 +267,8 @@ scalar calcProximityOfFeatureEdges
// Don't refine if the edges are connected to each other
if (!edgesConnected(e1, e2))
{
- scalar curDist
- = mag(pHit1.hitPoint() - pHit2.hitPoint());
+ scalar curDist =
+ mag(pHit1.hitPoint() - pHit2.hitPoint());
minDist = min(curDist, minDist);
}
@@ -503,353 +481,344 @@ int main(int argc, char *argv[])
"extract and write surface features to file"
);
argList::noParallel();
- argList::validArgs.append("surface");
- argList::validArgs.append("output set");
argList::addOption
(
- "includedAngle",
- "degrees",
- "construct feature set from included angle [0..180]"
- );
- argList::addOption
- (
- "set",
- "name",
- "use existing feature set from file"
- );
- argList::addOption
- (
- "minLen",
- "scalar",
- "remove features shorter than the specified cumulative length"
- );
- argList::addOption
- (
- "minElem",
- "int",
- "remove features with fewer than the specified number of edges"
- );
- argList::addOption
- (
- "subsetBox",
- "((x0 y0 z0)(x1 y1 z1))",
- "remove edges outside specified bounding box"
- );
- argList::addOption
- (
- "deleteBox",
- "((x0 y0 z0)(x1 y1 z1))",
- "remove edges within specified bounding box"
+ "dict",
+ "word",
+ "name of dictionary to provide feature extraction information"
);
- argList::addBoolOption
- (
- "writeObj",
- "write extendedFeatureEdgeMesh obj files"
- );
- argList::addBoolOption
- (
- "writeVTK",
- "write extendedFeatureEdgeMesh vtk files"
- );
- argList::addBoolOption
- (
- "closeness",
- "span to look for surface closeness"
- );
- argList::addOption
- (
- "featureProximity",
- "scalar",
- "distance to look for close features"
- );
- argList::addBoolOption
- (
- "manifoldEdgesOnly",
- "remove any non-manifold (open or more than two connected faces) edges"
- );
- argList::addOption
+
+# include "setRootCase.H"
+# include "createTime.H"
+
+ word dictName
(
- "plane",
- "(nx ny nz)(z0 y0 z0)",
- "use a plane to create feature edges for 2D meshing"
+ args.optionLookupOrDefault<word>("dict", "surfaceFeatureExtractDict")
);
-# ifdef ENABLE_CURVATURE
- argList::addBoolOption
+ Info<< "Reading " << dictName << nl << endl;
+
+ IOdictionary dict
(
- "curvature",
- "calculate curvature and closeness fields"
+ IOobject
+ (
+ dictName,
+ runTime.system(),
+ runTime,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE,
+ false
+ )
);
-# endif
+ forAllConstIter(dictionary, dict, iter)
+ {
+ dictionary surfaceDict = dict.subDict(iter().keyword());
-# include "setRootCase.H"
-# include "createTime.H"
+ const fileName surfFileName = iter().keyword();
+ const fileName sFeatFileName = surfFileName.lessExt().name();
- bool writeVTK = args.optionFound("writeVTK");
+ Info<< "Surface : " << surfFileName << nl << endl;
- bool writeObj = args.optionFound("writeObj");
+ const Switch writeVTK =
+ surfaceDict.lookupOrAddDefault<Switch>("writeVTK", "off");
+ const Switch writeObj =
+ surfaceDict.lookupOrAddDefault<Switch>("writeObj", "off");
+ const Switch writeFeatureEdgeMesh =
+ surfaceDict.lookupOrAddDefault<Switch>
+ (
+ "writeFeatureEdgeMesh",
+ "off"
+ );
- bool curvature = args.optionFound("curvature");
+ const Switch curvature =
+ surfaceDict.lookupOrAddDefault<Switch>("curvature", "off");
+ const Switch featureProximity =
+ surfaceDict.lookupOrAddDefault<Switch>("featureProximity", "off");
+ const Switch closeness =
+ surfaceDict.lookupOrAddDefault<Switch>("closeness", "off");
- if (curvature && env("FOAM_SIGFPE"))
- {
- WarningIn(args.executable())
- << "Detected floating point exception trapping (FOAM_SIGFPE)."
- << " This might give" << nl
- << " problems when calculating curvature on straight angles"
- << " (infinite curvature)" << nl
- << " Switch it off in case of problems." << endl;
- }
+ const word extractionMethod = surfaceDict.lookup("extractionMethod");
- Info<< "Feature line extraction is only valid on closed manifold surfaces."
- << endl;
+#ifndef ENABLE_CURVATURE
+ if (curvature)
+ {
+ WarningIn(args.executable())
+ << "Curvature calculation has been requested but "
+ << args.executable() << " has not " << nl
+ << " been compiled with CGAL. "
+ << "Skipping the curvature calculation." << endl;
+ }
+#else
+ if (curvature && env("FOAM_SIGFPE"))
+ {
+ WarningIn(args.executable())
+ << "Detected floating point exception trapping (FOAM_SIGFPE)."
+ << " This might give" << nl
+ << " problems when calculating curvature on straight angles"
+ << " (infinite curvature)" << nl
+ << " Switch it off in case of problems." << endl;
+ }
+#endif
- const fileName surfFileName = args[1];
- const fileName outFileName = args[2];
- Info<< "Surface : " << surfFileName << nl
- << "Output feature set : " << outFileName << nl
- << endl;
+ Info<< nl << "Feature line extraction is only valid on closed manifold "
+ << "surfaces." << endl;
- fileName sFeatFileName = surfFileName.lessExt().name();
+ // Read
+ // ~~~~
+ triSurface surf("constant/triSurface/" + surfFileName);
- // Read
- // ~~~~
+ Info<< "Statistics:" << endl;
+ surf.writeStats(Info);
+ Info<< endl;
- triSurface surf(surfFileName);
+ faceList faces(surf.size());
- Info<< "Statistics:" << endl;
- surf.writeStats(Info);
- Info<< endl;
+ forAll(surf, fI)
+ {
+ faces[fI] = surf[fI].triFaceFace();
+ }
- faceList faces(surf.size());
- forAll(surf, fI)
- {
- faces[fI] = surf[fI].triFaceFace();
- }
+ // Either construct features from surface & featureAngle or read set.
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- // Either construct features from surface&featureangle or read set.
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ surfaceFeatures set(surf);
- surfaceFeatures set(surf);
+ if (extractionMethod == "extractFromFile")
+ {
+ const fileName featureEdgeFile =
+ surfaceDict.subDict("extractFromFile").lookup
+ (
+ "featureEdgeFile"
+ );
- if (args.optionFound("set"))
- {
- const fileName setName = args["set"];
+ edgeMesh eMesh(featureEdgeFile);
- Info<< "Reading existing feature set from file " << setName << endl;
+ // Sometimes duplicate edges are present. Remove them.
+ eMesh.mergeEdges();
- set = surfaceFeatures(surf, setName);
- }
- else if (args.optionFound("includedAngle"))
- {
- const scalar includedAngle = args.optionRead<scalar>("includedAngle");
+ Info<< nl << "Reading existing feature edges from file "
+ << featureEdgeFile << endl;
+
+ set = surfaceFeatures(surf, eMesh.points(), eMesh.edges());
+ }
+ else if (extractionMethod == "extractFromSurface")
+ {
+ const scalar includedAngle =
+ readScalar
+ (
+ surfaceDict.subDict("extractFromSurface").lookup
+ (
+ "includedAngle"
+ )
+ );
+
+ Info<< nl << "Constructing feature set from included angle "
+ << includedAngle << endl;
- Info<< "Constructing feature set from included angle " << includedAngle
+ set = surfaceFeatures(surf, includedAngle);
+ }
+ else
+ {
+ FatalErrorIn(args.executable())
+ << "No initial feature set. Provide either one"
+ << " of extractFromFile (to read existing set)" << nl
+ << " or extractFromSurface (to construct new set from angle)"
+ << exit(FatalError);
+ }
+
+ Info<< nl
+ << "Initial feature set:" << nl
+ << " feature points : " << set.featurePoints().size() << nl
+ << " feature edges : " << set.featureEdges().size() << nl
+ << " of which" << nl
+ << " region edges : " << set.nRegionEdges() << nl
+ << " external edges : " << set.nExternalEdges() << nl
+ << " internal edges : " << set.nInternalEdges() << nl
<< endl;
- set = surfaceFeatures(surf, includedAngle);
- // Info<< nl << "Writing initial features" << endl;
- // set.write("initial.fSet");
- // set.writeObj("initial");
- }
- else
- {
- FatalErrorIn(args.executable())
- << "No initial feature set. Provide either one"
- << " of -set (to read existing set)" << nl
- << " or -includedAngle (to new set construct from angle)"
- << exit(FatalError);
- }
+ // Trim set
+ // ~~~~~~~~
- Info<< nl
- << "Initial feature set:" << nl
- << " feature points : " << set.featurePoints().size() << nl
- << " feature edges : " << set.featureEdges().size() << nl
- << " of which" << nl
- << " region edges : " << set.nRegionEdges() << nl
- << " external edges : " << set.nExternalEdges() << nl
- << " internal edges : " << set.nInternalEdges() << nl
- << endl;
+ if (surfaceDict.isDict("trimFeatures"))
+ {
+ dictionary trimDict = surfaceDict.subDict("trimFeatures");
+ scalar minLen =
+ trimDict.lookupOrAddDefault<scalar>("minLen", -GREAT);
- // Trim set
- // ~~~~~~~~
+ label minElem = trimDict.lookupOrAddDefault<label>("minElem", 0);
- scalar minLen = -GREAT;
- if (args.optionReadIfPresent("minLen", minLen))
- {
- Info<< "Removing features of length < " << minLen << endl;
- }
+ // Trim away small groups of features
+ if (minElem > 0 || minLen > 0)
+ {
+ Info<< "Removing features of length < "
+ << minLen << endl;
+ Info<< "Removing features with number of edges < "
+ << minElem << endl;
- label minElem = 0;
- if (args.optionReadIfPresent("minElem", minElem))
- {
- Info<< "Removing features with number of edges < " << minElem << endl;
- }
+ set.trimFeatures(minLen, minElem);
+ }
+ }
- // Trim away small groups of features
- if (minElem > 0 || minLen > 0)
- {
- set.trimFeatures(minLen, minElem);
- Info<< endl << "Removed small features" << endl;
- }
+ // Subset
+ // ~~~~~~
- // Subset
- // ~~~~~~
+ // Convert to marked edges, points
+ List<surfaceFeatures::edgeStatus> edgeStat(set.toStatus());
- // Convert to marked edges, points
- List<surfaceFeatures::edgeStatus> edgeStat(set.toStatus());
+ if (surfaceDict.isDict("subsetFeatures"))
+ {
+ dictionary subsetDict = surfaceDict.subDict("subsetFeatures");
- if (args.optionFound("subsetBox"))
- {
- treeBoundBox bb
- (
- args.optionLookup("subsetBox")()
- );
+ if (subsetDict.found("insideBox"))
+ {
+ treeBoundBox bb(subsetDict.lookup("insideBox")());
- Info<< "Removing all edges outside bb " << bb << endl;
- dumpBox(bb, "subsetBox.obj");
+ Info<< "Removing all edges outside bb " << bb << endl;
+ dumpBox(bb, "subsetBox.obj");
- deleteBox
- (
- surf,
- bb,
- false,
- edgeStat
- );
- }
- else if (args.optionFound("deleteBox"))
- {
- treeBoundBox bb
- (
- args.optionLookup("deleteBox")()
- );
+ deleteBox
+ (
+ surf,
+ bb,
+ false,
+ edgeStat
+ );
+ }
+ else if (subsetDict.found("outsideBox"))
+ {
+ treeBoundBox bb(subsetDict.lookup("outsideBox")());
- Info<< "Removing all edges inside bb " << bb << endl;
- dumpBox(bb, "deleteBox.obj");
+ Info<< "Removing all edges inside bb " << bb << endl;
+ dumpBox(bb, "deleteBox.obj");
- deleteBox
- (
- surf,
- bb,
- true,
- edgeStat
- );
- }
+ deleteBox
+ (
+ surf,
+ bb,
+ true,
+ edgeStat
+ );
+ }
- if (args.optionFound("manifoldEdgesOnly"))
- {
- Info<< "Removing all non-manifold edges" << endl;
+ const Switch manifoldEdges =
+ subsetDict.lookupOrAddDefault<Switch>("manifoldEdges", "no");
- forAll(edgeStat, edgeI)
- {
- if (surf.edgeFaces()[edgeI].size() != 2)
+ if (manifoldEdges)
{
- edgeStat[edgeI] = surfaceFeatures::NONE;
+ Info<< "Removing all non-manifold edges" << endl;
+
+ forAll(edgeStat, edgeI)
+ {
+ if (surf.edgeFaces()[edgeI].size() != 2)
+ {
+ edgeStat[edgeI] = surfaceFeatures::NONE;
+ }
+ }
}
- }
- }
- if (args.optionFound("plane"))
- {
- plane cutPlane(args.optionLookup("plane")());
+ if (subsetDict.found("plane"))
+ {
+ plane cutPlane(subsetDict.lookup("plane")());
- deleteEdges
- (
- surf,
- cutPlane,
- edgeStat
- );
+ deleteEdges
+ (
+ surf,
+ cutPlane,
+ edgeStat
+ );
- Info<< "Only edges that intersect the plane with normal "
- << cutPlane.normal() << " and base point " << cutPlane.refPoint()
- << " will be included as feature edges."<< endl;
- }
+ Info<< "Only edges that intersect the plane with normal "
+ << cutPlane.normal()
+ << " and base point " << cutPlane.refPoint()
+ << " will be included as feature edges."<< endl;
+ }
+ }
- surfaceFeatures newSet(surf);
- newSet.setFromStatus(edgeStat);
+ surfaceFeatures newSet(surf);
+ newSet.setFromStatus(edgeStat);
- //Info<< endl << "Writing trimmed features to " << outFileName << endl;
- //newSet.write(outFileName);
+ if (writeObj)
+ {
+ newSet.writeObj("final");
+ }
- // Info<< endl << "Writing edge objs." << endl;
- // newSet.writeObj("final");
+ Info<< nl
+ << "Final feature set after trimming and subsetting:" << nl
+ << " feature points : " << newSet.featurePoints().size() << nl
+ << " feature edges : " << newSet.featureEdges().size() << nl
+ << " of which" << nl
+ << " region edges : " << newSet.nRegionEdges() << nl
+ << " external edges : " << newSet.nExternalEdges() << nl
+ << " internal edges : " << newSet.nInternalEdges() << nl
+ << endl;
- Info<< nl
- << "Final feature set:" << nl
- << " feature points : " << newSet.featurePoints().size() << nl
- << " feature edges : " << newSet.featureEdges().size() << nl
- << " of which" << nl
- << " region edges : " << newSet.nRegionEdges() << nl
- << " external edges : " << newSet.nExternalEdges() << nl
- << " internal edges : " << newSet.nInternalEdges() << nl
- << endl;
+ // Extracting and writing a extendedFeatureEdgeMesh
+ extendedFeatureEdgeMesh feMesh
+ (
+ newSet,
+ runTime,
+ sFeatFileName + ".extendedFeatureEdgeMesh"
+ );
- // Extracting and writing a extendedFeatureEdgeMesh
- extendedFeatureEdgeMesh feMesh
- (
- newSet,
- runTime,
- sFeatFileName + ".extendedFeatureEdgeMesh"
- );
+ Info<< nl << "Writing extendedFeatureEdgeMesh to "
+ << feMesh.objectPath() << endl;
- Info<< nl << "Writing extendedFeatureEdgeMesh to " << feMesh.objectPath()
- << endl;
+ if (writeObj)
+ {
+ feMesh.writeObj(surfFileName.lessExt().name());
+ }
- if (writeObj)
- {
- feMesh.writeObj(surfFileName.lessExt().name());
- }
+ feMesh.write();
+
+ // Write a featureEdgeMesh for backwards compatibility
+ if (writeFeatureEdgeMesh)
+ {
+ featureEdgeMesh bfeMesh
+ (
+ IOobject
+ (
+ surfFileName.lessExt().name() + ".eMesh", // name
+ runTime.constant(), // instance
+ "triSurface",
+ runTime, // registry
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE,
+ false
+ ),
+ feMesh.points(),
+ feMesh.edges()
+ );
- feMesh.write();
+ Info<< nl << "Writing featureEdgeMesh to "
+ << bfeMesh.objectPath() << endl;
- // Write a featureEdgeMesh for backwards compatibility
- {
- featureEdgeMesh bfeMesh
+ bfeMesh.regIOobject::write();
+ }
+
+ triSurfaceMesh searchSurf
(
IOobject
(
- surfFileName.lessExt().name() + ".eMesh", // name
- runTime.constant(), // instance
+ sFeatFileName + ".closeness",
+ runTime.constant(),
"triSurface",
- runTime, // registry
+ runTime,
IOobject::NO_READ,
- IOobject::AUTO_WRITE,
- false
+ IOobject::NO_WRITE
),
- feMesh.points(),
- feMesh.edges()
+ surf
);
- Info<< nl << "Writing featureEdgeMesh to "
- << bfeMesh.objectPath() << endl;
-
- bfeMesh.regIOobject::write();
- }
-
- triSurfaceMesh searchSurf
- (
- IOobject
- (
- sFeatFileName + ".closeness",
- runTime.constant(),
- "triSurface",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- surf
- );
-
// Find close features
// // Dummy trim operation to mark features
@@ -917,72 +886,57 @@ int main(int argc, char *argv[])
// )
// );
- if (args.optionFound("closeness"))
- {
- Info<< nl << "Extracting internal and external closeness of surface."
- << endl;
-
- // Internal and external closeness
-
- // Prepare start and end points for intersection tests
+ if (closeness)
+ {
+ Info<< nl << "Extracting internal and external closeness of "
+ << "surface." << endl;
- const vectorField& normals = searchSurf.faceNormals();
+ // Internal and external closeness
- scalar span = searchSurf.bounds().mag();
+ // Prepare start and end points for intersection tests
- //args.optionReadIfPresent("closeness", span);
+ const vectorField& normals = searchSurf.faceNormals();
- scalar externalAngleTolerance = 10;
- scalar externalToleranceCosAngle = Foam::cos
- (
- degToRad(180 - externalAngleTolerance)
- );
+ scalar span = searchSurf.bounds().mag();
- scalar internalAngleTolerance = 45;
- scalar internalToleranceCosAngle = Foam::cos
- (
- degToRad(180 - internalAngleTolerance)
- );
+ scalar externalAngleTolerance = 10;
+ scalar externalToleranceCosAngle =
+ Foam::cos
+ (
+ degToRad(180 - externalAngleTolerance)
+ );
- Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle << nl
- << "internalToleranceCosAngle: " << internalToleranceCosAngle
- << endl;
+ scalar internalAngleTolerance = 45;
+ scalar internalToleranceCosAngle =
+ Foam::cos
+ (
+ degToRad(180 - internalAngleTolerance)
+ );
- // Info<< "span " << span << endl;
+ Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle
+ << nl
+ << "internalToleranceCosAngle: " << internalToleranceCosAngle
+ << endl;
- pointField start = searchSurf.faceCentres() - span*normals;
- pointField end = searchSurf.faceCentres() + span*normals;
- const pointField& faceCentres = searchSurf.faceCentres();
+ // Info<< "span " << span << endl;
- List<List<pointIndexHit> > allHitInfo;
+ pointField start = searchSurf.faceCentres() - span*normals;
+ pointField end = searchSurf.faceCentres() + span*normals;
+ const pointField& faceCentres = searchSurf.faceCentres();
- // Find all intersections (in order)
- searchSurf.findLineAll(start, end, allHitInfo);
+ List<List<pointIndexHit> > allHitInfo;
- scalarField internalCloseness(start.size(), GREAT);
- scalarField externalCloseness(start.size(), GREAT);
+ // Find all intersections (in order)
+ searchSurf.findLineAll(start, end, allHitInfo);
- forAll(allHitInfo, fI)
- {
- const List<pointIndexHit>& hitInfo = allHitInfo[fI];
+ scalarField internalCloseness(start.size(), GREAT);
+ scalarField externalCloseness(start.size(), GREAT);
- if (hitInfo.size() < 1)
+ forAll(allHitInfo, fI)
{
- drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
+ const List<pointIndexHit>& hitInfo = allHitInfo[fI];
- // FatalErrorIn(args.executable())
- // << "findLineAll did not hit its own face."
- // << exit(FatalError);
- }
- else if (hitInfo.size() == 1)
- {
- if (!hitInfo[0].hit())
- {
- // FatalErrorIn(args.executable())
- // << "findLineAll did not hit any face."
- // << exit(FatalError);
- }
- else if (hitInfo[0].index() != fI)
+ if (hitInfo.size() < 1)
{
drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
@@ -990,283 +944,339 @@ int main(int argc, char *argv[])
// << "findLineAll did not hit its own face."
// << exit(FatalError);
}
- }
- else
- {
- label ownHitI = -1;
-
- forAll(hitInfo, hI)
+ else if (hitInfo.size() == 1)
{
- // Find the hit on the triangle that launched the ray
-
- if (hitInfo[hI].index() == fI)
+ if (!hitInfo[0].hit())
{
- ownHitI = hI;
-
- break;
+ // FatalErrorIn(args.executable())
+ // << "findLineAll did not hit any face."
+ // << exit(FatalError);
}
- }
-
- if (ownHitI < 0)
- {
- drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
-
- // FatalErrorIn(args.executable())
- // << "findLineAll did not hit its own face."
- // << exit(FatalError);
- }
- else if (ownHitI == 0)
- {
- // There are no internal hits, the first hit is the closest
- // external hit
-
- if
- (
- (normals[fI] & normals[hitInfo[ownHitI + 1].index()])
- < externalToleranceCosAngle
- )
+ else if (hitInfo[0].index() != fI)
{
- externalCloseness[fI] = mag
+ drawHitProblem
(
- faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
+ fI,
+ surf,
+ start,
+ faceCentres,
+ end,
+ hitInfo
);
+
+ // FatalErrorIn(args.executable())
+ // << "findLineAll did not hit its own face."
+ // << exit(FatalError);
}
}
- else if (ownHitI == hitInfo.size() - 1)
+ else
{
- // There are no external hits, the last but one hit is the
- // closest internal hit
+ label ownHitI = -1;
- if
- (
- (normals[fI] & normals[hitInfo[ownHitI - 1].index()])
- < internalToleranceCosAngle
- )
+ forAll(hitInfo, hI)
+ {
+ // Find the hit on the triangle that launched the ray
+
+ if (hitInfo[hI].index() == fI)
+ {
+ ownHitI = hI;
+
+ break;
+ }
+ }
+
+ if (ownHitI < 0)
{
- internalCloseness[fI] = mag
+ drawHitProblem
(
- faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
+ fI,
+ surf,
+ start,
+ faceCentres,
+ end,
+ hitInfo
);
+
+ // FatalErrorIn(args.executable())
+ // << "findLineAll did not hit its own face."
+ // << exit(FatalError);
}
- }
- else
- {
- if
- (
- (normals[fI] & normals[hitInfo[ownHitI + 1].index()])
- < externalToleranceCosAngle
- )
+ else if (ownHitI == 0)
{
- externalCloseness[fI] = mag
+ // There are no internal hits, the first hit is the
+ // closest external hit
+
+ if
(
- faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
- );
+ (
+ normals[fI]
+ & normals[hitInfo[ownHitI + 1].index()]
+ )
+ < externalToleranceCosAngle
+ )
+ {
+ externalCloseness[fI] =
+ mag
+ (
+ faceCentres[fI]
+ - hitInfo[ownHitI + 1].hitPoint()
+ );
+ }
}
+ else if (ownHitI == hitInfo.size() - 1)
+ {
+ // There are no external hits, the last but one hit is
+ // the closest internal hit
- if
- (
- (normals[fI] & normals[hitInfo[ownHitI - 1].index()])
- < internalToleranceCosAngle
- )
+ if
+ (
+ (
+ normals[fI]
+ & normals[hitInfo[ownHitI - 1].index()]
+ )
+ < internalToleranceCosAngle
+ )
+ {
+ internalCloseness[fI] =
+ mag
+ (
+ faceCentres[fI]
+ - hitInfo[ownHitI - 1].hitPoint()
+ );
+ }
+ }
+ else
{
- internalCloseness[fI] = mag
+ if
(
- faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
- );
+ (
+ normals[fI]
+ & normals[hitInfo[ownHitI + 1].index()]
+ )
+ < externalToleranceCosAngle
+ )
+ {
+ externalCloseness[fI] =
+ mag
+ (
+ faceCentres[fI]
+ - hitInfo[ownHitI + 1].hitPoint()
+ );
+ }
+
+ if
+ (
+ (
+ normals[fI]
+ & normals[hitInfo[ownHitI - 1].index()]
+ )
+ < internalToleranceCosAngle
+ )
+ {
+ internalCloseness[fI] =
+ mag
+ (
+ faceCentres[fI]
+ - hitInfo[ownHitI - 1].hitPoint()
+ );
+ }
}
}
}
- }
- triSurfaceScalarField internalClosenessField
- (
- IOobject
+ triSurfaceScalarField internalClosenessField
(
- sFeatFileName + ".internalCloseness",
- runTime.constant(),
- "triSurface",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- surf,
- dimLength,
- internalCloseness
- );
+ IOobject
+ (
+ sFeatFileName + ".internalCloseness",
+ runTime.constant(),
+ "triSurface",
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ surf,
+ dimLength,
+ internalCloseness
+ );
- internalClosenessField.write();
+ internalClosenessField.write();
- triSurfaceScalarField externalClosenessField
- (
- IOobject
+ triSurfaceScalarField externalClosenessField
(
- sFeatFileName + ".externalCloseness",
- runTime.constant(),
- "triSurface",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- surf,
- dimLength,
- externalCloseness
- );
+ IOobject
+ (
+ sFeatFileName + ".externalCloseness",
+ runTime.constant(),
+ "triSurface",
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ surf,
+ dimLength,
+ externalCloseness
+ );
- externalClosenessField.write();
+ externalClosenessField.write();
- if (writeVTK)
- {
- vtkSurfaceWriter().write
- (
- runTime.constant()/"triSurface", // outputDir
- sFeatFileName, // surfaceName
- surf.points(),
- faces,
- "internalCloseness", // fieldName
- internalCloseness,
- false, // isNodeValues
- true // verbose
- );
+ if (writeVTK)
+ {
+ vtkSurfaceWriter().write
+ (
+ runTime.constant()/"triSurface", // outputDir
+ sFeatFileName, // surfaceName
+ surf.points(),
+ faces,
+ "internalCloseness", // fieldName
+ internalCloseness,
+ false, // isNodeValues
+ true // verbose
+ );
- vtkSurfaceWriter().write
- (
- runTime.constant()/"triSurface", // outputDir
- sFeatFileName, // surfaceName
- surf.points(),
- faces,
- "externalCloseness", // fieldName
- externalCloseness,
- false, // isNodeValues
- true // verbose
- );
+ vtkSurfaceWriter().write
+ (
+ runTime.constant()/"triSurface", // outputDir
+ sFeatFileName, // surfaceName
+ surf.points(),
+ faces,
+ "externalCloseness", // fieldName
+ externalCloseness,
+ false, // isNodeValues
+ true // verbose
+ );
+ }
}
- }
#ifdef ENABLE_CURVATURE
- if (args.optionFound("curvature"))
- {
- Info<< nl << "Extracting curvature of surface at the points." << endl;
+ if (curvature)
+ {
+ Info<< nl << "Extracting curvature of surface at the points."
+ << endl;
- scalarField k = calcCurvature(surf);
+ scalarField k = calcCurvature(surf);
// Modify the curvature values on feature edges and points to be zero.
- // forAll(newSet.featureEdges(), fEI)
- // {
- // const edge& e = surf.edges()[newSet.featureEdges()[fEI]];
- //
- // k[surf.meshPoints()[e.start()]] = 0.0;
- // k[surf.meshPoints()[e.end()]] = 0.0;
- // }
+ // forAll(newSet.featureEdges(), fEI)
+ // {
+ // const edge& e = surf.edges()[newSet.featureEdges()[fEI]];
+ //
+ // k[surf.meshPoints()[e.start()]] = 0.0;
+ // k[surf.meshPoints()[e.end()]] = 0.0;
+ // }
- triSurfacePointScalarField kField
- (
- IOobject
+ triSurfacePointScalarField kField
(
- sFeatFileName + ".curvature",
- runTime.constant(),
- "triSurface",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- surf,
- dimLength,
- k
- );
+ IOobject
+ (
+ sFeatFileName + ".curvature",
+ runTime.constant(),
+ "triSurface",
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ surf,
+ dimLength,
+ k
+ );
- kField.write();
+ kField.write();
- if (writeVTK)
- {
- vtkSurfaceWriter().write
- (
- runTime.constant()/"triSurface", // outputDir
- sFeatFileName, // surfaceName
- surf.points(),
- faces,
- "curvature", // fieldName
- k,
- true, // isNodeValues
- true // verbose
- );
+ if (writeVTK)
+ {
+ vtkSurfaceWriter().write
+ (
+ runTime.constant()/"triSurface", // outputDir
+ sFeatFileName, // surfaceName
+ surf.points(),
+ faces,
+ "curvature", // fieldName
+ k,
+ true, // isNodeValues
+ true // verbose
+ );
+ }
}
- }
#endif
- if (args.optionFound("featureProximity"))
- {
- Info<< nl << "Extracting proximity of close feature points and edges "
- << "to the surface" << endl;
+ if (featureProximity)
+ {
+ Info<< nl << "Extracting proximity of close feature points and "
+ << "edges to the surface" << endl;
- const scalar searchDistance =
- args.optionRead<scalar>("featureProximity");
+ const scalar searchDistance =
+ readScalar(surfaceDict.lookup("maxFeatureProximity"));
- const scalar radiusSqr = sqr(searchDistance);
+ const scalar radiusSqr = sqr(searchDistance);
- scalarField featureProximity(surf.size(), searchDistance);
+ scalarField featureProximity(surf.size(), searchDistance);
- forAll(surf, fI)
- {
- const triPointRef& tri = surf[fI].tri(surf.points());
- const point& triCentre = tri.circumCentre();
+ forAll(surf, fI)
+ {
+ const triPointRef& tri = surf[fI].tri(surf.points());
+ const point& triCentre = tri.circumCentre();
- List<pointIndexHit> hitList;
+ List<pointIndexHit> hitList;
- feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList);
+ feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList);
- featureProximity[fI] =
- calcProximityOfFeatureEdges
- (
- feMesh,
- hitList,
- featureProximity[fI]
- );
+ featureProximity[fI] =
+ calcProximityOfFeatureEdges
+ (
+ feMesh,
+ hitList,
+ featureProximity[fI]
+ );
- feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList);
+ feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList);
- featureProximity[fI] =
- calcProximityOfFeaturePoints
- (
- hitList,
- featureProximity[fI]
- );
- }
+ featureProximity[fI] =
+ calcProximityOfFeaturePoints
+ (
+ hitList,
+ featureProximity[fI]
+ );
+ }
- triSurfaceScalarField featureProximityField
- (
- IOobject
+ triSurfaceScalarField featureProximityField
(
- sFeatFileName + ".featureProximity",
- runTime.constant(),
- "triSurface",
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- surf,
- dimLength,
- featureProximity
- );
+ IOobject
+ (
+ sFeatFileName + ".featureProximity",
+ runTime.constant(),
+ "triSurface",
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ surf,
+ dimLength,
+ featureProximity
+ );
- featureProximityField.write();
+ featureProximityField.write();
- if (writeVTK)
- {
- vtkSurfaceWriter().write
- (
- runTime.constant()/"triSurface", // outputDir
- sFeatFileName, // surfaceName
- surf.points(),
- faces,
- "featureProximity", // fieldName
- featureProximity,
- false, // isNodeValues
- true // verbose
- );
+ if (writeVTK)
+ {
+ vtkSurfaceWriter().write
+ (
+ runTime.constant()/"triSurface", // outputDir
+ sFeatFileName, // surfaceName
+ surf.points(),
+ faces,
+ "featureProximity", // fieldName
+ featureProximity,
+ false, // isNodeValues
+ true // verbose
+ );
+ }
}
+
+ Info<< endl;
}
Info<< "End\n" << endl;
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
new file mode 100644
index 0000000000000000000000000000000000000000..3621d64a293ffa8a1a2daf3be441818036555e92
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
@@ -0,0 +1,86 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object surfaceFeatureExtractDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+surface1.stl
+{
+ // extractFromFile || extractFromSurface
+ extractionMethod extractFromFile;
+
+ extractFromFile
+ {
+ // Load from an existing feature edge file
+ featureEdgeFile "constant/triSurface/featureEdges.nas";
+ }
+
+ trimFeatures
+ {
+ // Remove features with fewer than the specified number of edges
+ minElem 0;
+
+ // Remove features shorter than the specified cumulative length
+ minLen 0.0;
+ }
+
+ subsetFeatures
+ {
+ // Use a plane to select feature edges
+ // (normal)(basePoint)
+ plane (1 0 0)(0 0 0);
+
+ // Select feature edges using a box
+ // (minPt)(maxPt)
+ insideBox (0 0 0)(1 1 1);
+ outsideBox (0 0 0)(1 1 1);
+
+ // Remove any non-manifold (open or > 2 connected faces) edges
+ manifoldEdges no;
+ }
+
+ // Output the curvature of the surface
+ curvature no;
+
+ // Output the proximity of feature points and edges to each other
+ featureProximity no;
+ // The maximum search distance to use when looking for other feature
+ // points and edges
+ maxFeatureProximity 1;
+
+ // Out put the closeness of surface elements to other surface elements.
+ closeness no;
+
+ // Write options
+ writeVTK no;
+ writeObj yes;
+ writeFeatureEdgeMesh no;
+}
+
+
+surface2.nas
+{
+ extractionMethod extractFromSurface;
+
+ extractFromSurface
+ {
+ // Mark edges whose adjacent surface normals are at an angle less
+ // than includedAngle as features
+ // - 0 : selects no edges
+ // - 180: selects all edges
+ includedAngle 120;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
index 024f7d27ecf4203c1684e3295bc7e20ad52692fc..f16dd2595c523183104a52a8330d62588e8f2758 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
@@ -2832,7 +2832,7 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
) const
{
label nearestShapeI = -1;
- point nearestPoint;
+ point nearestPoint = vector::zero;
if (nodes_.size())
{
@@ -2847,10 +2847,6 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
nearestPoint
);
}
- else
- {
- nearestPoint = vector::zero;
- }
return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
}
diff --git a/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C b/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C
index 93196c0e7b537c046089811f91cf5aaeba6d194f..c043ae3c0bb38bba14596fcb06f91afebb4b012e 100644
--- a/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C
+++ b/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C
@@ -113,6 +113,17 @@ bool Foam::fileFormats::NASedgeFormat::read
// discard groupID
dynEdges.append(e);
}
+ else if (cmd == "PLOTEL")
+ {
+ edge e;
+
+ // label groupId = readLabel(IStringStream(line.substr(16,8))());
+ e[0] = readLabel(IStringStream(line.substr(16,8))());
+ e[1] = readLabel(IStringStream(line.substr(24,8))());
+
+ // discard groupID
+ dynEdges.append(e);
+ }
else if (cmd == "GRID")
{
label index = readLabel(IStringStream(line.substr(8,8))());
diff --git a/src/edgeMesh/edgeMesh.C b/src/edgeMesh/edgeMesh.C
index 113112fa288200c568776d3324b77052df627d0e..bc63e3d117a66e28586db39a78cc81e8d289da4e 100644
--- a/src/edgeMesh/edgeMesh.C
+++ b/src/edgeMesh/edgeMesh.C
@@ -28,6 +28,7 @@ License
#include "addToRunTimeSelectionTable.H"
#include "addToMemberFunctionSelectionTable.H"
#include "ListOps.H"
+#include "EdgeMap.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -327,10 +328,7 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
}
// Compact using a hashtable and commutative hash of edge.
- HashTable<label, edge, Hash<edge> > edgeToLabel
- (
- 2*edges_.size()
- );
+ EdgeMap<label> edgeToLabel(2*edges_.size());
label newEdgeI = 0;
@@ -349,13 +347,7 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
edges_.setSize(newEdgeI);
- for
- (
- HashTable<label, edge, Hash<edge> >::const_iterator iter =
- edgeToLabel.begin();
- iter != edgeToLabel.end();
- ++iter
- )
+ forAllConstIter(EdgeMap<label>, edgeToLabel, iter)
{
edges_[iter()] = iter.key();
}
@@ -363,4 +355,35 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
}
+void Foam::edgeMesh::mergeEdges()
+{
+ EdgeMap<label> existingEdges(2*edges_.size());
+
+ label curEdgeI = 0;
+ forAll(edges_, edgeI)
+ {
+ const edge& e = edges_[edgeI];
+
+ if (existingEdges.insert(e, curEdgeI))
+ {
+ curEdgeI++;
+ }
+ }
+
+ if (debug)
+ {
+ Info<< "Merging duplicate edges: "
+ << edges_.size() - existingEdges.size()
+ << " edges will be deleted." << endl;
+ }
+
+ edges_.setSize(existingEdges.size());
+
+ forAllConstIter(EdgeMap<label>, existingEdges, iter)
+ {
+ edges_[iter()] = iter.key();
+ }
+}
+
+
// ************************************************************************* //
diff --git a/src/edgeMesh/edgeMesh.H b/src/edgeMesh/edgeMesh.H
index baa6142815464b4872c5b1f52a09ec927f45741f..66574e9fac6866fac894d3fc9478ea130a270c05 100644
--- a/src/edgeMesh/edgeMesh.H
+++ b/src/edgeMesh/edgeMesh.H
@@ -246,6 +246,9 @@ public:
//- Merge common points (points within mergeDist)
void mergePoints(const scalar mergeDist);
+ //- Merge similar edges
+ void mergeEdges();
+
// Write
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 3652bea89a44b4caf93d16ac97cd4e3ff58eb7f4..0e2de76c07e9afe4c38740a43dfcbf49a9079ec8 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -124,6 +124,7 @@ $(derivedFvPatchFields)/mappedFixedPushedInternalValue/mappedFixedPushedInternal
$(derivedFvPatchFields)/mappedFixedValue/mappedFixedValueFvPatchFields.C
$(derivedFvPatchFields)/mappedVelocityFluxFixedValue/mappedVelocityFluxFixedValueFvPatchField.C
$(derivedFvPatchFields)/mappedFlowRate/mappedFlowRateFvPatchVectorField.C
+$(derivedFvPatchFields)/fixedMean/fixedMeanFvPatchFields.C
$(derivedFvPatchFields)/fan/fanFvPatchFields.C
$(derivedFvPatchFields)/fanPressure/fanPressureFvPatchScalarField.C
$(derivedFvPatchFields)/buoyantPressure/buoyantPressureFvPatchScalarField.C
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C
new file mode 100644
index 0000000000000000000000000000000000000000..e59f98910eb8a1f23cdaf39b7f48f59205b09a73
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixedMeanFvPatchField.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+ const fvPatch& p,
+ const DimensionedField<Type, volMesh>& iF
+)
+:
+ fixedValueFvPatchField<Type>(p, iF),
+ meanValue_(pTraits<Type>::zero)
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+ const fixedMeanFvPatchField<Type>& ptf,
+ const fvPatch& p,
+ const DimensionedField<Type, volMesh>& iF,
+ const fvPatchFieldMapper& mapper
+)
+:
+ fixedValueFvPatchField<Type>(ptf, p, iF, mapper),
+ meanValue_(ptf.meanValue_)
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+ const fvPatch& p,
+ const DimensionedField<Type, volMesh>& iF,
+ const dictionary& dict
+)
+:
+ fixedValueFvPatchField<Type>(p, iF, dict),
+ meanValue_(pTraits<Type>(dict.lookup("meanValue")))
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+ const fixedMeanFvPatchField<Type>& ptf
+)
+:
+ fixedValueFvPatchField<Type>(ptf),
+ meanValue_(ptf.meanValue_)
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+ const fixedMeanFvPatchField<Type>& ptf,
+ const DimensionedField<Type, volMesh>& iF
+)
+:
+ fixedValueFvPatchField<Type>(ptf, iF),
+ meanValue_(ptf.meanValue_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class Type>
+void fixedMeanFvPatchField<Type>::updateCoeffs()
+{
+ if (this->updated())
+ {
+ return;
+ }
+
+ Field<Type> newValues(this->patchInternalField());
+
+ Type meanValuePsi =
+ gSum(this->patch().magSf()*newValues)
+ /gSum(this->patch().magSf());
+
+ if (mag(meanValue_) > SMALL && mag(meanValuePsi)/mag(meanValue_) > 0.5)
+ {
+ newValues *= mag(meanValue_)/mag(meanValuePsi);
+ }
+ else
+ {
+ newValues += (meanValue_ - meanValuePsi);
+ }
+
+ this->operator==(newValues);
+
+ fixedValueFvPatchField<Type>::updateCoeffs();
+}
+
+
+template<class Type>
+void fixedMeanFvPatchField<Type>::write(Ostream& os) const
+{
+ fvPatchField<Type>::write(os);
+ os.writeKeyword("meanValue") << meanValue_ << token::END_STATEMENT << nl;
+ this->writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.H
new file mode 100644
index 0000000000000000000000000000000000000000..b59969a967210d0a31a76751e57fbf2996cb6301
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.H
@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::fixedMeanFvPatchField
+
+Description
+ Extrapolates field to the patch using the near-cell values and adjusts
+ the distribution to match the specified meanValue.
+
+SourceFiles
+ fixedMeanFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedMeanFvPatchField_H
+#define fixedMeanFvPatchField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class fixedMeanFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class fixedMeanFvPatchField
+:
+ public fixedValueFvPatchField<Type>
+{
+
+protected:
+
+ // Protected data
+
+ //- MeanValue value the field is adjusted to maintain
+ Type meanValue_;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("fixedMean");
+
+
+ // Constructors
+
+ //- Construct from patch and internal field
+ fixedMeanFvPatchField
+ (
+ const fvPatch&,
+ const DimensionedField<Type, volMesh>&
+ );
+
+ //- Construct from patch, internal field and dictionary
+ fixedMeanFvPatchField
+ (
+ const fvPatch&,
+ const DimensionedField<Type, volMesh>&,
+ const dictionary&
+ );
+
+ //- Construct by mapping given fixedMeanFvPatchField
+ // onto a new patch
+ fixedMeanFvPatchField
+ (
+ const fixedMeanFvPatchField<Type>&,
+ const fvPatch&,
+ const DimensionedField<Type, volMesh>&,
+ const fvPatchFieldMapper&
+ );
+
+ //- Construct as copy
+ fixedMeanFvPatchField
+ (
+ const fixedMeanFvPatchField<Type>&
+ );
+
+ //- Construct and return a clone
+ virtual tmp<fvPatchField<Type> > clone() const
+ {
+ return tmp<fvPatchField<Type> >
+ (
+ new fixedMeanFvPatchField<Type>(*this)
+ );
+ }
+
+ //- Construct as copy setting internal field reference
+ fixedMeanFvPatchField
+ (
+ const fixedMeanFvPatchField<Type>&,
+ const DimensionedField<Type, volMesh>&
+ );
+
+ //- Construct and return a clone setting internal field reference
+ virtual tmp<fvPatchField<Type> > clone
+ (
+ const DimensionedField<Type, volMesh>& iF
+ ) const
+ {
+ return tmp<fvPatchField<Type> >
+ (
+ new fixedMeanFvPatchField<Type>(*this, iF)
+ );
+ }
+
+
+ // Member functions
+
+ // Evaluation functions
+
+ //- Update the coefficients associated with the patch field
+ virtual void updateCoeffs();
+
+ //- Write
+ virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+# include "fixedMeanFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.C
new file mode 100644
index 0000000000000000000000000000000000000000..e8e269c8a5b4689fdf26d1306f6b5cf0ce59f834
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.C
@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixedMeanFvPatchFields.H"
+#include "volMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makePatchFields(fixedMean);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.H b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..64996c3a682db52980833172f8c9104cc9862fee
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.H
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedMeanFvPatchFields_H
+#define fixedMeanFvPatchFields_H
+
+#include "fixedMeanFvPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeFieldTypedefs(fixedMean);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFieldsFwd.H
similarity index 50%
rename from applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C
rename to src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFieldsFwd.H
index 63e2b904bf5adff04fc1efad940e17d14ad52ffb..9f04dc965bfe98be0cb1e99fcc22a644c4613a1f 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFieldsFwd.H
@@ -23,67 +23,28 @@ License
\*---------------------------------------------------------------------------*/
-#include "buildCGALPolyhedron.H"
+#ifndef fixedMeanFvPatchFieldsFwd_H
+#define fixedMeanFvPatchFieldsFwd_H
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+#include "fieldTypes.H"
-Foam::buildCGALPolyhedron::buildCGALPolyhedron
-(
- const Foam::triSurface& surf
-)
-:
- CGAL::Modifier_base<HalfedgeDS>(),
- surf_(surf)
-{}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-Foam::buildCGALPolyhedron::~buildCGALPolyhedron()
-{}
-
-
-// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
-
-void Foam::buildCGALPolyhedron::operator()
-(
- HalfedgeDS& hds
-)
+namespace Foam
{
- typedef HalfedgeDS::Traits Traits;
- typedef Traits::Point_3 Point;
-
- // Postcondition: `hds' is a valid polyhedral surface.
- CGAL::Polyhedron_incremental_builder_3<HalfedgeDS> B(hds, false);
-
- B.begin_surface
- (
- surf_.points().size(), // n points
- surf_.size(), // n facets
- 2*surf_.edges().size() // n halfedges
- );
- forAll(surf_.points(), pI)
- {
- const Foam::point& p = surf_.points()[pI];
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
- B.add_vertex(Point(p.x(), p.y(), p.z()));
- }
+template<class Type> class fixedMeanFvPatchField;
- forAll(surf_, fI)
- {
- B.begin_facet();
+makePatchTypeFieldTypedefs(fixedMean);
- B.add_vertex_to_facet(surf_[fI][0]);
- B.add_vertex_to_facet(surf_[fI][1]);
- B.add_vertex_to_facet(surf_[fI][2]);
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
- B.end_facet();
- }
+} // End namespace Foam
- B.end_surface();
-}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#endif
// ************************************************************************* //
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
index e39eae5704a61cbc53438060ac1a8ca77154144e..465734ff879bd4aa2507b05fe3779cb1d4984ea9 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
@@ -33,11 +33,15 @@ License
#include "OFstream.H"
#include "IFstream.H"
#include "unitConversion.H"
+#include "EdgeMap.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::surfaceFeatures, 0);
+const Foam::scalar Foam::surfaceFeatures::parallelTolerance =
+ sin(degToRad(1.0));
+
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -170,7 +174,10 @@ void Foam::surfaceFeatures::setFromStatus(const List<edgeStatus>& edgeStat)
//construct feature points where more than 2 feature edges meet
-void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
+void Foam::surfaceFeatures::calcFeatPoints
+(
+ const List<edgeStatus>& edgeStat
+)
{
DynamicList<label> featurePoints(surf_.nPoints()/1000);
@@ -200,6 +207,60 @@ void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
}
+void Foam::surfaceFeatures::classifyFeatureAngles
+(
+ const labelListList& edgeFaces,
+ List<edgeStatus>& edgeStat,
+ const scalar minCos
+) const
+{
+ const vectorField& faceNormals = surf_.faceNormals();
+ const pointField& points = surf_.points();
+
+ forAll(edgeFaces, edgeI)
+ {
+ const labelList& eFaces = edgeFaces[edgeI];
+
+ if (eFaces.size() != 2)
+ {
+ // Non-manifold. What to do here? Is region edge? external edge?
+ edgeStat[edgeI] = REGION;
+ }
+ else
+ {
+ label face0 = eFaces[0];
+ label face1 = eFaces[1];
+
+ if (surf_[face0].region() != surf_[face1].region())
+ {
+ edgeStat[edgeI] = REGION;
+ }
+ else
+ {
+ if ((faceNormals[face0] & faceNormals[face1]) < minCos)
+ {
+
+ // Check if convex or concave by looking at angle
+ // between face centres and normal
+ vector f0Tof1 =
+ surf_[face1].centre(points)
+ - surf_[face0].centre(points);
+
+ if ((f0Tof1 & faceNormals[face0]) > 0.0)
+ {
+ edgeStat[edgeI] = INTERNAL;
+ }
+ else
+ {
+ edgeStat[edgeI] = EXTERNAL;
+ }
+ }
+ }
+ }
+ }
+}
+
+
// Returns next feature edge connected to pointI with correct value.
Foam::label Foam::surfaceFeatures::nextFeatEdge
(
@@ -436,6 +497,80 @@ Foam::surfaceFeatures::surfaceFeatures
}
+Foam::surfaceFeatures::surfaceFeatures
+(
+ const triSurface& surf,
+ const pointField& points,
+ const edgeList& edges,
+ const scalar mergeTol
+)
+:
+ surf_(surf),
+ featurePoints_(0),
+ featureEdges_(0),
+ externalStart_(0),
+ internalStart_(0)
+{
+ // Match edge mesh edges with the triSurface edges
+
+ const labelListList& surfEdgeFaces = surf_.edgeFaces();
+ const edgeList& surfEdges = surf_.edges();
+
+ scalar mergeTolSqr = sqr(mergeTol);
+
+ EdgeMap<label> dynFeatEdges(2*edges.size());
+ DynamicList<labelList> dynFeatureEdgeFaces(edges.size());
+
+ labelList edgeLabel;
+
+ nearestFeatEdge
+ (
+ edges,
+ points,
+ mergeTolSqr,
+ edgeLabel // label of surface edge or -1
+ );
+
+ label count = 0;
+ forAll(edgeLabel, sEdgeI)
+ {
+ const label sEdge = edgeLabel[sEdgeI];
+
+ if (sEdge == -1)
+ {
+ continue;
+ }
+
+ dynFeatEdges.insert(surfEdges[sEdge], count++);
+ dynFeatureEdgeFaces.append(surfEdgeFaces[sEdge]);
+ }
+
+ // Find whether an edge is external or internal
+ List<edgeStatus> edgeStat(dynFeatEdges.size(), NONE);
+
+ classifyFeatureAngles(dynFeatureEdgeFaces, edgeStat, GREAT);
+
+ // Transfer the edge status to a list encompassing all edges in the surface
+ // so that calcFeatPoints can be used.
+ List<edgeStatus> allEdgeStat(surf_.nEdges(), NONE);
+
+ forAll(allEdgeStat, eI)
+ {
+ EdgeMap<label>::const_iterator iter = dynFeatEdges.find(surfEdges[eI]);
+
+ if (iter != dynFeatEdges.end())
+ {
+ allEdgeStat[eI] = edgeStat[iter()];
+ }
+ }
+
+ edgeStat.clear();
+ dynFeatEdges.clear();
+
+ setFromStatus(allEdgeStat);
+}
+
+
//- Construct as copy
Foam::surfaceFeatures::surfaceFeatures(const surfaceFeatures& sf)
:
@@ -496,54 +631,10 @@ void Foam::surfaceFeatures::findFeatures(const scalar includedAngle)
{
scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
- const labelListList& edgeFaces = surf_.edgeFaces();
- const vectorField& faceNormals = surf_.faceNormals();
- const pointField& points = surf_.points();
-
// Per edge whether is feature edge.
List<edgeStatus> edgeStat(surf_.nEdges(), NONE);
- forAll(edgeFaces, edgeI)
- {
- const labelList& eFaces = edgeFaces[edgeI];
-
- if (eFaces.size() != 2)
- {
- // Non-manifold. What to do here? Is region edge? external edge?
- edgeStat[edgeI] = REGION;
- }
- else
- {
- label face0 = eFaces[0];
- label face1 = eFaces[1];
-
- if (surf_[face0].region() != surf_[face1].region())
- {
- edgeStat[edgeI] = REGION;
- }
- else
- {
- if ((faceNormals[face0] & faceNormals[face1]) < minCos)
- {
-
- // Check if convex or concave by looking at angle
- // between face centres and normal
- vector f0Tof1 =
- surf_[face1].centre(points)
- - surf_[face0].centre(points);
-
- if ((f0Tof1 & faceNormals[face0]) > 0.0)
- {
- edgeStat[edgeI] = INTERNAL;
- }
- else
- {
- edgeStat[edgeI] = EXTERNAL;
- }
- }
- }
- }
- }
+ classifyFeatureAngles(surf_.edgeFaces(), edgeStat, minCos);
setFromStatus(edgeStat);
}
@@ -1147,6 +1238,9 @@ void Foam::surfaceFeatures::nearestSurfEdge
const pointField& localPoints = surf_.localPoints();
+ treeBoundBox searchDomain(localPoints);
+ searchDomain.inflate(0.1);
+
indexedOctree<treeDataEdge> ppTree
(
treeDataEdge
@@ -1156,7 +1250,7 @@ void Foam::surfaceFeatures::nearestSurfEdge
localPoints,
selectedEdges
), // all information needed to do geometric checks
- treeBoundBox(localPoints), // overall search domain
+ searchDomain, // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
@@ -1216,6 +1310,8 @@ void Foam::surfaceFeatures::nearestSurfEdge
pointOnEdge.setSize(selectedSampleEdges.size());
pointOnFeature.setSize(selectedSampleEdges.size());
+ treeBoundBox searchDomain(surf_.localPoints());
+
indexedOctree<treeDataEdge> ppTree
(
treeDataEdge
@@ -1224,11 +1320,11 @@ void Foam::surfaceFeatures::nearestSurfEdge
surf_.edges(),
surf_.localPoints(),
selectedEdges
- ), // all information needed to do geometric checks
- treeBoundBox(surf_.localPoints()), // overall search domain
- 8, // maxLevel
- 10, // leafsize
- 3.0 // duplicity
+ ), // all information needed to do geometric checks
+ searchDomain, // overall search domain
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
);
forAll(selectedSampleEdges, i)
@@ -1262,6 +1358,67 @@ void Foam::surfaceFeatures::nearestSurfEdge
}
+void Foam::surfaceFeatures::nearestFeatEdge
+(
+ const edgeList& edges,
+ const pointField& points,
+ scalar searchSpanSqr, // Search span
+ labelList& edgeLabel
+) const
+{
+ edgeLabel = labelList(surf_.nEdges(), -1);
+
+ treeBoundBox searchDomain(points);
+ searchDomain.inflate(0.1);
+
+ indexedOctree<treeDataEdge> ppTree
+ (
+ treeDataEdge
+ (
+ false,
+ edges,
+ points,
+ identity(edges.size())
+ ), // all information needed to do geometric checks
+ searchDomain, // overall search domain
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
+ );
+
+ const edgeList& surfEdges = surf_.edges();
+ const pointField& surfLocalPoints = surf_.localPoints();
+
+ forAll(surfEdges, edgeI)
+ {
+ const edge& sample = surfEdges[edgeI];
+
+ const point& startPoint = surfLocalPoints[sample.start()];
+ const point& midPoint = sample.centre(surfLocalPoints);
+
+ pointIndexHit infoMid = ppTree.findNearest
+ (
+ midPoint,
+ searchSpanSqr
+ );
+
+ if (infoMid.hit())
+ {
+ const vector surfEdgeDir = midPoint - startPoint;
+
+ const edge& featEdge = edges[infoMid.index()];
+ const vector featEdgeDir = featEdge.vec(points);
+
+ // Check that the edges are nearly parallel
+ if (mag(surfEdgeDir ^ featEdgeDir) < parallelTolerance)
+ {
+ edgeLabel[edgeI] = edgeI;
+ }
+ }
+ }
+}
+
+
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void Foam::surfaceFeatures::operator=(const surfaceFeatures& rhs)
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
index 7f7101b662a04046d8888d0b698f09735812141d..0fca9798a2e312c613ace576f1ba12b19c443069 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
@@ -95,6 +95,11 @@ private:
{}
};
+ // Static data
+
+ //- Tolerance for determining whether two vectors are parallel
+ static const scalar parallelTolerance;
+
// Private data
@@ -130,6 +135,14 @@ private:
//- Construct feature points where more than 2 feature edges meet
void calcFeatPoints(const List<edgeStatus>&);
+ //- Classify the angles of the feature edges
+ void classifyFeatureAngles
+ (
+ const labelListList& edgeFaces,
+ List<edgeStatus>& edgeStat,
+ const scalar minCos
+ ) const;
+
//- Choose next unset feature edge.
label nextFeatEdge
(
@@ -186,6 +199,15 @@ public:
//- Construct from file
surfaceFeatures(const triSurface&, const fileName& fName);
+ //- Construct from pointField and edgeList (edgeMesh)
+ surfaceFeatures
+ (
+ const triSurface&,
+ const pointField& points,
+ const edgeList& edges,
+ const scalar mergeTol = 1e-6
+ );
+
//- Construct as copy
surfaceFeatures(const surfaceFeatures&);
@@ -327,9 +349,8 @@ public:
pointField& edgePoint
) const;
-
//- Find nearest surface edge (out of selectedEdges) for each
- // sample edge.
+ // sample edge.
// Sets:
// - edgeLabel : label of surface edge.
// - pointOnEdge : exact position of nearest point on edge.
@@ -347,6 +368,16 @@ public:
pointField& pointOnFeature // point on sample edge
) const;
+ //- Find nearest feature edge to each surface edge. Uses the
+ // mid-point of the surface edges.
+ void nearestFeatEdge
+ (
+ const edgeList& edges,
+ const pointField& points,
+ scalar searchSpanSqr,
+ labelList& edgeLabel
+ ) const;
+
// Write