mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
BUG: meshRefinement: feature edge refinement did not visit all edges
This commit is contained in:
mattijs
@ -623,6 +623,149 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
|
||||
}
|
||||
|
||||
|
||||
// Split faces
|
||||
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitFaces
|
||||
(
|
||||
const labelList& splitFaces,
|
||||
const labelPairList& splits
|
||||
)
|
||||
{
|
||||
polyTopoChange meshMod(mesh_);
|
||||
|
||||
forAll(splitFaces, i)
|
||||
{
|
||||
label faceI = splitFaces[i];
|
||||
const face& f = mesh_.faces()[faceI];
|
||||
|
||||
|
||||
// Split as start and end index in face
|
||||
const labelPair& split = splits[i];
|
||||
|
||||
label nVerts = split[1]-split[0];
|
||||
if (nVerts < 0)
|
||||
{
|
||||
nVerts += f.size();
|
||||
}
|
||||
nVerts += 1;
|
||||
|
||||
|
||||
// Split into f0, f1
|
||||
face f0(nVerts);
|
||||
|
||||
label fp = split[0];
|
||||
forAll(f0, i)
|
||||
{
|
||||
f0[i] = f[fp];
|
||||
fp = f.fcIndex(fp);
|
||||
}
|
||||
|
||||
face f1(f.size()-f0.size()+2);
|
||||
fp = split[1];
|
||||
forAll(f1, i)
|
||||
{
|
||||
f1[i] = f[fp];
|
||||
fp = f.fcIndex(fp);
|
||||
}
|
||||
|
||||
|
||||
// Determine face properties
|
||||
label own = mesh_.faceOwner()[faceI];
|
||||
label nei = -1;
|
||||
label patchI = -1;
|
||||
if (faceI >= mesh_.nInternalFaces())
|
||||
{
|
||||
patchI = mesh_.boundaryMesh().whichPatch(faceI);
|
||||
}
|
||||
else
|
||||
{
|
||||
nei = mesh_.faceNeighbour()[faceI];
|
||||
}
|
||||
|
||||
label zoneI = mesh_.faceZones().whichZone(faceI);
|
||||
bool zoneFlip = false;
|
||||
if (zoneI != -1)
|
||||
{
|
||||
const faceZone& fz = mesh_.faceZones()[zoneI];
|
||||
zoneFlip = fz.flipMap()[fz.whichFace(faceI)];
|
||||
}
|
||||
|
||||
|
||||
Pout<< "face:" << faceI << " verts:" << f
|
||||
<< " split into f0:" << f0
|
||||
<< " f1:" << f1 << endl;
|
||||
|
||||
// Change/add faces
|
||||
meshMod.modifyFace
|
||||
(
|
||||
f0, // modified face
|
||||
faceI, // label of face
|
||||
own, // owner
|
||||
nei, // neighbour
|
||||
false, // face flip
|
||||
patchI, // patch for face
|
||||
zoneI, // zone for face
|
||||
zoneFlip // face flip in zone
|
||||
);
|
||||
meshMod.addFace
|
||||
(
|
||||
f1, // modified face
|
||||
own, // owner
|
||||
nei, // neighbour
|
||||
-1, // master point
|
||||
-1, // master edge
|
||||
faceI, // master face
|
||||
false, // face flip
|
||||
patchI, // patch for face
|
||||
zoneI, // zone for face
|
||||
zoneFlip // face flip in zone
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
|
||||
// Change the mesh (no inflation)
|
||||
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
|
||||
|
||||
// Update fields
|
||||
mesh_.updateMesh(map);
|
||||
|
||||
// Move mesh (since morphing might not do this)
|
||||
if (map().hasMotionPoints())
|
||||
{
|
||||
mesh_.movePoints(map().preMotionPoints());
|
||||
}
|
||||
else
|
||||
{
|
||||
// Delete mesh volumes. No other way to do this?
|
||||
mesh_.clearOut();
|
||||
}
|
||||
|
||||
// Reset the instance for if in overwrite mode
|
||||
mesh_.setInstance(timeName());
|
||||
setInstance(mesh_.facesInstance());
|
||||
|
||||
// Update local mesh data
|
||||
const labelList& oldToNew = map().reverseFaceMap();
|
||||
labelList newSplitFaces(renumber(oldToNew, splitFaces));
|
||||
// Add added faces (every splitFaces becomes two faces)
|
||||
label sz = newSplitFaces.size();
|
||||
newSplitFaces.setSize(2*sz);
|
||||
forAll(map().faceMap(), faceI)
|
||||
{
|
||||
label oldFaceI = map().faceMap()[faceI];
|
||||
if (oldToNew[oldFaceI] != faceI)
|
||||
{
|
||||
// Added face
|
||||
newSplitFaces[sz++] = faceI;
|
||||
}
|
||||
}
|
||||
|
||||
updateMesh(map, newSplitFaces);
|
||||
|
||||
return map;
|
||||
}
|
||||
|
||||
|
||||
//// Determine for multi-processor regions the lowest numbered cell on the
|
||||
//// lowest numbered processor.
|
||||
//void Foam::meshRefinement::getCoupledRegionMaster
|
||||
|
||||
@ -866,6 +866,13 @@ public:
|
||||
const point& keepPoint
|
||||
);
|
||||
|
||||
//- Split faces into two
|
||||
autoPtr<mapPolyMesh> splitFaces
|
||||
(
|
||||
const labelList& splitFaces,
|
||||
const labelPairList& splits
|
||||
);
|
||||
|
||||
//- Update local numbering for mesh redistribution
|
||||
void distribute(const mapDistributePolyMesh&);
|
||||
|
||||
@ -1004,6 +1011,26 @@ public:
|
||||
//- Do any one of above IO functions. flag is combination of
|
||||
// writeFlag values.
|
||||
void write(const label flag, const fileName&) const;
|
||||
|
||||
//- Helper: calculate average
|
||||
template<class T>
|
||||
static T gAverage
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const PackedBoolList& isMasterElem,
|
||||
const UList<T>& values
|
||||
);
|
||||
|
||||
//- Helper: calculate average over selected elements
|
||||
template<class T>
|
||||
static T gAverage
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const PackedBoolList& isMasterElem,
|
||||
const labelList& meshPoints,
|
||||
const UList<T>& values
|
||||
);
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
@ -267,7 +267,12 @@ void Foam::meshRefinement::markFeatureCellLevel
|
||||
|
||||
|
||||
// Find all start cells of features. Is done by tracking from keepPoint.
|
||||
Cloud<trackedParticle> cloud(mesh_, IDLList<trackedParticle>());
|
||||
Cloud<trackedParticle> startPointCloud
|
||||
(
|
||||
mesh_,
|
||||
"startPointCloud",
|
||||
IDLList<trackedParticle>()
|
||||
);
|
||||
|
||||
|
||||
// Features are identical on all processors. Number them so we know
|
||||
@ -315,7 +320,7 @@ void Foam::meshRefinement::markFeatureCellLevel
|
||||
}
|
||||
|
||||
// Non-manifold point. Create particle.
|
||||
cloud.addParticle
|
||||
startPointCloud.addParticle
|
||||
(
|
||||
new trackedParticle
|
||||
(
|
||||
@ -359,7 +364,7 @@ void Foam::meshRefinement::markFeatureCellLevel
|
||||
}
|
||||
|
||||
// Non-manifold point. Create particle.
|
||||
cloud.addParticle
|
||||
startPointCloud.addParticle
|
||||
(
|
||||
new trackedParticle
|
||||
(
|
||||
@ -384,12 +389,13 @@ void Foam::meshRefinement::markFeatureCellLevel
|
||||
maxFeatureLevel = labelList(mesh_.nCells(), -1);
|
||||
|
||||
// Database to pass into trackedParticle::move
|
||||
trackedParticle::trackingData td(cloud, maxFeatureLevel);
|
||||
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.
|
||||
cloud.move(td, GREAT);
|
||||
startPointCloud.move(td, GREAT);
|
||||
|
||||
|
||||
// Reset level
|
||||
maxFeatureLevel = -1;
|
||||
@ -403,8 +409,64 @@ void Foam::meshRefinement::markFeatureCellLevel
|
||||
featureEdgeVisited[featI] = 0u;
|
||||
}
|
||||
|
||||
|
||||
Cloud<trackedParticle> cloud
|
||||
(
|
||||
mesh_,
|
||||
"featureCloud",
|
||||
IDLList<trackedParticle>()
|
||||
);
|
||||
|
||||
forAllIter(Cloud<trackedParticle>, startPointCloud, iter)
|
||||
{
|
||||
const trackedParticle& startTp = iter();
|
||||
|
||||
label featI = startTp.i();
|
||||
label pointI = startTp.j();
|
||||
|
||||
const featureEdgeMesh& featureMesh = features_[featI];
|
||||
const labelList& pEdges = featureMesh.pointEdges()[pointI];
|
||||
|
||||
// Now shoot particles down all pEdges.
|
||||
forAll(pEdges, pEdgeI)
|
||||
{
|
||||
label edgeI = pEdges[pEdgeI];
|
||||
|
||||
if (featureEdgeVisited[featI].set(edgeI, 1u))
|
||||
{
|
||||
// Unvisited edge. Make the particle go to the other point
|
||||
// on the edge.
|
||||
|
||||
const edge& e = featureMesh.edges()[edgeI];
|
||||
label otherPointI = e.otherVertex(pointI);
|
||||
|
||||
trackedParticle* tp(new trackedParticle(startTp));
|
||||
tp->end() = featureMesh.points()[otherPointI];
|
||||
tp->j() = otherPointI;
|
||||
|
||||
if (debug&meshRefinement::FEATURESEEDS)
|
||||
{
|
||||
Pout<< "Adding particle for point:" << pointI
|
||||
<< " coord:" << tp->position()
|
||||
<< " feature:" << featI
|
||||
<< " to track to:" << tp->end()
|
||||
<< endl;
|
||||
}
|
||||
|
||||
cloud.addParticle(tp);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
startPointCloud.clear();
|
||||
|
||||
|
||||
while (true)
|
||||
{
|
||||
// Track all particles to their end position.
|
||||
cloud.move(td, GREAT);
|
||||
|
||||
|
||||
label nParticles = 0;
|
||||
|
||||
// Make particle follow edge.
|
||||
@ -460,10 +522,23 @@ void Foam::meshRefinement::markFeatureCellLevel
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
// Track all particles to their end position.
|
||||
cloud.move(td, GREAT);
|
||||
}
|
||||
|
||||
|
||||
|
||||
//if (debug&meshRefinement::FEATURESEEDS)
|
||||
//{
|
||||
// forAll(maxFeatureLevel, cellI)
|
||||
// {
|
||||
// if (maxFeatureLevel[cellI] != -1)
|
||||
// {
|
||||
// Pout<< "Feature went through cell:" << cellI
|
||||
// << " coord:" << mesh_.cellCentres()[cellI]
|
||||
// << " leve:" << maxFeatureLevel[cellI]
|
||||
// << endl;
|
||||
// }
|
||||
// }
|
||||
//}
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -57,6 +57,108 @@ template<class T> void meshRefinement::updateList
|
||||
}
|
||||
|
||||
|
||||
template<class T>
|
||||
T meshRefinement::gAverage
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const PackedBoolList& isMasterElem,
|
||||
const UList<T>& values
|
||||
)
|
||||
{
|
||||
if (values.size() != isMasterElem.size())
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"meshRefinement::gAverage\n"
|
||||
"(\n"
|
||||
" const polyMesh&,\n"
|
||||
" const PackedBoolList& isMasterElem,\n"
|
||||
" const UList<T>& values\n"
|
||||
")\n"
|
||||
) << "Number of elements in list " << values.size()
|
||||
<< " does not correspond to number of elements in isMasterElem "
|
||||
<< isMasterElem.size()
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
T sum = pTraits<T>::zero;
|
||||
label n = 0;
|
||||
|
||||
forAll(values, i)
|
||||
{
|
||||
if (isMasterElem[i])
|
||||
{
|
||||
sum += values[i];
|
||||
n++;
|
||||
}
|
||||
}
|
||||
|
||||
reduce(sum, sumOp<T>());
|
||||
reduce(n, sumOp<label>());
|
||||
|
||||
if (n > 0)
|
||||
{
|
||||
return sum/n;
|
||||
}
|
||||
else
|
||||
{
|
||||
return pTraits<T>::max;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class T>
|
||||
T meshRefinement::gAverage
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const PackedBoolList& isMasterElem,
|
||||
const labelList& meshPoints,
|
||||
const UList<T>& values
|
||||
)
|
||||
{
|
||||
if (values.size() != meshPoints.size())
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"meshRefinement::gAverage\n"
|
||||
"(\n"
|
||||
" const polyMesh&,\n"
|
||||
" const labelList&,\n"
|
||||
" const PackedBoolList& isMasterElem,\n"
|
||||
" const UList<T>& values\n"
|
||||
")\n"
|
||||
) << "Number of elements in list " << values.size()
|
||||
<< " does not correspond to number of elements in meshPoints "
|
||||
<< meshPoints.size()
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
T sum = pTraits<T>::zero;
|
||||
label n = 0;
|
||||
|
||||
forAll(values, i)
|
||||
{
|
||||
if (isMasterElem[meshPoints[i]])
|
||||
{
|
||||
sum += values[i];
|
||||
n++;
|
||||
}
|
||||
}
|
||||
|
||||
reduce(sum, sumOp<T>());
|
||||
reduce(n, sumOp<label>());
|
||||
|
||||
if (n > 0)
|
||||
{
|
||||
return sum/n;
|
||||
}
|
||||
else
|
||||
{
|
||||
return pTraits<T>::max;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Compare two lists over all boundary faces
|
||||
template<class T>
|
||||
void meshRefinement::testSyncBoundaryFaceList
|
||||
|
||||
Reference in New Issue
Block a user