diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
index a6db1ff0b9c5d00ed487ea8d8485cc945d1d0974..0aa74be6c908020c934f2e7d1818fd59f3bb6332 100644
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
@@ -51,7 +51,6 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
Foam::labelList Foam::medialAxisMeshMover::getFixedValueBCs
(
const pointVectorField& fld
@@ -118,63 +117,6 @@ Foam::medialAxisMeshMover::getPatch
}
-void Foam::medialAxisMeshMover::calculateEdgeWeights
-(
- const PackedBoolList& isMasterEdge,
- const labelList& meshEdges,
- const labelList& meshPoints,
- const edgeList& edges,
- scalarField& edgeWeights,
- scalarField& invSumWeight
-) const
-{
- const pointField& pts = mesh().points();
-
- // Calculate edgeWeights and inverse sum of edge weights
- edgeWeights.setSize(meshEdges.size());
- invSumWeight.setSize(meshPoints.size());
-
- forAll(edges, edgeI)
- {
- const edge& e = edges[edgeI];
- scalar eMag = max
- (
- VSMALL,
- mag
- (
- pts[meshPoints[e[1]]]
- - pts[meshPoints[e[0]]]
- )
- );
- edgeWeights[edgeI] = 1.0/eMag;
- }
-
- // Sum per point all edge weights
- weightedSum
- (
- mesh(),
- isMasterEdge,
- meshEdges,
- meshPoints,
- edges,
- edgeWeights,
- scalarField(meshPoints.size(), 1.0), // data
- invSumWeight
- );
-
- // Inplace invert
- forAll(invSumWeight, pointI)
- {
- scalar w = invSumWeight[pointI];
-
- if (w > 0.0)
- {
- invSumWeight[pointI] = 1.0/w;
- }
- }
-}
-
-
void Foam::medialAxisMeshMover::smoothPatchNormals
(
const label nSmoothDisp,
@@ -193,8 +135,9 @@ void Foam::medialAxisMeshMover::smoothPatchNormals
scalarField edgeWeights(meshEdges.size());
scalarField invSumWeight(meshPoints.size());
- calculateEdgeWeights
+ meshRefinement::calculateEdgeWeights
(
+ mesh(),
isMasterEdge,
meshEdges,
meshPoints,
@@ -207,7 +150,7 @@ void Foam::medialAxisMeshMover::smoothPatchNormals
vectorField average;
for (label iter = 0; iter < nSmoothDisp; iter++)
{
- weightedSum
+ meshRefinement::weightedSum
(
mesh(),
isMasterEdge,
@@ -284,8 +227,9 @@ void Foam::medialAxisMeshMover::smoothNormals
scalarField edgeWeights(meshEdges.size());
scalarField invSumWeight(meshPoints.size());
- calculateEdgeWeights
+ meshRefinement::calculateEdgeWeights
(
+ mesh(),
isMasterEdge,
meshEdges,
meshPoints,
@@ -297,7 +241,7 @@ void Foam::medialAxisMeshMover::smoothNormals
vectorField average;
for (label iter = 0; iter < nSmoothDisp; iter++)
{
- weightedSum
+ meshRefinement::weightedSum
(
mesh(),
isMasterEdge,
@@ -1059,8 +1003,9 @@ void Foam::medialAxisMeshMover::minSmoothField
scalarField edgeWeights(meshEdges.size());
scalarField invSumWeight(meshPoints.size());
- calculateEdgeWeights
+ meshRefinement::calculateEdgeWeights
(
+ mesh(),
isMasterEdge,
meshEdges,
meshPoints,
@@ -1075,7 +1020,7 @@ void Foam::medialAxisMeshMover::minSmoothField
for (label iter = 0; iter < nSmoothDisp; iter++)
{
scalarField average(pp.nPoints());
- weightedSum
+ meshRefinement::weightedSum
(
mesh(),
isMasterEdge,
@@ -1534,8 +1479,9 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
// Calculate inverse sum of weights
scalarField edgeWeights(meshEdges.size());
scalarField invSumWeight(meshPoints.size());
- calculateEdgeWeights
+ meshRefinement::calculateEdgeWeights
(
+ mesh(),
isMasterEdge,
meshEdges,
meshPoints,
@@ -1554,7 +1500,7 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
for (label iter = 0; iter < nSmoothDisp; iter++)
{
- weightedSum
+ meshRefinement::weightedSum
(
mesh(),
isMasterEdge,
@@ -1575,7 +1521,7 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
}
}
- weightedSum
+ meshRefinement::weightedSum
(
mesh(),
isMasterEdge,
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
index b402ecba259731069ccc03c8f02b8be9a87f2fb4..1c270de274fd857db5df10b3e594e1acba7fe876 100644
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
@@ -105,31 +105,6 @@ class medialAxisMeshMover
const labelList&
);
- //- Weighted sum (over all subset of mesh points) by
- // summing contribution from (master) edges
- template<class Type>
- static void weightedSum
- (
- const polyMesh& mesh,
- const PackedBoolList& isMasterEdge,
- const labelList& meshEdges,
- const labelList& meshPoints,
- const edgeList& edges,
- const scalarField& edgeWeights,
- const Field<Type>& data,
- Field<Type>& sum
- );
-
- void calculateEdgeWeights
- (
- const PackedBoolList& isMasterEdge,
- const labelList& meshEdges,
- const labelList& meshPoints,
- const edgeList& edges,
- scalarField& edgeWeights,
- scalarField& invSumWeight
- ) const;
-
// Calculation of medial axis information
@@ -314,12 +289,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#ifdef NoRepository
-# include "medialAxisMeshMoverTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
#endif
// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C
deleted file mode 100644
index edd053f63ff270de357ec0ca2d7c253ab51fb038..0000000000000000000000000000000000000000
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C
+++ /dev/null
@@ -1,93 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "medialAxisMeshMover.H"
-#include "syncTools.H"
-
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-template<class Type>
-void Foam::medialAxisMeshMover::weightedSum
-(
- const polyMesh& mesh,
- const PackedBoolList& isMasterEdge,
- const labelList& meshEdges,
- const labelList& meshPoints,
- const edgeList& edges,
- const scalarField& edgeWeights,
- const Field<Type>& pointData,
- Field<Type>& sum
-)
-{
- if
- (
- mesh.nEdges() != isMasterEdge.size()
- || edges.size() != meshEdges.size()
- || edges.size() != edgeWeights.size()
- || meshPoints.size() != pointData.size()
- )
- {
- FatalErrorIn("medialAxisMeshMover::weightedSum(..)")
- << "Inconsistent sizes for edge or point data:"
- << " meshEdges:" << meshEdges.size()
- << " isMasterEdge:" << isMasterEdge.size()
- << " edgeWeights:" << edgeWeights.size()
- << " edges:" << edges.size()
- << " pointData:" << pointData.size()
- << " meshPoints:" << meshPoints.size()
- << abort(FatalError);
- }
-
- sum.setSize(meshPoints.size());
- sum = pTraits<Type>::zero;
-
- forAll(edges, edgeI)
- {
- if (isMasterEdge.get(meshEdges[edgeI]) == 1)
- {
- const edge& e = edges[edgeI];
-
- scalar eWeight = edgeWeights[edgeI];
-
- label v0 = e[0];
- label v1 = e[1];
-
- sum[v0] += eWeight*pointData[v1];
- sum[v1] += eWeight*pointData[v0];
- }
- }
-
- syncTools::syncPointList
- (
- mesh,
- meshPoints,
- sum,
- plusEqOp<Type>(),
- pTraits<Type>::zero // null value
- );
-}
-
-
-// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index d93a5ca6d37b501bf12a060f503b2adda2c2d2e1..8d41eebd2540a279e5c3c95d5154624da9808725 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -1949,6 +1949,64 @@ void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
}
+void Foam::meshRefinement::calculateEdgeWeights
+(
+ const polyMesh& mesh,
+ const PackedBoolList& isMasterEdge,
+ const labelList& meshEdges,
+ const labelList& meshPoints,
+ const edgeList& edges,
+ scalarField& edgeWeights,
+ scalarField& invSumWeight
+)
+{
+ const pointField& pts = mesh.points();
+
+ // Calculate edgeWeights and inverse sum of edge weights
+ edgeWeights.setSize(meshEdges.size());
+ invSumWeight.setSize(meshPoints.size());
+
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+ scalar eMag = max
+ (
+ VSMALL,
+ mag
+ (
+ pts[meshPoints[e[1]]]
+ - pts[meshPoints[e[0]]]
+ )
+ );
+ edgeWeights[edgeI] = 1.0/eMag;
+ }
+
+ // Sum per point all edge weights
+ weightedSum
+ (
+ mesh,
+ isMasterEdge,
+ meshEdges,
+ meshPoints,
+ edges,
+ edgeWeights,
+ scalarField(meshPoints.size(), 1.0), // data
+ invSumWeight
+ );
+
+ // Inplace invert
+ forAll(invSumWeight, pointI)
+ {
+ scalar w = invSumWeight[pointI];
+
+ if (w > 0.0)
+ {
+ invSumWeight[pointI] = 1.0/w;
+ }
+ }
+}
+
+
Foam::label Foam::meshRefinement::appendPatch
(
fvMesh& mesh,
@@ -2167,9 +2225,7 @@ Foam::label Foam::meshRefinement::addMeshedPatch
// );
// Store
- label sz = meshedPatches_.size();
- meshedPatches_.setSize(sz+1);
- meshedPatches_[sz] = name;
+ meshedPatches_.append(name);
return patchI;
}
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index f8d6258ba7c88deb888272674a583525308a19ba..57114bd7404e6358535df9f598a2de00f01fe84f 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -729,6 +729,33 @@ public:
//- Helper function: check that face zones are synced
static void checkCoupledFaceZones(const polyMesh&);
+ //- Helper: calculate edge weights (1/length)
+ static void calculateEdgeWeights
+ (
+ const polyMesh& mesh,
+ const PackedBoolList& isMasterEdge,
+ const labelList& meshEdges,
+ const labelList& meshPoints,
+ const edgeList& edges,
+ scalarField& edgeWeights,
+ scalarField& invSumWeight
+ );
+
+ //- Helper: weighted sum (over all subset of mesh points) by
+ // summing contribution from (master) edges
+ template<class Type>
+ static void weightedSum
+ (
+ const polyMesh& mesh,
+ const PackedBoolList& isMasterEdge,
+ const labelList& meshEdges,
+ const labelList& meshPoints,
+ const edgeList& edges,
+ const scalarField& edgeWeights,
+ const Field<Type>& data,
+ Field<Type>& sum
+ );
+
// Refinement
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index a7084e9fc6ef204bd6e1c658ea5bf617e7fe44cd..701c0dbbc077d1f58f5b92cff33265c6cefb711f 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -458,10 +458,18 @@ void Foam::meshRefinement::markFeatureCellLevel
// Database to pass into trackedParticle::move
trackedParticle::trackingData td(startPointCloud, maxFeatureLevel);
+
// Track all particles to their end position (= starting feature point)
// Note that the particle might have started on a different processor
// so this will transfer across nicely until we can start tracking proper.
scalar maxTrackLen = 2.0*mesh_.bounds().mag();
+
+ if (debug&meshRefinement::FEATURESEEDS)
+ {
+ Pout<< "Tracking " << startPointCloud.size()
+ << " particles over distance " << maxTrackLen
+ << " to find the starting cell" << endl;
+ }
startPointCloud.move(td, maxTrackLen);
@@ -485,6 +493,11 @@ void Foam::meshRefinement::markFeatureCellLevel
IDLList<trackedParticle>()
);
+ if (debug&meshRefinement::FEATURESEEDS)
+ {
+ Pout<< "Constructing cloud for cell marking" << endl;
+ }
+
forAllIter(Cloud<trackedParticle>, startPointCloud, iter)
{
const trackedParticle& startTp = iter();
@@ -532,11 +545,14 @@ void Foam::meshRefinement::markFeatureCellLevel
while (true)
{
// Track all particles to their end position.
+ if (debug&meshRefinement::FEATURESEEDS)
+ {
+ Pout<< "Tracking " << cloud.size()
+ << " particles over distance " << maxTrackLen
+ << " to mark cells" << endl;
+ }
cloud.move(td, maxTrackLen);
-
- label nParticles = 0;
-
// Make particle follow edge.
forAllIter(Cloud<trackedParticle>, cloud, iter)
{
@@ -578,15 +594,15 @@ void Foam::meshRefinement::markFeatureCellLevel
// seeded. Delete particle.
cloud.deleteParticle(tp);
}
- else
- {
- // Keep particle
- nParticles++;
- }
}
- reduce(nParticles, sumOp<label>());
- if (nParticles == 0)
+
+ if (debug&meshRefinement::FEATURESEEDS)
+ {
+ Pout<< "Remaining particles " << cloud.size() << endl;
+ }
+
+ if (returnReduce(cloud.size(), sumOp<label>()) == 0)
{
break;
}
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
index 45fcee4d5d793861002b045680ef7763618c8005..528e812197f2da2fd1870d9a5b43bdb04be5042e 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
@@ -26,16 +26,12 @@ License
#include "meshRefinement.H"
#include "fvMesh.H"
#include "globalIndex.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
+#include "syncTools.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Add a T entry
-template<class T> void meshRefinement::updateList
+template<class T> void Foam::meshRefinement::updateList
(
const labelList& newToOld,
const T& nullValue,
@@ -59,7 +55,7 @@ template<class T> void meshRefinement::updateList
template<class T>
-T meshRefinement::gAverage
+T Foam::meshRefinement::gAverage
(
const polyMesh& mesh,
const PackedBoolList& isMasterElem,
@@ -109,7 +105,7 @@ T meshRefinement::gAverage
template<class T>
-T meshRefinement::gAverage
+T Foam::meshRefinement::gAverage
(
const polyMesh& mesh,
const PackedBoolList& isMasterElem,
@@ -162,7 +158,7 @@ T meshRefinement::gAverage
// Compare two lists over all boundary faces
template<class T>
-void meshRefinement::testSyncBoundaryFaceList
+void Foam::meshRefinement::testSyncBoundaryFaceList
(
const scalar tol,
const string& msg,
@@ -222,7 +218,7 @@ void meshRefinement::testSyncBoundaryFaceList
// Print list sorted by coordinates. Used for comparing non-parallel v.s.
// parallel operation
template<class T>
-void meshRefinement::collectAndPrint
+void Foam::meshRefinement::collectAndPrint
(
const UList<point>& points,
const UList<T>& data
@@ -268,7 +264,11 @@ void meshRefinement::collectAndPrint
//template<class T, class Mesh>
template<class GeoField>
-void meshRefinement::addPatchFields(fvMesh& mesh, const word& patchFieldType)
+void Foam::meshRefinement::addPatchFields
+(
+ fvMesh& mesh,
+ const word& patchFieldType
+)
{
HashTable<GeoField*> flds
(
@@ -298,7 +298,11 @@ void meshRefinement::addPatchFields(fvMesh& mesh, const word& patchFieldType)
// Reorder patch field
template<class GeoField>
-void meshRefinement::reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
+void Foam::meshRefinement::reorderPatchFields
+(
+ fvMesh& mesh,
+ const labelList& oldToNew
+)
{
HashTable<GeoField*> flds
(
@@ -316,7 +320,11 @@ void meshRefinement::reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
template<class Enum>
-int meshRefinement::readFlags(const Enum& namedEnum, const wordList& words)
+int Foam::meshRefinement::readFlags
+(
+ const Enum& namedEnum,
+ const wordList& words
+)
{
int flags = 0;
@@ -330,8 +338,66 @@ int meshRefinement::readFlags(const Enum& namedEnum, const wordList& words)
}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+template<class Type>
+void Foam::meshRefinement::weightedSum
+(
+ const polyMesh& mesh,
+ const PackedBoolList& isMasterEdge,
+ const labelList& meshEdges,
+ const labelList& meshPoints,
+ const edgeList& edges,
+ const scalarField& edgeWeights,
+ const Field<Type>& pointData,
+ Field<Type>& sum
+)
+{
+ if
+ (
+ mesh.nEdges() != isMasterEdge.size()
+ || edges.size() != meshEdges.size()
+ || edges.size() != edgeWeights.size()
+ || meshPoints.size() != pointData.size()
+ )
+ {
+ FatalErrorIn("medialAxisMeshMover::weightedSum(..)")
+ << "Inconsistent sizes for edge or point data:"
+ << " meshEdges:" << meshEdges.size()
+ << " isMasterEdge:" << isMasterEdge.size()
+ << " edgeWeights:" << edgeWeights.size()
+ << " edges:" << edges.size()
+ << " pointData:" << pointData.size()
+ << " meshPoints:" << meshPoints.size()
+ << abort(FatalError);
+ }
+
+ sum.setSize(meshPoints.size());
+ sum = pTraits<Type>::zero;
+
+ forAll(edges, edgeI)
+ {
+ if (isMasterEdge.get(meshEdges[edgeI]) == 1)
+ {
+ const edge& e = edges[edgeI];
+
+ scalar eWeight = edgeWeights[edgeI];
+
+ label v0 = e[0];
+ label v1 = e[1];
+
+ sum[v0] += eWeight*pointData[v1];
+ sum[v1] += eWeight*pointData[v0];
+ }
+ }
+
+ syncTools::syncPointList
+ (
+ mesh,
+ meshPoints,
+ sum,
+ plusEqOp<Type>(),
+ pTraits<Type>::zero // null value
+ );
+}
-} // End namespace Foam
// ************************************************************************* //