exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Commit details:

Browse source 
- Added initial support for parallel dynamicTopoFvMesh
 - Conservative remapping support (serial / parallel)
 - Layer addition/removal support (partially tested / not parallel)
 - Dynamic parallel load-balancing support
This commit is contained in:
Sandeep Menon 2011-05-17 15:31:43 -04:00
parent 4a271dbf06
commit dc43e9152c
52 changed files with 30104 additions and 4971 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

48
src/dynamicMesh/dynamicFvMesh/Make/files
View file

@ -20,6 +20,38 @@ $(solidBodyMotionFunctions)/translation/translation.C
mixerGgiFvMesh/mixerGgiFvMesh.C
turboFvMesh/turboFvMesh.C
eMesh = dynamicTopoFvMesh/eMesh
$(eMesh)/eMesh.C
$(eMesh)/eMeshDemandDrivenData.C
$(eMesh)/eBoundaryMesh/eBoundaryMesh.C
ePatches = $(eMesh)/ePatches
$(ePatches)/ePatch/ePatch.C
$(ePatches)/ePatch/newEPatch.C
dynamicTopoFvMesh/dynamicTopoFvMesh.C
dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C
dynamicTopoFvMesh/dynamicTopoFvMeshReOrder.C
dynamicTopoFvMesh/dynamicTopoFvMeshCoupled.C
dynamicTopoFvMesh/dynamicTopoFvMeshMapping.C
dynamicTopoFvMesh/edgeSwap.C
dynamicTopoFvMesh/edgeBisect.C
dynamicTopoFvMesh/edgeCollapse.C
dynamicTopoFvMesh/coupleMap.C
convexSetAlgorithm = dynamicTopoFvMesh/convexSetAlgorithm
$(convexSetAlgorithm)/convexSetAlgorithm.C
$(convexSetAlgorithm)/faceSetAlgorithm.C
$(convexSetAlgorithm)/cellSetAlgorithm.C
fieldMapping = dynamicTopoFvMesh/fieldMapping
$(fieldMapping)/topoMapper.C
$(fieldMapping)/fluxCorrector.C
$(fieldMapping)/topoCellMapper.C
$(fieldMapping)/topoPatchMapper.C
$(fieldMapping)/topoSurfaceMapper.C
tetMetrics = dynamicTopoFvMesh/tetMetrics
$(tetMetrics)/tetMetric.C
$(tetMetrics)/tetMetrics.C
@ -27,20 +59,4 @@ $(tetMetrics)/tetMetrics.C
lengthScaleEstimator = dynamicTopoFvMesh/lengthScaleEstimator
$(lengthScaleEstimator)/lengthScaleEstimator.C
eMesh = dynamicTopoFvMesh/eMesh
$(eMesh)/eMesh.C
$(eMesh)/eMeshDemandDrivenData.C
$(eMesh)/eBoundaryMesh/eBoundaryMesh.C
ePatches = $(eMesh)/ePatches
$(ePatches)/ePatch/ePatch.C
$(ePatches)/ePatch/newEPatch.C
dynamicTopoFvMesh/meshOps.C
dynamicTopoFvMesh/dynamicTopoFvMesh.C
dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C
dynamicTopoFvMesh/dynamicTopoFvMeshReOrder.C
dynamicTopoFvMesh/dynamicTopoFvMeshMapping.C
dynamicTopoFvMesh/edgeSwap.C
dynamicTopoFvMesh/edgeBisect.C
dynamicTopoFvMesh/edgeCollapse.C
LIB = $(FOAM_LIBBIN)/libdynamicFvMesh

2
src/dynamicMesh/dynamicFvMesh/Make/options
View file

@ -3,6 +3,7 @@ EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/decompositionMethods/decompositionMethods/lnInclude \
-I$(LIB_SRC)/dynamicMesh/meshMotion/tetDecompositionMotionSolver/lnInclude \
-I$(LIB_SRC)/tetDecompositionFiniteElement/lnInclude \
$(WM_DECOMP_INC) \
@ -15,6 +16,7 @@ LIB_LIBS = \
-lmeshTools \
-ldynamicMesh \
-lfiniteVolume \
-ldecompositionMethods \
$(WM_DECOMP_LIBS) \
-lfvMotionSolver \
-lRBFMotionSolver \

123
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/changeMap.H
View file

@ -41,109 +41,140 @@ SourceFiles
#ifndef changeMap_H
#define changeMap_H
#include "objectMap.H"
#include "labelList.H"
#include "FixedList.H"
#include "dictionary.H"
#include "DynamicList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
class changeMap;
// * * * * * * * * Forward declaration of friend fuctions * * * * * * * * * //
Ostream& operator<<(Ostream&, const changeMap&);
/*---------------------------------------------------------------------------*\
Class changeMap Declaration
\*---------------------------------------------------------------------------*/
class changeMap
:
public dictionary
{
// Sliver type index.
// Type 1: Sliver
// Type 2: Cap
// Type 3: Spade
// Type 4: Wedge
// Entity index
label index_;
// Coupled patch index
label pIndex_;
// Type
label type_;
// Data specific to sliver-type cells
label firstEdge_;
label secondEdge_;
// Data specific to cap-type cells
label apexPoint_;
label opposingFace_;
// Entities that were added during the operation.
// - Add a master mapping for slaves.
// - Masters get a map of -1.
List<FixedList<label,2> > addedPoints_;
List<FixedList<label,2> > addedEdges_;
List<FixedList<label,2> > addedFaces_;
List<FixedList<label,2> > addedCells_;
DynamicList<objectMap> addedPoints_;
DynamicList<objectMap> addedEdges_;
DynamicList<objectMap> addedFaces_;
DynamicList<objectMap> addedCells_;
// Entities that were removed during the operation
DynamicList<label> removedPoints_;
DynamicList<label> removedEdges_;
DynamicList<label> removedFaces_;
DynamicList<label> removedCells_;
public:
// Constructor
changeMap()
:
index_(-1),
pIndex_(-1),
type_(-1),
firstEdge_(-1),
secondEdge_(-1),
apexPoint_(-1),
opposingFace_(-1),
addedPoints_(0),
addedEdges_(0),
addedFaces_(0),
addedCells_(0)
addedPoints_(5),
addedEdges_(5),
addedFaces_(5),
addedCells_(5),
removedPoints_(5),
removedEdges_(5),
removedFaces_(5),
removedCells_(5)
{}
//- Access
// Entity index
inline label& index();
inline label index() const;
// Coupled patch index
inline label& patchIndex();
inline label patchIndex() const;
// Type
inline label& type();
inline label type() const;
// For sliver-type cells, opposite edges
// are identified for removal.
inline label& firstEdge();
inline label& secondEdge();
// For cap-type cells, the face requiring splitting
// is identified for removal.
inline label& apexPoint();
inline label& opposingFace();
// Added entities
inline void addPoint
(
const label pIndex,
const label master = -1
const labelList& master = labelList()
);
inline void addEdge
(
const label eIndex,
const label master = -1
const labelList& master = labelList()
);
inline void addFace
(
const label fIndex,
const label master = -1
const labelList& master = labelList()
);
inline void addCell
(
const label cIndex,
const label master = -1
const labelList& master = labelList()
);
// Return the list of added entities
inline const List<FixedList<label,2> >& addedPointList() const;
inline const List<FixedList<label,2> >& addedEdgeList() const;
inline const List<FixedList<label,2> >& addedFaceList() const;
inline const List<FixedList<label,2> >& addedCellList() const;
inline const List<objectMap>& addedPointList() const;
inline const List<objectMap>& addedEdgeList() const;
inline const List<objectMap>& addedFaceList() const;
inline const List<objectMap>& addedCellList() const;
// Removed entities
inline void removePoint(const label pIndex);
inline void removeEdge(const label eIndex);
inline void removeFace(const label fIndex);
inline void removeCell(const label cIndex);
// Return the list of removed entities
inline const labelList& removedPointList() const;
inline const labelList& removedEdgeList() const;
inline const labelList& removedFaceList() const;
inline const labelList& removedCellList() const;
//- Edit
// Clear existing lists
inline void clear();
//- Operators
inline void operator=(const changeMap& rhs);
//- IOstream Operators
inline friend Ostream& operator<<(Ostream&, const changeMap&);
};

280
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/changeMapI.H
View file

@ -22,6 +22,17 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
changeMap
Description
Accumulate topology change statistics
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
namespace Foam
@ -29,6 +40,32 @@ namespace Foam
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Entity index
inline label& changeMap::index()
{
return index_;
}
inline label changeMap::index() const
{
return index_;
}
// Coupled patch index
inline label& changeMap::patchIndex()
{
return pIndex_;
}
inline label changeMap::patchIndex() const
{
return pIndex_;
}
// Type
inline label& changeMap::type()
{
@ -42,140 +79,241 @@ inline label changeMap::type() const
}
// For sliver-type cells, opposite edges
// are identified for removal.
inline label& changeMap::firstEdge()
{
return firstEdge_;
}
inline label& changeMap::secondEdge()
{
return secondEdge_;
}
// For cap-type cells, the face requiring splitting
// is identified for removal.
inline label& changeMap::apexPoint()
{
return apexPoint_;
}
inline label& changeMap::opposingFace()
{
return opposingFace_;
}
// Added entities
inline void changeMap::addPoint
(
const label pIndex,
const label master
const labelList& master
)
{
label curSize = addedPoints_.size();
addedPoints_.setSize(curSize + 1);
addedPoints_[curSize][0] = pIndex;
addedPoints_[curSize][1] = master;
addedPoints_.append(objectMap(pIndex, master));
}
inline void changeMap::addEdge
(
const label eIndex,
const label master
const labelList& master
)
{
label curSize = addedEdges_.size();
addedEdges_.setSize(curSize + 1);
addedEdges_[curSize][0] = eIndex;
addedEdges_[curSize][1] = master;
addedEdges_.append(objectMap(eIndex, master));
}
inline void changeMap::addFace
(
const label fIndex,
const label master
const labelList& master
)
{
label curSize = addedFaces_.size();
addedFaces_.setSize(curSize + 1);
addedFaces_[curSize][0] = fIndex;
addedFaces_[curSize][1] = master;
addedFaces_.append(objectMap(fIndex, master));
}
inline void changeMap::addCell
(
const label cIndex,
const label master
const labelList& master
)
{
label curSize = addedCells_.size();
addedCells_.setSize(curSize + 1);
addedCells_[curSize][0] = cIndex;
addedCells_[curSize][1] = master;
addedCells_.append(objectMap(cIndex, master));
}
// Return an added point
inline const List<FixedList<label,2> >&
inline const List<objectMap>&
changeMap::addedPointList() const
{
return addedPoints_;
}
// Return the list of added entities
inline const List<FixedList<label,2> >&
inline const List<objectMap>&
changeMap::addedEdgeList() const
{
return addedEdges_;
}
inline const List<FixedList<label,2> >&
inline const List<objectMap>&
changeMap::addedFaceList() const
{
return addedFaces_;
}
inline const List<FixedList<label,2> >&
inline const List<objectMap>&
changeMap::addedCellList() const
{
return addedCells_;
}
inline void changeMap::removePoint(const label pIndex)
{
removedPoints_.append(pIndex);
}
inline void changeMap::removeEdge(const label eIndex)
{
removedEdges_.append(eIndex);
}
inline void changeMap::removeFace(const label fIndex)
{
removedFaces_.append(fIndex);
}
inline void changeMap::removeCell(const label cIndex)
{
removedCells_.append(cIndex);
}
inline const labelList&
changeMap::removedPointList() const
{
return removedPoints_;
}
inline const labelList&
changeMap::removedEdgeList() const
{
return removedEdges_;
}
inline const labelList&
changeMap::removedFaceList() const
{
return removedFaces_;
}
inline const labelList&
changeMap::removedCellList() const
{
return removedCells_;
}
// Clear existing lists
inline void changeMap::clear()
{
// Clear base dictionary
dictionary::clear();
index_ = -1;
pIndex_ = -1;
type_ = -1;
// Clear added entities
addedPoints_.clear();
addedEdges_.clear();
addedFaces_.clear();
addedCells_.clear();
// Clear removed entities
removedPoints_.clear();
removedEdges_.clear();
removedFaces_.clear();
removedCells_.clear();
}
inline void changeMap::operator=(const changeMap& rhs)
{
// Copy base dictionary
dictionary::operator=(rhs);
index_ = rhs.index_;
pIndex_ = rhs.pIndex_;
type_ = rhs.type_;
firstEdge_ = rhs.firstEdge_;
secondEdge_ = rhs.secondEdge_;
apexPoint_ = rhs.apexPoint_;
opposingFace_ = rhs.opposingFace_;
// Copy maps
addedPoints_.setSize(rhs.addedPoints_.size());
addedPoints_ = rhs.addedPoints_;
addedEdges_ = rhs.addedEdges_;
addedFaces_ = rhs.addedFaces_;
addedCells_ = rhs.addedCells_;
forAll(addedPoints_, indexI)
{
addedPoints_[indexI].index() = rhs.addedPoints_[indexI].index();
addedPoints_[indexI].masterObjects() =
(
rhs.addedPoints_[indexI].masterObjects()
);
}
addedEdges_.setSize(rhs.addedEdges_.size());
forAll(addedEdges_, indexI)
{
addedEdges_[indexI].index() = rhs.addedEdges_[indexI].index();
addedEdges_[indexI].masterObjects() =
(
rhs.addedEdges_[indexI].masterObjects()
);
}
addedFaces_.setSize(rhs.addedFaces_.size());
forAll(addedFaces_, indexI)
{
addedFaces_[indexI].index() = rhs.addedFaces_[indexI].index();
addedFaces_[indexI].masterObjects() =
(
rhs.addedFaces_[indexI].masterObjects()
);
}
addedCells_.setSize(rhs.addedCells_.size());
forAll(addedCells_, indexI)
{
addedCells_[indexI].index() = rhs.addedCells_[indexI].index();
addedCells_[indexI].masterObjects() =
(
rhs.addedCells_[indexI].masterObjects()
);
}
removedPoints_ = rhs.removedPoints_;
removedEdges_ = rhs.removedEdges_;
removedFaces_ = rhs.removedFaces_;
removedCells_ = rhs.removedCells_;
}
inline Ostream& operator<<(Ostream& os, const changeMap& cm)
{
// Write base dictionary
os << static_cast<const dictionary&>(cm) << endl;
os << " index: " << cm.index_ << nl
<< " patchIndex: " << cm.pIndex_ << nl
<< " type: " << cm.type_ << nl
<< " addedPoints: " << cm.addedPoints_ << nl
<< " addedEdges: " << cm.addedEdges_ << nl
<< " addedFaces: " << cm.addedFaces_ << nl
<< " addedCells: " << cm.addedCells_ << nl
<< " removedPoints: " << cm.removedPoints_ << nl
<< " removedEdges: " << cm.removedEdges_ << nl
<< " removedFaces: " << cm.removedFaces_ << nl
<< " removedCells: " << cm.removedCells_ << endl;
// Check state of IOstream
os.check("Ostream& operator<<(Ostream&, changeMap&)");
return os;
}
} // End namespace Foam
// ************************************************************************* //

470
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/cellSetAlgorithm.C Normal file
View file

@ -0,0 +1,470 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
cellSetAlgorithm
Description
Class for convexSetAlgorithms pertaining to cells
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#include "cellSetAlgorithm.H"
#include "meshOps.H"
#include "triFace.H"
#include "polyMesh.H"
#include "tetIntersection.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
cellSetAlgorithm::cellSetAlgorithm
(
const polyMesh& mesh,
const pointField& newPoints,
const UList<edge>& newEdges,
const UList<face>& newFaces,
const UList<cell>& newCells,
const UList<label>& newOwner,
const UList<label>& newNeighbour
)
:
convexSetAlgorithm
(
mesh,
newPoints,
newEdges,
newFaces,
newCells,
newOwner,
newNeighbour
)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void cellSetAlgorithm::computeNormFactor(const label index) const
{
// Clear existing fields
parents_.clear();
centres_.clear();
weights_.clear();
// Compute volume / centre (using refNorm_ as centre)
meshOps::cellCentreAndVolume
(
index,
newPoints_,
newFaces_,
newCells_,
newOwner_,
refNorm_,
normFactor_
);
// Compute a bounding box around the cell
box_ = boundBox(newCells_[index].points(newFaces_, newPoints_), false);
vector minToXb = (box_.min() - box_.midpoint());
vector maxToXb = (box_.max() - box_.midpoint());
// Scale it by a bit
box_.min() += (1.5 * minToXb);
box_.max() += (1.5 * maxToXb);
}
// Check whether the bounding box contains the entity
bool cellSetAlgorithm::contains(const label index) const
{
// Fetch old points
const labelListList& cellPoints = mesh_.cellPoints();
const labelList& checkCell = cellPoints[index];
// Check if the bounding box contains any of the supplied points
forAll(checkCell, pointI)
{
if (box_.contains(mesh_.points()[checkCell[pointI]]))
{
return true;
}
}
return false;
}
// Compute intersection
bool cellSetAlgorithm::computeIntersection
(
const label newIndex,
const label oldIndex,
const label offset,
bool output
) const
{
bool intersects = false;
const pointField& newPoints = newPoints_;
const pointField& oldPoints = mesh_.points();
const cell& newCell = newCells_[newIndex];
const cell& oldCell = mesh_.cells()[oldIndex];
// Check if decomposition is necessary
if (oldCell.size() > 4 || newCell.size() > 4)
{
// Decompose new / old cells
DynamicList<FixedList<point, 4> > clippingTets(15);
DynamicList<FixedList<point, 4> > subjectTets(15);
label ntOld = 0, ntNew = 0;
vector oldCentre = vector::zero, newCentre = vector::zero;
// Configure tets from oldCell
forAll(oldCell, faceI)
{
const face& oldFace = mesh_.faces()[oldCell[faceI]];
vector fCentre = oldFace.centre(oldPoints);
if (oldFace.size() == 3)
{
// Add a new entry
subjectTets.append(FixedList<point, 4>(vector::zero));
subjectTets[ntOld][0] = oldPoints[oldFace[0]];
subjectTets[ntOld][1] = oldPoints[oldFace[1]];
subjectTets[ntOld][2] = oldPoints[oldFace[2]];
ntOld++;
}
else
{
forAll(oldFace, pI)
{
// Add a new entry
subjectTets.append(FixedList<point, 4>(vector::zero));
subjectTets[ntOld][0] = oldPoints[oldFace[pI]];
subjectTets[ntOld][1] = oldPoints[oldFace.nextLabel(pI)];
subjectTets[ntOld][2] = fCentre;
ntOld++;
}
}
oldCentre += fCentre;
}
// Configure tets from newCell
forAll(newCell, faceI)
{
const face& newFace = newFaces_[newCell[faceI]];
vector fCentre = newFace.centre(newPoints);
if (newFace.size() == 3)
{
// Add a new entry
clippingTets.append(FixedList<point, 4>(vector::zero));
clippingTets[ntNew][0] = newPoints[newFace[0]];
clippingTets[ntNew][1] = newPoints[newFace[1]];
clippingTets[ntNew][2] = newPoints[newFace[2]];
ntNew++;
}
else
{
forAll(newFace, pI)
{
// Add a new entry
clippingTets.append(FixedList<point, 4>(vector::zero));
clippingTets[ntNew][0] = newPoints[newFace[pI]];
clippingTets[ntNew][1] = newPoints[newFace.nextLabel(pI)];
clippingTets[ntNew][2] = fCentre;
ntNew++;
}
}
newCentre += fCentre;
}
oldCentre /= oldCell.size();
newCentre /= newCell.size();
// Fill-in last points for all tets
forAll(subjectTets, i)
{
subjectTets[i][3] = oldCentre;
}
forAll(clippingTets, i)
{
clippingTets[i][3] = newCentre;
}
// Accumulate volume / centroid over all intersections
bool foundIntersect = false;
scalar totalVolume = 0.0;
vector totalCentre = vector::zero;
// Loop through all clipping tets
forAll(clippingTets, i)
{
// Initialize the intersector
tetIntersection intersector(clippingTets[i]);
// Test for intersection and evaluate
// against all subject tets
forAll(subjectTets, j)
{
intersects = intersector.evaluate(subjectTets[j]);
if (intersects)
{
scalar volume = 0.0;
vector centre = vector::zero;
// Fetch volume and centre
intersector.getVolumeAndCentre(volume, centre);
// Accumulate volume / centroid
totalVolume += volume;
totalCentre += (volume * centre);
foundIntersect = true;
if (output)
{
writeVTK
(
"polyhedronIntersect_"
+ Foam::name(newIndex)
+ '_'
+ Foam::name(oldIndex)
+ '_' + Foam::name(i)
+ '_' + Foam::name(j),
intersector.getIntersection()
);
}
}
}
}
// Size-up the internal lists
if (foundIntersect && !output)
{
// Normalize centre
totalCentre /= totalVolume + VSMALL;
// Normalize and check if this is worth it
if (mag(totalVolume/normFactor_) > SMALL)
{
// Size-up the internal lists
meshOps::sizeUpList((oldIndex - offset), parents_);
meshOps::sizeUpList(totalVolume, weights_);
meshOps::sizeUpList(totalCentre, centres_);
intersects = true;
}
else
{
intersects = false;
}
}
}
else
{
// Configure points for clipping tetrahedron
FixedList<point, 4> clippingTet(vector::zero);
// Configure the clipping tetrahedron.
const face& firstNewFace = newFaces_[newCell[0]];
const face& secondNewFace = newFaces_[newCell[1]];
// Find the isolated point
label fourthNewPoint =
(
meshOps::findIsolatedPoint
(
firstNewFace,
secondNewFace
)
);
// Fill in points
clippingTet[0] = newPoints[firstNewFace[0]];
clippingTet[1] = newPoints[firstNewFace[1]];
clippingTet[2] = newPoints[firstNewFace[2]];
clippingTet[3] = newPoints[fourthNewPoint];
// Configure points for subject tetrahedron
FixedList<point, 4> subjectTet(vector::zero);
// Configure the subject tetrahedron.
const face& firstOldFace = mesh_.faces()[oldCell[0]];
const face& secondOldFace = mesh_.faces()[oldCell[1]];
// Find the isolated point
label fourthOldPoint =
(
meshOps::findIsolatedPoint
(
firstOldFace,
secondOldFace
)
);
// Fill in points
subjectTet[0] = oldPoints[firstOldFace[0]];
subjectTet[1] = oldPoints[firstOldFace[1]];
subjectTet[2] = oldPoints[firstOldFace[2]];
subjectTet[3] = oldPoints[fourthOldPoint];
// Initialize the intersector
tetIntersection intersector(clippingTet);
// Test for intersection and evaluate
intersects = intersector.evaluate(subjectTet);
if (intersects)
{
scalar volume = 0.0;
vector centre = vector::zero;
// Fetch volume and centre
intersector.getVolumeAndCentre(volume, centre);
// Normalize and check if this is worth it
if (mag(volume/normFactor_) > SMALL)
{
if (output)
{
writeVTK
(
"tetIntersectNew_"
+ Foam::name(newIndex),
List<FixedList<point, 4> >(1, clippingTet)
);
writeVTK
(
"tetIntersectOld_"
+ Foam::name(newIndex)
+ '_'
+ Foam::name(oldIndex),
List<FixedList<point, 4> >(1, subjectTet)
);
writeVTK
(
"tetIntersect_"
+ Foam::name(newIndex)
+ '_'
+ Foam::name(oldIndex),
intersector.getIntersection()
);
}
// Size-up the internal lists
meshOps::sizeUpList((oldIndex - offset), parents_);
meshOps::sizeUpList(volume, weights_);
meshOps::sizeUpList(centre, centres_);
}
else
{
intersects = false;
}
}
}
return intersects;
}
//- Write out tets as a VTK
void cellSetAlgorithm::writeVTK
(
const word& name,
const List<FixedList<point, 4> >& tetList
) const
{
// Fill up all points
label pI = 0;
List<cell> allTets(tetList.size(), cell(4));
pointField allPoints(4 * tetList.size());
forAll(tetList, tetI)
{
allTets[tetI][0] = pI;
allPoints[pI++] = tetList[tetI][0];
allTets[tetI][1] = pI;
allPoints[pI++] = tetList[tetI][1];
allTets[tetI][2] = pI;
allPoints[pI++] = tetList[tetI][2];
allTets[tetI][3] = pI;
allPoints[pI++] = tetList[tetI][3];
}
// Write out in cell-to-node addressing
meshOps::writeVTK
(
mesh_,
name,
identity(tetList.size()),
3,
allPoints,
List<edge>(0),
List<face>(0),
allTets,
List<label>(0)
);
}
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

128
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/cellSetAlgorithm.H Normal file
View file

@ -0,0 +1,128 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
cellSetAlgorithm
Description
Class for convexSetAlgorithms pertaining to cells
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
cellSetAlgorithm.C
\*---------------------------------------------------------------------------*/
#ifndef cellSetAlgorithm_H
#define cellSetAlgorithm_H
#include "convexSetAlgorithm.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class cellSetAlgorithm Declaration
\*---------------------------------------------------------------------------*/
class cellSetAlgorithm
:
public convexSetAlgorithm
{
// Private Member Functions
//- Disallow default bitwise copy construct
cellSetAlgorithm(const cellSetAlgorithm&);
//- Disallow default bitwise assignment
void operator=(const cellSetAlgorithm&);
//- Write out tets as a VTK
void writeVTK
(
const word&,
const List<FixedList<point, 4> >&
) const;
public:
//- Constructor
cellSetAlgorithm
(
const polyMesh& mesh,
const pointField& newPoints,
const UList<edge>& newEdges,
const UList<face>& newFaces,
const UList<cell>& newCells,
const UList<label>& newOwner,
const UList<label>& newNeighbour
);
//- Destructor
virtual ~cellSetAlgorithm()
{}
//- Virtual member functions
// Dimensions of the algorithm
virtual label dimension() const
{
return 3;
}
// Compute normFactor
virtual void computeNormFactor(const label index) const;
// Check whether the bounding box contains the entity
virtual bool contains(const label index) const;
// Compute intersections
virtual bool computeIntersection
(
const label newIndex,
const label oldIndex,
const label offset,
bool output
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

565
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/convexSetAlgorithm.C Normal file
View file

@ -0,0 +1,565 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
convexSetAlgorithm
Description
Base class for convexSetAlgorithms
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#include "Time.H"
#include "IOMap.H"
#include "meshOps.H"
#include "polyMesh.H"
#include "objectMap.H"
#include "edgeIOList.H"
#include "cellIOList.H"
#include "convexSetAlgorithm.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
convexSetAlgorithm::convexSetAlgorithm
(
const polyMesh& mesh,
const pointField& newPoints,
const UList<edge>& newEdges,
const UList<face>& newFaces,
const UList<cell>& newCells,
const UList<label>& newOwner,
const UList<label>& newNeighbour
)
:
nOldPoints_(mesh.nPoints()),
mesh_(mesh),
newPoints_(newPoints),
newEdges_(newEdges),
newFaces_(newFaces),
newCells_(newCells),
newOwner_(newOwner),
newNeighbour_(newNeighbour)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Obtain map weighting factors
void convexSetAlgorithm::computeWeights
(
const label index,
const label offset,
const labelList& mapCandidates,
const labelListList& oldNeighbourList,
labelList& parents,
scalarField& weights,
vectorField& centres,
bool output
)
{
if (parents.size() || weights.size() || centres.size())
{
FatalErrorIn
(
"\n\n"
"void convexSetAlgorithm::computeWeights\n"
"(\n"
" const label index,\n"
" const label offset,\n"
" const labelList& mapCandidates,\n"
" const labelListList& oldNeighbourList,\n"
" labelList& parents,\n"
" scalarField& weights,\n"
" vectorField& centres,\n"
" bool output\n"
")\n"
)
<< " Addressing has already been calculated." << nl
<< " Index: " << index << nl
<< " Offset: " << offset << nl
<< " Type: " << (dimension() == 2 ? "Face" : "Cell") << nl
<< " mapCandidates: " << mapCandidates << nl
<< " Parents: " << parents << nl
<< " Weights: " << weights << nl
<< " Centres: " << centres << nl
<< abort(FatalError);
}
// Do nothing for empty lists
if (mapCandidates.empty() || oldNeighbourList.empty())
{
return;
}
bool changed;
label nAttempts = 0, nIntersects = 0;
// Calculate the algorithm normFactor
computeNormFactor(index);
// Maintain a check-list
labelHashSet checked, skipped;
// Loop and add intersections until nothing changes
do
{
// Reset flag
changed = false;
// Fetch the set of candidates
labelList checkList;
if (nAttempts == 0)
{
checkList = mapCandidates;
}
else
{
checkList = checked.toc();
}
forAll(checkList, indexI)
{
labelList checkEntities;
if (nAttempts == 0)
{
checkEntities = labelList(1, checkList[indexI] - offset);
}
else
{
checkEntities = oldNeighbourList[checkList[indexI]];
}
forAll(checkEntities, entityI)
{
label checkEntity = checkEntities[entityI];
// Skip if this is already
// on the checked / skipped list
if
(
(checked.found(checkEntity)) ||
(skipped.found(checkEntity))
)
{
continue;
}
bool intersect =
(
computeIntersection
(
index,
checkEntity + offset,
offset,
output
)
);
if (intersect)
{
nIntersects++;
if (!checked.found(checkEntity))
{
checked.insert(checkEntity);
}
changed = true;
}
else
{
// Add to the skipped list
if (!skipped.found(checkEntity))
{
skipped.insert(checkEntity);
}
}
}
}
if (nAttempts == 0 && !changed)
{
// Need to setup a rescue mechanism.
labelHashSet rescue;
forAll(mapCandidates, cI)
{
if (!rescue.found(mapCandidates[cI] - offset))
{
rescue.insert(mapCandidates[cI] - offset);
}
}
for (label level = 0; level < 10; level++)
{
labelList initList = rescue.toc();
forAll(initList, fI)
{
const labelList& ff = oldNeighbourList[initList[fI]];
forAll(ff, entityI)
{
if (!rescue.found(ff[entityI]))
{
rescue.insert(ff[entityI]);
}
}
}
}
labelList finalList = rescue.toc();
forAll(finalList, entityI)
{
label checkEntity = finalList[entityI];
bool intersect =
(
computeIntersection
(
index,
checkEntity + offset,
offset,
output
)
);
if (intersect)
{
nIntersects++;
if (!checked.found(checkEntity))
{
checked.insert(checkEntity);
}
changed = true;
break;
}
}
if (!changed)
{
// No point in continuing further...
break;
}
}
nAttempts++;
// Break out if we're taking too long
if (nAttempts > 20)
{
break;
}
} while (changed);
// Normalize weights
normalize(false);
// Populate lists
populateLists(parents, centres, weights);
}
// Output an entity as a VTK file
void convexSetAlgorithm::writeVTK
(
const word& name,
const label entity,
const label primitiveType,
const bool useOldConnectivity
) const
{
writeVTK
(
name,
labelList(1, entity),
primitiveType,
useOldConnectivity
);
}
// Output a list of entities as a VTK file
void convexSetAlgorithm::writeVTK
(
const word& name,
const labelList& cList,
const label primitiveType,
const bool useOldConnectivity
) const
{
if (useOldConnectivity)
{
const polyMesh& mesh = this->mesh_;
meshOps::writeVTK
(
mesh,
name,
cList,
primitiveType,
mesh.points(),
mesh.edges(),
mesh.faces(),
mesh.cells(),
mesh.faceOwner()
);
}
else
{
meshOps::writeVTK
(
this->mesh_,
name,
cList,
primitiveType,
newPoints_,
newEdges_,
newFaces_,
newCells_,
newOwner_
);
}
}
bool convexSetAlgorithm::consistent(const scalar tolerance) const
{
if (weights_.size())
{
scalar normError = mag(1.0 - (sum(weights_)/normFactor_));
if (normError < tolerance)
{
return true;
}
else
{
return false;
}
}
return false;
}
// Return the normFactor
scalar convexSetAlgorithm::normFactor() const
{
return normFactor_;
}
// Normalize stored weights
void convexSetAlgorithm::normalize(bool normSum) const
{
if (normSum)
{
if (weights_.size())
{
weights_ /= sum(weights_);
}
}
else
{
if (weights_.size())
{
weights_ /= normFactor_;
}
}
}
// Extract weights and centres to lists
void convexSetAlgorithm::populateLists
(
labelList& parents,
vectorField& centres,
scalarField& weights
) const
{
// Clear inputs
parents.clear();
centres.clear();
weights.clear();
if (weights_.size())
{
parents = parents_;
centres = centres_;
weights = weights_;
}
}
// Write out connectivity information to disk
bool convexSetAlgorithm::write() const
{
pointIOField
(
IOobject
(
"newPoints",
mesh_.time().timeName(),
"convexSetAlgorithm",
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
newPoints_
).writeObject
(
IOstream::BINARY,
IOstream::currentVersion,
IOstream::COMPRESSED
);
edgeIOList
(
IOobject
(
"newEdges",
mesh_.time().timeName(),
"convexSetAlgorithm",
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
newEdges_
).writeObject
(
IOstream::BINARY,
IOstream::currentVersion,
IOstream::COMPRESSED
);
faceIOList
(
IOobject
(
"newFaces",
mesh_.time().timeName(),
"convexSetAlgorithm",
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
newFaces_
).writeObject
(
IOstream::BINARY,
IOstream::currentVersion,
IOstream::COMPRESSED
);
cellIOList
(
IOobject
(
"newCells",
mesh_.time().timeName(),
"convexSetAlgorithm",
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
newCells_
).writeObject
(
IOstream::BINARY,
IOstream::currentVersion,
IOstream::COMPRESSED
);
labelIOList
(
IOobject
(
"newOwner",
mesh_.time().timeName(),
"convexSetAlgorithm",
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
newOwner_
).writeObject
(
IOstream::BINARY,
IOstream::currentVersion,
IOstream::COMPRESSED
);
labelIOList
(
IOobject
(
"newNeighbour",
mesh_.time().timeName(),
"convexSetAlgorithm",
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
newNeighbour_
).writeObject
(
IOstream::BINARY,
IOstream::currentVersion,
IOstream::COMPRESSED
);
return true;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

198
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/convexSetAlgorithm.H Normal file
View file

@ -0,0 +1,198 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
convexSetAlgorithm
Description
Base class for convexSetAlgorithms
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
convexSetAlgorithm.C
\*---------------------------------------------------------------------------*/
#ifndef convexSetAlgorithm_H
#define convexSetAlgorithm_H
#include "Map.H"
#include "label.H"
#include "edgeList.H"
#include "faceList.H"
#include "cellList.H"
#include "boundBox.H"
#include "objectMap.H"
#include "vectorField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class polyMesh;
/*---------------------------------------------------------------------------*\
Class convexSetAlgorithm Declaration
\*---------------------------------------------------------------------------*/
class convexSetAlgorithm
{
protected:
// Protected data
const label nOldPoints_;
//- References to old-level connectivity
// [Before topo-changes, at old point-positions]
const polyMesh& mesh_;
//- References to new-level connectivity
// [After topo-changes, at old point-positions]
const pointField& newPoints_;
const UList<edge>& newEdges_;
const UList<face>& newFaces_;
const UList<cell>& newCells_;
const UList<label>& newOwner_;
const UList<label>& newNeighbour_;
//- Entity parents
mutable labelList parents_;
//- Internal data members
mutable boundBox box_;
mutable vector refNorm_;
mutable scalar normFactor_;
mutable vectorField centres_;
mutable scalarField weights_;
public:
//- Constructor
// Construct from components
convexSetAlgorithm
(
const polyMesh& mesh,
const pointField& newPoints,
const UList<edge>& newEdges,
const UList<face>& newFaces,
const UList<cell>& newCells,
const UList<label>& newOwner,
const UList<label>& newNeighbour
);
//- Destructor
virtual ~convexSetAlgorithm()
{}
//- Member functions
// Dimensions of the algorithm
virtual label dimension() const = 0;
// Return true if accumulated weights are consistent
virtual bool consistent(const scalar tolerance) const;
// Return the normFactor
virtual scalar normFactor() const;
// Normalize stored weights
virtual void normalize(bool normSum) const;
// Check whether the bounding box contains the entity
virtual bool contains(const label index) const = 0;
// Extract weights and centres to lists
virtual void populateLists
(
labelList& parents,
vectorField& centres,
scalarField& weights
) const;
// Compute normFactor
virtual void computeNormFactor(const label index) const = 0;
// Compute intersections
virtual bool computeIntersection
(
const label newIndex,
const label oldIndex,
const label offset,
bool output
) const = 0;
// Obtain map weighting factors
virtual void computeWeights
(
const label index,
const label offset,
const labelList& mapCandidates,
const labelListList& oldNeighbourList,
labelList& parents,
scalarField& weights,
vectorField& centres,
bool output = false
);
// Write out connectivity information to disk
bool write() const;
// Output an entity as a VTK file
void writeVTK
(
const word& name,
const label entity,
const label primitiveType = 3,
const bool useOldConnectivity = false
) const;
// Output a list of entities as a VTK file
void writeVTK
(
const word& name,
const labelList& cList,
const label primitiveType = 3,
const bool useOldConnectivity = false
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

415
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/faceSetAlgorithm.C Normal file
View file

@ -0,0 +1,415 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
faceSetAlgorithm
Description
Class for convexSetAlgorithms pertaining to faces
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#include "faceSetAlgorithm.H"
#include "meshOps.H"
#include "triFace.H"
#include "polyMesh.H"
#include "triIntersection.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
faceSetAlgorithm::faceSetAlgorithm
(
const polyMesh& mesh,
const pointField& newPoints,
const UList<edge>& newEdges,
const UList<face>& newFaces,
const UList<cell>& newCells,
const UList<label>& newOwner,
const UList<label>& newNeighbour
)
:
convexSetAlgorithm
(
mesh,
newPoints,
newEdges,
newFaces,
newCells,
newOwner,
newNeighbour
)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void faceSetAlgorithm::computeNormFactor(const label index) const
{
// Clear existing fields
parents_.clear();
centres_.clear();
weights_.clear();
// Compute refNorm and normFactor
refNorm_ = newFaces_[index].normal(newPoints_);
normFactor_ = mag(refNorm_);
// Normalize for later use
refNorm_ /= normFactor_ + VSMALL;
// Compute a bounding box around the face
box_ = boundBox(newFaces_[index].points(newPoints_), false);
vector minToXb = (box_.min() - box_.midpoint());
vector maxToXb = (box_.max() - box_.midpoint());
// Scale it by a bit
box_.min() += (1.5 * minToXb);
box_.max() += (1.5 * maxToXb);
}
// Check whether the bounding box contains the entity
bool faceSetAlgorithm::contains(const label index) const
{
// Fetch old face
const face& checkFace = mesh_.faces()[index];
// Check if the bounding box contains any of the supplied points
forAll(checkFace, pointI)
{
if (box_.contains(mesh_.points()[checkFace[pointI]]))
{
return true;
}
}
return false;
}
// Compute intersection
bool faceSetAlgorithm::computeIntersection
(
const label newIndex,
const label oldIndex,
const label offset,
bool output
) const
{
bool intersects = false;
const pointField& newPoints = newPoints_;
const pointField& oldPoints = mesh_.points();
const face& newFace = newFaces_[newIndex];
const face& oldFace = mesh_.faces()[oldIndex];
// Compute the normal for the old face
vector oldNorm = oldFace.normal(oldPoints);
// Normalize
oldNorm /= mag(oldNorm) + VSMALL;
if ((oldNorm & refNorm_) < 0.0)
{
// Opposite face orientation. Skip it.
return false;
}
// Check if decomposition is necessary
if (oldFace.size() > 3 || newFace.size() > 3)
{
// Decompose new / old faces
DynamicList<FixedList<point, 3> > clippingTris(15);
DynamicList<FixedList<point, 3> > subjectTris(15);
label ntOld = 0, ntNew = 0;
vector oldCentre = vector::zero, newCentre = vector::zero;
if (oldFace.size() == 3)
{
// Add a new entry
subjectTris.append(FixedList<point, 3>(vector::zero));
subjectTris[ntOld][0] = oldPoints[oldFace[0]];
subjectTris[ntOld][1] = oldPoints[oldFace[1]];
subjectTris[ntOld][2] = oldPoints[oldFace[2]];
ntOld++;
}
else
{
// Configure tris from oldFace
vector ofCentre = oldFace.centre(oldPoints);
forAll(oldFace, pI)
{
// Add a new entry
subjectTris.append(FixedList<point, 3>(vector::zero));
subjectTris[ntOld][0] = oldPoints[oldFace[pI]];
subjectTris[ntOld][1] = oldPoints[oldFace.nextLabel(pI)];
subjectTris[ntOld][2] = ofCentre;
ntOld++;
}
}
if (newFace.size() == 3)
{
// Add a new entry
clippingTris.append(FixedList<point, 3>(vector::zero));
clippingTris[ntNew][0] = newPoints[newFace[0]];
clippingTris[ntNew][1] = newPoints[newFace[1]];
clippingTris[ntNew][2] = newPoints[newFace[2]];
ntNew++;
}
else
{
// Configure tris from newFace
vector nfCentre = newFace.centre(newPoints);
forAll(newFace, pI)
{
// Add a new entry
clippingTris.append(FixedList<point, 3>(vector::zero));
clippingTris[ntNew][0] = newPoints[newFace[pI]];
clippingTris[ntNew][1] = newPoints[newFace.nextLabel(pI)];
clippingTris[ntNew][2] = nfCentre;
ntNew++;
}
}
// Accumulate area / centroid over all intersections
bool foundIntersect = false;
scalar totalArea = 0.0;
vector totalCentre = vector::zero;
// Loop through all clipping tris
forAll(clippingTris, i)
{
// Initialize the intersector
triIntersection intersector(clippingTris[i]);
// Test for intersection and evaluate
// against all subject tris
forAll(subjectTris, j)
{
intersects = intersector.evaluate(subjectTris[j]);
if (intersects)
{
scalar area = 0.0;
vector centre = vector::zero;
// Fetch area and centre
intersector.getAreaAndCentre(area, centre);
// Accumulate area / centroid
totalArea += area;
totalCentre += (area * centre);
foundIntersect = true;
if (output)
{
writeVTK
(
"polygonIntersect_"
+ Foam::name(newIndex)
+ '_'
+ Foam::name(oldIndex)
+ '_' + Foam::name(i)
+ '_' + Foam::name(j),
intersector.getIntersection()
);
}
}
}
}
// Size-up the internal lists
if (foundIntersect && !output)
{
// Normalize centre
totalCentre /= totalArea + VSMALL;
// Normalize and check if this is worth it
if (mag(totalArea/normFactor_) > SMALL)
{
// Size-up the internal lists
meshOps::sizeUpList((oldIndex - offset), parents_);
meshOps::sizeUpList(totalArea, weights_);
meshOps::sizeUpList(totalCentre, centres_);
intersects = true;
}
else
{
intersects = false;
}
}
}
else
{
// Configure points for clipping triangle
FixedList<point, 3> clippingTri(vector::zero);
// Fill in points
clippingTri[0] = newPoints[newFace[0]];
clippingTri[1] = newPoints[newFace[1]];
clippingTri[2] = newPoints[newFace[2]];
// Configure points for subject triangle
FixedList<point, 3> subjectTri(vector::zero);
// Fill in points
subjectTri[0] = oldPoints[oldFace[0]];
subjectTri[1] = oldPoints[oldFace[1]];
subjectTri[2] = oldPoints[oldFace[2]];
// Initialize the intersector
triIntersection intersector(clippingTri);
// Test for intersection and evaluate
intersects = intersector.evaluate(subjectTri);
if (intersects)
{
scalar area;
vector centre;
// Fetch area and centre
intersector.getAreaAndCentre(area, centre);
// Normalize and check if this is worth it
if (mag(area/normFactor_) > SMALL)
{
if (output)
{
writeVTK
(
"triIntersectNew_"
+ Foam::name(newIndex),
List<FixedList<point, 3> >(1, clippingTri)
);
writeVTK
(
"triIntersectOld_"
+ Foam::name(newIndex)
+ '_'
+ Foam::name(oldIndex),
List<FixedList<point, 3> >(1, subjectTri)
);
writeVTK
(
"triIntersect_"
+ Foam::name(newIndex)
+ '_'
+ Foam::name(oldIndex),
intersector.getIntersection()
);
}
// Size-up the internal lists
meshOps::sizeUpList((oldIndex - offset), parents_);
meshOps::sizeUpList(area, weights_);
meshOps::sizeUpList(centre, centres_);
}
else
{
intersects = false;
}
}
}
return intersects;
}
//- Write out tris as a VTK
void faceSetAlgorithm::writeVTK
(
const word& name,
const List<FixedList<point, 3> >& triList
) const
{
// Fill up all points
label pI = 0;
List<face> allTris(triList.size(), face(3));
pointField allPoints(3 * triList.size());
forAll(triList, triI)
{
allTris[triI][0] = pI;
allPoints[pI++] = triList[triI][0];
allTris[triI][1] = pI;
allPoints[pI++] = triList[triI][1];
allTris[triI][2] = pI;
allPoints[pI++] = triList[triI][2];
}
// Write out in face-to-node addressing
meshOps::writeVTK
(
mesh_,
name,
identity(triList.size()),
2,
allPoints,
List<edge>(0),
allTris,
List<cell>(0),
List<label>(0)
);
}
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

128
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/faceSetAlgorithm.H Normal file
View file

@ -0,0 +1,128 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
faceSetAlgorithm
Description
Class for convexSetAlgorithms pertaining to faces
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
faceSetAlgorithm.C
\*---------------------------------------------------------------------------*/
#ifndef faceSetAlgorithm_H
#define faceSetAlgorithm_H
#include "convexSetAlgorithm.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class faceSetAlgorithm Declaration
\*---------------------------------------------------------------------------*/
class faceSetAlgorithm
:
public convexSetAlgorithm
{
// Private Member Functions
//- Disallow default bitwise copy construct
faceSetAlgorithm(const faceSetAlgorithm&);
//- Disallow default bitwise assignment
void operator=(const faceSetAlgorithm&);
//- Write out tris as a VTK
void writeVTK
(
const word&,
const List<FixedList<point, 3> >&
) const;
public:
//- Constructor
faceSetAlgorithm
(
const polyMesh& mesh,
const pointField& newPoints,
const UList<edge>& newEdges,
const UList<face>& newFaces,
const UList<cell>& newCells,
const UList<label>& newOwner,
const UList<label>& newNeighbour
);
//- Destructor
virtual ~faceSetAlgorithm()
{}
//- Virtual member functions
// Dimensions of the algorithm
virtual label dimension() const
{
return 2;
}
// Compute normFactor
virtual void computeNormFactor(const label index) const;
// Check whether the bounding box contains the entity
virtual bool contains(const label index) const;
// Compute intersections
virtual bool computeIntersection
(
const label newIndex,
const label oldIndex,
const label offset,
bool output
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

145
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/tetIntersection.H Normal file
View file

@ -0,0 +1,145 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
tetIntersection
Description
Implementation of the tetrahedron intersection algorithm given in:
D.H. Eberly, 3D Game Engine Design: A Practical Approach to Real-time
Computer Graphics, Morgan Kaufmann, 2001.
Geometric Tools, LLC
Distributed under the Boost Software License, Version 1.0.
http://www.boost.org/LICENSE_1_0.txt
Implemented by
Sandeep Menon
University of Massachusetts Amherst
SourceFiles
tetIntersectionI.H
\*---------------------------------------------------------------------------*/
#ifndef tetIntersection_H
#define tetIntersection_H
#include "point.H"
#include "label.H"
#include "Tuple2.H"
#include "FixedList.H"
#include "DynamicList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class tetIntersection Declaration
\*---------------------------------------------------------------------------*/
class tetIntersection
{
// Private data
//- Const reference to clipping tetrahedron
const FixedList<point, 4>& clipTet_;
//- Hessian-normal plane definition
typedef Tuple2<vector, scalar> hPlane;
FixedList<hPlane, 4> clipPlanes_;
//- Magnitude of clipping tetrahedron
scalar clipTetMag_;
//- Tetrahedra used as temporaries
DynamicList<FixedList<point, 4> > inside_;
//- All intersection tets
DynamicList<FixedList<point, 4> > allTets_;
// Private Member Functions
//- Disallow default bitwise copy construct
tetIntersection(const tetIntersection&);
//- Disallow default bitwise assignment
void operator=(const tetIntersection&);
//- Compute clip-planes
inline void computeClipPlanes();
//- Split and decompose
inline void splitAndDecompose
(
const label tetPlaneIndex,
FixedList<point, 4>& tetra,
DynamicList<FixedList<point, 4> >& decompTets
) const;
public:
// Constructors
//- Construct from components
inline tetIntersection(const FixedList<point, 4>& clipTet);
// Destructor
inline ~tetIntersection();
// Member Functions
//- Return magnitude of clipping tetrahedron
inline scalar clipTetMag() const;
//- Evaluate for intersections against input tetrahedron
inline bool evaluate(const FixedList<point, 4>& subjectTet);
//- Return intersections
inline const DynamicList<FixedList<point, 4> >& getIntersection() const;
//- Evaluate and return volume / centroid
inline void getVolumeAndCentre(scalar& volume, vector& centre) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "tetIntersectionI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

387
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/tetIntersectionI.H Normal file
View file

@ -0,0 +1,387 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Implemented by
Sandeep Menon
University of Massachusetts Amherst
\*---------------------------------------------------------------------------*/
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Compute clip-planes
inline void tetIntersection::computeClipPlanes()
{
// Define edge vectors
vector edge10 = clipTet_[1] - clipTet_[0];
vector edge20 = clipTet_[2] - clipTet_[0];
vector edge30 = clipTet_[3] - clipTet_[0];
vector edge21 = clipTet_[2] - clipTet_[1];
vector edge31 = clipTet_[3] - clipTet_[1];
// Cross-products
clipPlanes_[0].first() = (edge20 ^ edge10);
clipPlanes_[1].first() = (edge10 ^ edge30);
clipPlanes_[2].first() = (edge30 ^ edge20);
clipPlanes_[3].first() = (edge21 ^ edge31);
// Normalize
clipPlanes_[0].first() /= mag(clipPlanes_[0].first()) + VSMALL;
clipPlanes_[1].first() /= mag(clipPlanes_[1].first()) + VSMALL;
clipPlanes_[2].first() /= mag(clipPlanes_[2].first()) + VSMALL;
clipPlanes_[3].first() /= mag(clipPlanes_[3].first()) + VSMALL;
// Compute magnitude of clipping tetrahedron
clipTetMag_ = (1.0 / 6.0) * (edge10 & clipPlanes_[3].first());
if (clipTetMag_ < 0.0)
{
// Reverse normal directions
clipPlanes_[0].first() = -clipPlanes_[0].first();
clipPlanes_[1].first() = -clipPlanes_[1].first();
clipPlanes_[2].first() = -clipPlanes_[2].first();
clipPlanes_[3].first() = -clipPlanes_[3].first();
// Reverse sign
clipTetMag_ = mag(clipTetMag_);
}
// Determine plane constants
clipPlanes_[0].second() = (clipTet_[0] & clipPlanes_[0].first());
clipPlanes_[1].second() = (clipTet_[1] & clipPlanes_[1].first());
clipPlanes_[2].second() = (clipTet_[2] & clipPlanes_[2].first());
clipPlanes_[3].second() = (clipTet_[3] & clipPlanes_[3].first());
}
// Split and decompose
inline void tetIntersection::splitAndDecompose
(
const label tetPlaneIndex,
FixedList<point, 4>& tetra,
DynamicList<FixedList<point, 4> >& decompTets
) const
{
FixedList<scalar, 4> C;
FixedList<vector, 4> tmpTetra;
FixedList<label, 4> pos, neg, zero;
label i = 0, nPos = 0, nNeg = 0, nZero = 0;
// Fetch reference to plane
const hPlane& tetPlane = clipPlanes_[tetPlaneIndex];
for (i = 0; i < 4; ++i)
{
// Compute distance to plane
C[i] = (tetra[i] & tetPlane.first()) - tetPlane.second();
if (C[i] > 0.0)
{
pos[nPos++] = i;
}
else
if (C[i] < 0.0)
{
neg[nNeg++] = i;
}
else
{
zero[nZero++] = i;
}
}
if (nNeg == 0)
{
return;
}
if (nPos == 0)
{
decompTets.append(tetra);
return;
}
// Tetrahedron is split by plane. Determine how it is split and how to
// decompose the negative-side portion into tetrahedra (6 cases).
scalar w0, w1, invCDiff;
vector intp[4];
if (nPos == 3)
{
// +++-
for (i = 0; i < nPos; ++i)
{
invCDiff = (1.0 / (C[pos[i]] - C[neg[0]]));
w0 = -C[neg[0]] * invCDiff;
w1 = +C[pos[i]] * invCDiff;
tetra[pos[i]] = (w0 * tetra[pos[i]]) + (w1 * tetra[neg[0]]);
}
decompTets.append(tetra);
}
else
if (nPos == 2)
{
if (nNeg == 2)
{
// ++--
for (i = 0; i < nPos; ++i)
{
invCDiff = (1.0 / (C[pos[i]] - C[neg[0]]));
w0 = -C[neg[0]] * invCDiff;
w1 = +C[pos[i]] * invCDiff;
intp[i] = (w0 * tetra[pos[i]]) + (w1 * tetra[neg[0]]);
}
for (i = 0; i < nNeg; ++i)
{
invCDiff = (1.0 / (C[pos[i]] - C[neg[1]]));
w0 = -C[neg[1]] * invCDiff;
w1 = +C[pos[i]] * invCDiff;
intp[i+2] = (w0 * tetra[pos[i]]) + (w1 * tetra[neg[1]]);
}
tetra[pos[0]] = intp[2];
tetra[pos[1]] = intp[1];
decompTets.append(tetra);
tmpTetra[0] = tetra[neg[1]];
tmpTetra[1] = intp[3];
tmpTetra[2] = intp[2];
tmpTetra[3] = intp[1];
decompTets.append(tmpTetra);
tmpTetra[0] = tetra[neg[0]];
tmpTetra[1] = intp[0];
tmpTetra[2] = intp[1];
tmpTetra[3] = intp[2];
decompTets.append(tmpTetra);
}
else
{
// ++-0
for (i = 0; i < nPos; ++i)
{
invCDiff = (1.0 / (C[pos[i]] - C[neg[0]]));
w0 = -C[neg[0]] * invCDiff;
w1 = +C[pos[i]] * invCDiff;
tetra[pos[i]] = (w0 * tetra[pos[i]]) + (w1 * tetra[neg[0]]);
}
decompTets.append(tetra);
}
}
else
if (nPos == 1)
{
if (nNeg == 3)
{
// +---
for (i = 0; i < nNeg; ++i)
{
invCDiff = (1.0 / (C[pos[0]] - C[neg[i]]));
w0 = -C[neg[i]] * invCDiff;
w1 = +C[pos[0]] * invCDiff;
intp[i] = (w0 * tetra[pos[0]]) + (w1 * tetra[neg[i]]);
}
tetra[pos[0]] = intp[0];
decompTets.append(tetra);
tmpTetra[0] = intp[0];
tmpTetra[1] = tetra[neg[1]];
tmpTetra[2] = tetra[neg[2]];
tmpTetra[3] = intp[1];
decompTets.append(tmpTetra);
tmpTetra[0] = tetra[neg[2]];
tmpTetra[1] = intp[1];
tmpTetra[2] = intp[2];
tmpTetra[3] = intp[0];
decompTets.append(tmpTetra);
}
else
if (nNeg == 2)
{
// +--0
for (i = 0; i < nNeg; ++i)
{
invCDiff = (1.0 / (C[pos[0]] - C[neg[i]]));
w0 = -C[neg[i]] * invCDiff;
w1 = +C[pos[0]] * invCDiff;
intp[i] = (w0 * tetra[pos[0]]) + (w1 * tetra[neg[i]]);
}
tetra[pos[0]] = intp[0];
decompTets.append(tetra);
tmpTetra[0] = intp[1];
tmpTetra[1] = tetra[zero[0]];
tmpTetra[2] = tetra[neg[1]];
tmpTetra[3] = intp[0];
decompTets.append(tmpTetra);
}
else
{
// +-00
invCDiff = (1.0 / (C[pos[0]] - C[neg[0]]));
w0 = -C[neg[0]] * invCDiff;
w1 = +C[pos[0]] * invCDiff;
tetra[pos[0]] = (w0 * tetra[pos[0]]) + (w1 * tetra[neg[0]]);
decompTets.append(tetra);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
inline tetIntersection::tetIntersection(const FixedList<point, 4>& clipTet)
:
clipTet_(clipTet),
clipTetMag_(0.0),
inside_(10),
allTets_(10)
{
// Pre-compute clipping planes
computeClipPlanes();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
inline tetIntersection::~tetIntersection()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Return magnitude of clipping tetrahedron
inline scalar tetIntersection::clipTetMag() const
{
return clipTetMag_;
}
// Evaluate for intersections
inline bool tetIntersection::evaluate(const FixedList<point, 4>& subjectTet)
{
// Clear lists
inside_.clear();
allTets_.clear();
// Add initial tetrahedron to list
allTets_.append(subjectTet);
// Clip second tetrahedron against planes of clipping tetrahedron
for (label i = 0; i < 4; i++)
{
forAll(allTets_, tetI)
{
splitAndDecompose(i, allTets_[tetI], inside_);
}
// Prep for next clipping plane
allTets_ = inside_;
inside_.clear();
}
return (allTets_.size() > 0);
}
// Return intersections
inline const DynamicList<FixedList<point, 4> >&
tetIntersection::getIntersection() const
{
return allTets_;
}
//- Evaluate and return volume / centroid
inline void tetIntersection::getVolumeAndCentre
(
scalar& volume,
vector& centre
) const
{
volume = 0.0;
centre = vector::zero;
forAll(allTets_, tetI)
{
const FixedList<point, 4>& t = allTets_[tetI];
// Calculate volume (no check for orientation)
scalar tV =
(
Foam::mag
(
(1.0/6.0) *
(
((t[1] - t[0]) ^ (t[2] - t[0])) & (t[3] - t[0])
)
)
);
// Calculate centroid
vector tC = (0.25 * (t[0] + t[1] + t[2] + t[3]));
volume += tV;
centre += (tV * tC);
}
centre /= volume + VSMALL;
}
}
// ************************************************************************* //

134
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/triIntersection.H Normal file
View file

@ -0,0 +1,134 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
triIntersection
Description
Implementation of the triangle intersection algorithm
Implemented by
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
triIntersectionI.H
\*---------------------------------------------------------------------------*/
#ifndef triIntersection_H
#define triIntersection_H
#include "point.H"
#include "label.H"
#include "Tuple2.H"
#include "FixedList.H"
#include "DynamicList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class triIntersection Declaration
\*---------------------------------------------------------------------------*/
class triIntersection
{
// Private data
//- Const reference to clipping triangle
const FixedList<point, 3>& clipTri_;
//- Clip triangle normal
vector tNorm_;
//- Hessian-normal plane definition
typedef Tuple2<vector, scalar> hPlane;
FixedList<hPlane, 3> clipPlanes_;
//- Triangles used as temporaries
DynamicList<FixedList<point, 3> > inside_;
//- All intersection tris
DynamicList<FixedList<point, 3> > allTris_;
// Private Member Functions
//- Disallow default bitwise copy construct
triIntersection(const triIntersection&);
//- Disallow default bitwise assignment
void operator=(const triIntersection&);
//- Compute clip-planes
inline void computeClipPlanes();
//- Split and decompose
inline void splitAndDecompose
(
const label triPlaneIndex,
FixedList<point, 3>& tri,
DynamicList<FixedList<point, 3> >& decompTris
) const;
public:
// Constructors
//- Construct from components
inline triIntersection(const FixedList<point, 3>& clipTri);
// Destructor
inline ~triIntersection();
// Member Functions
//- Evaluate for intersections against input triangle
inline bool evaluate(const FixedList<point, 3>& subjectTri);
//- Return intersections
inline const DynamicList<FixedList<point, 3> >& getIntersection() const;
//- Evaluate and return area / centroid
inline void getAreaAndCentre(scalar& area, vector& centre) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "triIntersectionI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

269
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/triIntersectionI.H Normal file
View file

@ -0,0 +1,269 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Implemented by
Sandeep Menon
University of Massachusetts Amherst
\*---------------------------------------------------------------------------*/
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Compute clip-planes
inline void triIntersection::computeClipPlanes()
{
// Face normal
tNorm_ = (clipTri_[1] - clipTri_[0]) ^ (clipTri_[2] - clipTri_[0]);
// Normalize
tNorm_ /= mag(tNorm_) + VSMALL;
// Edge normals
clipPlanes_[0].first() = ((clipTri_[1] - clipTri_[0]) ^ tNorm_);
clipPlanes_[1].first() = ((clipTri_[2] - clipTri_[1]) ^ tNorm_);
clipPlanes_[2].first() = ((clipTri_[0] - clipTri_[2]) ^ tNorm_);
// Normalize
clipPlanes_[0].first() /= mag(clipPlanes_[0].first()) + VSMALL;
clipPlanes_[1].first() /= mag(clipPlanes_[1].first()) + VSMALL;
clipPlanes_[2].first() /= mag(clipPlanes_[2].first()) + VSMALL;
// Determine plane constants
clipPlanes_[0].second() = (clipTri_[0] & clipPlanes_[0].first());
clipPlanes_[1].second() = (clipTri_[1] & clipPlanes_[1].first());
clipPlanes_[2].second() = (clipTri_[2] & clipPlanes_[2].first());
}
// Split and decompose
inline void triIntersection::splitAndDecompose
(
const label triPlaneIndex,
FixedList<point, 3>& tri,
DynamicList<FixedList<point, 3> >& decompTris
) const
{
FixedList<scalar, 3> C;
FixedList<vector, 3> tmpTri;
FixedList<label, 3> pos, neg, zero;
label i = 0, nPos = 0, nNeg = 0, nZero = 0;
// Fetch reference to plane
const hPlane& triPlane = clipPlanes_[triPlaneIndex];
for (i = 0; i < 3; ++i)
{
// Compute distance to plane
C[i] = (tri[i] & triPlane.first()) - triPlane.second();
if (C[i] > 0.0)
{
pos[nPos++] = i;
}
else
if (C[i] < 0.0)
{
neg[nNeg++] = i;
}
else
{
zero[nZero++] = i;
}
}
if (nNeg == 0)
{
return;
}
if (nPos == 0)
{
decompTris.append(tri);
return;
}
// Triangle is split by plane. Determine how it is split and how to
// decompose the negative-side portion into triangles (3 cases).
scalar w0, w1, invCDiff;
vector intp[3];
if (nPos == 2)
{
// ++-
for (i = 0; i < nPos; ++i)
{
invCDiff = (1.0 / (C[pos[i]] - C[neg[0]]));
w0 = -C[neg[0]] * invCDiff;
w1 = +C[pos[i]] * invCDiff;
tri[pos[i]] = (w0 * tri[pos[i]]) + (w1 * tri[neg[0]]);
}
decompTris.append(tri);
}
else
if (nPos == 1)
{
if (nNeg == 2)
{
// +--
for (i = 0; i < nNeg; ++i)
{
invCDiff = (1.0 / (C[pos[0]] - C[neg[i]]));
w0 = -C[neg[i]] * invCDiff;
w1 = +C[pos[0]] * invCDiff;
intp[i] = (w0 * tri[pos[0]]) + (w1 * tri[neg[i]]);
}
tri[pos[0]] = intp[0];
decompTris.append(tri);
tmpTri[0] = intp[0];
tmpTri[1] = tri[neg[1]];
tmpTri[2] = intp[1];
decompTris.append(tmpTri);
}
else
{
// +-0
invCDiff = (1.0 / (C[pos[0]] - C[neg[0]]));
w0 = -C[neg[0]] * invCDiff;
w1 = +C[pos[0]] * invCDiff;
tri[pos[0]] = (w0 * tri[pos[0]]) + (w1 * tri[neg[0]]);
decompTris.append(tri);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
inline triIntersection::triIntersection(const FixedList<point, 3>& clipTri)
:
clipTri_(clipTri),
inside_(10),
allTris_(10)
{
// Pre-compute clipping planes
computeClipPlanes();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
inline triIntersection::~triIntersection()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Evaluate for intersections
inline bool triIntersection::evaluate(const FixedList<point, 3>& subjectTri)
{
// Clear lists
inside_.clear();
allTris_.clear();
// Project subject triangle to clipping triangle plane
FixedList<point, 3> sTP(subjectTri);
for (label i = 0; i < 3; i++)
{
// Relative to zeroth face-point
sTP[i] -= clipTri_[0];
// Project to face
sTP[i] = clipTri_[0] + (sTP[i] - (sTP[i] & tNorm_)*tNorm_);
}
// Add initial triangle to list
allTris_.append(sTP);
// Clip second triangle against planes of clipping triangle
for (label i = 0; i < 3; i++)
{
forAll(allTris_, triI)
{
splitAndDecompose(i, allTris_[triI], inside_);
}
// Prep for next clipping plane
allTris_ = inside_;
inside_.clear();
}
return (allTris_.size() > 0);
}
// Return intersections
inline const DynamicList<FixedList<point, 3> >&
triIntersection::getIntersection() const
{
return allTris_;
}
//- Evaluate and return area / centroid
inline void triIntersection::getAreaAndCentre
(
scalar& area,
vector& centre
) const
{
area = 0.0;
centre = vector::zero;
forAll(allTris_, triI)
{
const FixedList<point, 3>& t = allTris_[triI];
// Calculate area (no check for orientation)
scalar tA = Foam::mag(0.5 * ((t[1] - t[0]) ^ (t[2] - t[0])));
// Calculate centroid
vector tC = (1.0 / 3.0) * (t[0] + t[1] + t[2]);
area += tA;
centre += (tA * tC);
}
centre /= area + VSMALL;
}
}
// ************************************************************************* //

636
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/coupleMap.C Normal file
View file

@ -0,0 +1,636 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
coupleMap
Description
Implementation of the coupleMap class
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#include "coupleMap.H"
#include "boolList.H"
#include "demandDrivenData.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const char* coupleMap::names[coupleMap::INVALID + 1] =
{
"BISECTION",
"COLLAPSE_FIRST",
"COLLAPSE_SECOND",
"COLLAPSE_MIDPOINT",
"REMOVE_CELL",
"MOVE_POINT",
"CONVERT_PATCH",
"INVALID"
};
// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
const char* coupleMap::asText(const opType oType)
{
return names[oType];
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
coupleMap::coupleMap
(
const IOobject& io,
const bool twoDMesh,
const bool isLocal,
const bool isSend,
const label patchIndex,
const label masterIndex,
const label slaveIndex
)
:
regIOobject(io),
twoDMesh_(twoDMesh),
isLocal_(isLocal),
isSend_(isSend),
patchIndex_(patchIndex),
masterIndex_(masterIndex),
slaveIndex_(slaveIndex),
nEntities_(-1),
edgesPtr_(NULL),
facesPtr_(NULL),
faceEdgesPtr_(NULL)
{
if
(
(io.readOpt() == IOobject::MUST_READ)
|| (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
)
{
// Construct an Istream and read from disk.
readData(readStream(typeName));
close();
}
}
// Construct as copy
coupleMap::coupleMap(const coupleMap& cm)
:
regIOobject(cm, true),
twoDMesh_(cm.twoDMesh_),
isLocal_(cm.isLocal_),
isSend_(cm.isSend_),
patchIndex_(cm.patchIndex_),
masterIndex_(cm.masterIndex_),
slaveIndex_(cm.slaveIndex_),
nEntities_(cm.nEntities_),
edgesPtr_(NULL),
facesPtr_(NULL),
faceEdgesPtr_(NULL)
{
if
(
(cm.readOpt() == IOobject::MUST_READ)
|| (cm.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
)
{
// Construct an Istream and read from disk.
readData(readStream(typeName));
close();
}
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
coupleMap::~coupleMap()
{
clearMaps();
clearBuffers();
clearAddressing();
}
// * * * * * * * * * * * * * * * Private Functions * * * * * * * * * * * * * //
void coupleMap::clearAddressing() const
{
deleteDemandDrivenData(edgesPtr_);
deleteDemandDrivenData(facesPtr_);
deleteDemandDrivenData(faceEdgesPtr_);
}
void coupleMap::makeEdges() const
{
// It is an error to attempt to recalculate
// if the pointer is already set
if (edgesPtr_)
{
FatalErrorIn("void coupleMap::makeEdges() const")
<< "Edges have already been calculated."
<< abort(FatalError);
}
label nEdges = nEntities(coupleMap::EDGE);
const labelList& eBuffer = entityBuffer(coupleMap::EDGE);
edgesPtr_ = new edgeList(nEdges);
edgeList& edges = *edgesPtr_;
forAll(edges, edgeI)
{
edges[edgeI][0] = eBuffer[(2*edgeI)+0];
edges[edgeI][1] = eBuffer[(2*edgeI)+1];
}
}
void coupleMap::makeFaces() const
{
// It is an error to attempt to recalculate
// if the pointer is already set
if (facesPtr_ || faceEdgesPtr_)
{
FatalErrorIn("void coupleMap::makeFaces() const")
<< "Faces have already been calculated."
<< abort(FatalError);
}
label nFaces = nEntities(coupleMap::FACE);
const labelList& fBuffer = entityBuffer(coupleMap::FACE);
const labelList& feBuffer = entityBuffer(coupleMap::FACE_EDGE);
const labelList& nfeBuffer = entityBuffer(coupleMap::NFE_BUFFER);
facesPtr_ = new faceList(nFaces);
faceEdgesPtr_ = new labelListList(nFaces);
faceList& faces = *facesPtr_;
labelListList& faceEdges = *faceEdgesPtr_;
label sumNFE = 0;
forAll(faces, faceI)
{
face& f = faces[faceI];
labelList& fe = faceEdges[faceI];
// Fetch the buffer value for 2D meshes
label nfe = twoDMesh_ ? nfeBuffer[faceI] : 3;
// Size up the lists
f.setSize(nfe, -1);
fe.setSize(nfe, -1);
for (label p = 0; p < nfe; p++)
{
f[p] = fBuffer[sumNFE + p];
fe[p] = feBuffer[sumNFE + p];
}
sumNFE += nfe;
}
if (sumNFE != nEntities(coupleMap::NFE_SIZE))
{
FatalErrorIn("void coupleMap::makeFaces() const")
<< " Mismatched buffer." << nl
<< " sumNFE: " << sumNFE << nl
<< " NFE_SIZE: " << nEntities(coupleMap::NFE_SIZE) << nl
<< abort(FatalError);
}
}
void coupleMap::makeFaceMap() const
{
// It is an error to attempt to recalculate
// if the map is already calculated
if (faceMap_.size())
{
FatalErrorIn("void coupleMap::makeFaceMap() const")
<< "faceMap has already been calculated."
<< abort(FatalError);
}
const Map<label>& mapFaces = entityMap(coupleMap::FACE);
// Size the list
faceMap_.setSize(mapFaces.size(), -1);
// Fill-in entries
forAllConstIter(Map<label>, mapFaces, fIter)
{
faceMap_[fIter.key()] = fIter();
}
}
void coupleMap::makeCellMap() const
{
// It is an error to attempt to recalculate
// if the map is already calculated
if (cellMap_.size())
{
FatalErrorIn("void coupleMap::makeCellMap() const")
<< "cellMap has already been calculated."
<< abort(FatalError);
}
const Map<label>& mapCells = entityMap(coupleMap::CELL);
// Size the list
cellMap_.setSize(mapCells.size(), -1);
// Fill-in entries
forAllConstIter(Map<label>, mapCells, cIter)
{
cellMap_[cIter.key()] = cIter();
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
pointField& coupleMap::pointBuffer() const
{
return pointBuffer_;
}
pointField& coupleMap::oldPointBuffer() const
{
return oldPointBuffer_;
}
labelList& coupleMap::subMeshPoints() const
{
return subMeshPoints_;
}
List<labelPair>& coupleMap::globalProcPoints() const
{
return globalProcPoints_;
}
void coupleMap::allocateBuffers() const
{
forAll(nEntities_, entityI)
{
if (nEntities_[entityI] < 0)
{
FatalErrorIn("void coupleMap::allocateBuffers() const")
<< " Entity sizes are not valid." << nl
<< " nEntities: " << nEntities_
<< abort(FatalError);
}
}
// Size up point buffers
pointBuffer().setSize(nEntities(POINT));
oldPointBuffer().setSize(nEntities(POINT));
// Size up connectivity buffers
entityBuffer(POINT).setSize
(
nEntities(GLOBAL_POINT)
- nEntities(SHARED_POINT)
);
// Set edge buffer
entityBuffer(EDGE).setSize(2 * nEntities(EDGE));
// Set addressing buffers
entityBuffer(OWNER).setSize(nEntities(FACE));
entityBuffer(NEIGHBOUR).setSize(nEntities(INTERNAL_FACE));
// Size up boundary buffers
entityBuffer(FACE_STARTS).setSize(nEntities(NBDY));
entityBuffer(FACE_SIZES).setSize(nEntities(NBDY));
entityBuffer(EDGE_STARTS).setSize(nEntities(NBDY));
entityBuffer(EDGE_SIZES).setSize(nEntities(NBDY));
entityBuffer(PATCH_ID).setSize(nEntities(NBDY));
// nFaceEdges buffer is required only for 2D,
// due to a mix of triangle / quad faces
if (twoDMesh_)
{
entityBuffer(NFE_BUFFER).setSize(nEntities(FACE));
}
// Allocate for variable size face-lists
entityBuffer(FACE).setSize(nEntities(NFE_SIZE));
entityBuffer(FACE_EDGE).setSize(nEntities(NFE_SIZE));
}
label coupleMap::findSlave
(
const label eType,
const label Index
) const
{
Map<label>::const_iterator it = entityMap_[eType].find(Index);
if (it == entityMap_[eType].end())
{
return -1;
}
else
{
return it();
}
}
label coupleMap::findMaster
(
const label eType,
const label Index
) const
{
Map<label>::const_iterator it = reverseEntityMap_[eType].find(Index);
if (it == reverseEntityMap_[eType].end())
{
return -1;
}
else
{
return it();
}
}
void coupleMap::removeSlave
(
const label eType,
const label Index
) const
{
Map<label>::iterator it = reverseEntityMap_[eType].find(Index);
if (it != reverseEntityMap_[eType].end())
{
reverseEntityMap_[eType].erase(it);
}
}
void coupleMap::removeMaster
(
const label eType,
const label Index
) const
{
Map<label>::iterator it = entityMap_[eType].find(Index);
if (it != entityMap_[eType].end())
{
entityMap_[eType].erase(it);
}
}
void coupleMap::mapSlave
(
const label eType,
const label master,
const label slave
) const
{
entityMap_[eType].set(master, slave);
}
void coupleMap::mapMaster
(
const label eType,
const label slave,
const label master
) const
{
reverseEntityMap_[eType].set(slave, master);
}
void coupleMap::pushOperation
(
const label index,
const opType oType,
const point& newPoint,
const point& oldPoint
) const
{
entityIndices_.setSize(entityIndices_.size() + 1, index);
entityOperations_.setSize(entityOperations_.size() + 1, oType);
if
(
oType == coupleMap::MOVE_POINT ||
oType == coupleMap::CONVERT_PATCH
)
{
moveNewPoints_.setSize(moveNewPoints_.size() + 1, newPoint);
moveOldPoints_.setSize(moveOldPoints_.size() + 1, oldPoint);
}
}
void coupleMap::transferMaps
(
const label eType,
Map<label>& newEntityMap,
Map<label>& newReverseEntityMap
) const
{
entityMap_[eType].transfer(newEntityMap);
reverseEntityMap_[eType].transfer(newReverseEntityMap);
}
void coupleMap::clearMaps() const
{
faceMap_.clear();
cellMap_.clear();
forAll(entityMap_, mapI)
{
entityMap_[mapI].clear();
reverseEntityMap_[mapI].clear();
}
}
void coupleMap::clearBuffers() const
{
pointBuffer_.clear();
oldPointBuffer_.clear();
subMeshPoints_.clear();
globalProcPoints_.clear();
forAll(entityBuffer_, bufferI)
{
entityBuffer_[bufferI].clear();
}
entityIndices_.clear();
entityOperations_.clear();
moveNewPoints_.clear();
moveOldPoints_.clear();
}
label coupleMap::nInternalFaces() const
{
return nEntities(INTERNAL_FACE);
}
const labelList& coupleMap::owner() const
{
return entityBuffer(OWNER);
}
const labelList& coupleMap::neighbour() const
{
return entityBuffer(NEIGHBOUR);
}
const edgeList& coupleMap::edges() const
{
if (!edgesPtr_)
{
makeEdges();
}
return *edgesPtr_;
}
const faceList& coupleMap::faces() const
{
if (!facesPtr_)
{
makeFaces();
}
return *facesPtr_;
}
const labelListList& coupleMap::faceEdges() const
{
if (!faceEdgesPtr_)
{
makeFaces();
}
return *faceEdgesPtr_;
}
const labelList& coupleMap::faceMap() const
{
if (faceMap_.empty())
{
makeFaceMap();
}
return faceMap_;
}
const labelList& coupleMap::cellMap() const
{
if (cellMap_.empty())
{
makeCellMap();
}
return cellMap_;
}
bool coupleMap::readData(Istream& is)
{
Map<label> tmpMap(is);
entityMap(coupleMap::POINT).transfer(tmpMap);
// Prepare the reversePointMap as well.
const Map<label>& pMap = entityMap(coupleMap::POINT);
Map<label>& rpMap = reverseEntityMap(coupleMap::POINT);
forAllConstIter(Map<label>, pMap, pIter)
{
rpMap.set(pIter(), pIter.key());
}
return !is.bad();
}
bool coupleMap::writeData(Ostream& os) const
{
// Only write-out point-map information
// to avoid geometric checking.
// The rest can be constructed topologically.
return (os << entityMap(coupleMap::POINT)).good();;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void coupleMap::operator=(const coupleMap& rhs)
{
// Check for assignment to self
if (this == &rhs)
{
FatalErrorIn("coupleMap::operator=(const coupleMap&)")
<< "Attempted assignment to self"
<< abort(FatalError);
}
}
} // End namespace Foam
// ************************************************************************* //

375
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/coupleMap.H Normal file
View file

@ -0,0 +1,375 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
coupleMap
Description
Coupled patch information with registry support.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
coupleMapI.H
coupleMap.C
\*---------------------------------------------------------------------------*/
#ifndef coupleMap_H
#define coupleMap_H
#include "regIOobject.H"
#include "labelList.H"
#include "pointField.H"
#include "fieldTypes.H"
#include "edgeList.H"
#include "faceList.H"
#include "cellList.H"
#include "labelPair.H"
#include "Map.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class coupleMap Declaration
\*---------------------------------------------------------------------------*/
class coupleMap
:
public regIOobject
{
public:
// Public data types
//- Enumerants for entity-types
enum entityType
{
// Buffers and sizes
POINT = 0,
EDGE = 1,
FACE = 2,
// Sizes only
CELL = 3,
INTERNAL_EDGE = 4,
INTERNAL_FACE = 5,
SHARED_POINT = 6,
GLOBAL_POINT = 7,
NFE_SIZE = 8,
NBDY = 9,
// Buffers only
OWNER = 3,
NEIGHBOUR = 4,
FACE_EDGE = 5,
NFE_BUFFER = 6,
FACE_STARTS = 7,
FACE_SIZES = 8,
EDGE_STARTS = 9,
EDGE_SIZES = 10,
PATCH_ID = 11
};
//- Enumerants for operations
enum opType
{
BISECTION = 0,
COLLAPSE_FIRST = 1,
COLLAPSE_SECOND = 2,
COLLAPSE_MIDPOINT = 3,
REMOVE_CELL = 4,
MOVE_POINT = 5,
CONVERT_PATCH = 6,
INVALID
};
private:
// Private data
mutable bool twoDMesh_;
// Flags for coupled entities
mutable bool isLocal_;
mutable bool isSend_;
// Patch index for coupled entities
mutable label patchIndex_;
// Master / slave indices
mutable label masterIndex_;
mutable label slaveIndex_;
// Point buffers
mutable pointField pointBuffer_;
mutable pointField oldPointBuffer_;
// List of points that are required
// during patch sub-mesh creation
mutable labelList subMeshPoints_;
mutable List<labelPair> globalProcPoints_;
// Entity sizes (as specified by the entityType enumerant)
mutable FixedList<label,10> nEntities_;
// Maps for entities
mutable FixedList<Map<label>,4> entityMap_;
mutable FixedList<Map<label>,4> reverseEntityMap_;
// Entity Buffers (as specified by the entityType enumerant)
mutable FixedList<labelList,12> entityBuffer_;
// List of entity indices with topological operations
mutable labelList entityIndices_;
// List of operations performed on entities
mutable List<opType> entityOperations_;
// List of point-locations to move points to
mutable pointField moveNewPoints_;
mutable pointField moveOldPoints_;
// Addressing for field mapping
mutable labelList faceMap_;
mutable labelList cellMap_;
//- Demand-driven connectivity data.
mutable edgeList* edgesPtr_;
mutable faceList* facesPtr_;
mutable labelListList* faceEdgesPtr_;
void makeEdges() const;
void makeFaces() const;
void makeFaceMap() const;
void makeCellMap() const;
void clearAddressing() const;
//- Disallow default bitwise assignment
void operator=(const coupleMap&);
public:
// Static data members
//- The set of names corresponding to the operationType enumeration
// Includes an extra entry for "invalid".
static const char* names[INVALID + 1];
//- Runtime type information
TypeName("coupleMap");
// Static Member Functions
//- Return a text representation of an opType
static const char* asText(const opType);
// Constructors
//- Construct from components
coupleMap
(
const IOobject& io,
const bool twoDMesh,
const bool isLocal,
const bool isSend,
const label patchIndex,
const label masterIndex,
const label slaveIndex
);
//- Construct as copy
coupleMap(const coupleMap&);
// Destructor
virtual ~coupleMap();
//- Interpolation functions
//- Interpolate point field
template<class Type>
tmp<Field<Type> > pointInterpolate
(
const Map<label>& mPointMap,
const Map<label>& sPointMap,
const Field<Type>& pf,
bool reverse = false
) const;
template<class Type>
tmp<Field<Type> > pointInterpolate
(
const Map<label>& mPointMap,
const Map<label>& sPointMap,
const tmp<Field<Type> >& tpf,
bool reverse = false
) const;
//- Interpolate face field
template<class Type>
tmp<Field<Type> > faceInterpolate
(
const label mStart,
const label sStart,
const Field<Type>& pf,
bool reverse = false
) const;
template<class Type>
tmp<Field<Type> > faceInterpolate
(
const label mStart,
const label sStart,
const tmp<Field<Type> >& tpf,
bool reverse = false
) const;
//- Access
inline label& patchIndex() const;
inline label masterIndex() const;
inline label slaveIndex() const;
inline bool isTwoDMesh() const;
inline bool isLocal() const;
inline bool isProcessor() const;
inline bool isSend() const;
inline bool isRecv() const;
pointField& pointBuffer() const;
pointField& oldPointBuffer() const;
labelList& subMeshPoints() const;
List<labelPair>& globalProcPoints() const;
void allocateBuffers() const;
label findSlave
(
const label eType,
const label Index
) const;
label findMaster
(
const label eType,
const label Index
) const;
void removeSlave
(
const label eType,
const label Index
) const;
void removeMaster
(
const label eType,
const label Index
) const;
void mapSlave
(
const label eType,
const label master,
const label slave
) const;
void mapMaster
(
const label eType,
const label slave,
const label master
) const;
inline FixedList<label,10>& nEntities() const;
inline label& nEntities(const label eType) const;
inline Map<label>& entityMap(const label eType) const;
inline Map<label>& reverseEntityMap(const label eType) const;
inline FixedList<labelList,12>& entityBuffer() const;
inline labelList& entityBuffer(const label eType) const;
inline labelList& entityIndices() const;
inline List<opType>& entityOperations() const;
inline pointField& moveNewPoints() const;
inline pointField& moveOldPoints() const;
void pushOperation
(
const label index,
const opType oType,
const point& newPoint = vector::zero,
const point& oldPoint = vector::zero
) const;
//- Demand-driven data
label nInternalFaces() const;
const labelList& owner() const;
const labelList& neighbour() const;
const edgeList& edges() const;
const faceList& faces() const;
const labelListList& faceEdges() const;
const labelList& faceMap() const;
const labelList& cellMap() const;
void transferMaps
(
const label eType,
Map<label>& newEntityMap,
Map<label>& newReverseEntityMap
) const;
void clearMaps() const;
void clearBuffers() const;
bool readData(Istream&);
bool writeData(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#include "coupleMapI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

288
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/coupleMapI.H Normal file
View file

@ -0,0 +1,288 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Interpolate point field
template<class Type>
tmp<Field<Type> >
coupleMap::pointInterpolate
(
const Map<label>& mPointMap,
const Map<label>& sPointMap,
const Field<Type>& pf,
bool reverse
) const
{
const Map<label>& pointMap = entityMap(coupleMap::POINT);
if (pf.size() != pointMap.size())
{
FatalErrorIn
(
"coupleMap::pointInterpolate"
"(const Field<Type> pf)"
) << "given field does not correspond to patch. Patch size: "
<< pointMap.size() << " field size: " << pf.size()
<< abort(FatalError);
}
tmp<Field<Type> > tresult
(
new Field<Type>
(
pointMap.size(),
pTraits<Type>::zero
)
);
Field<Type>& result = tresult();
if (reverse)
{
forAllConstIter(Map<label>, pointMap, pIter)
{
result[mPointMap[pIter.key()]] = pf[sPointMap[pIter()]];
}
}
else
{
forAllConstIter(Map<label>, pointMap, pIter)
{
result[sPointMap[pIter()]] = pf[mPointMap[pIter.key()]];
}
}
return tresult;
}
template<class Type>
tmp<Field<Type> >
coupleMap::pointInterpolate
(
const Map<label>& mPointMap,
const Map<label>& sPointMap,
const tmp<Field<Type> >& tpf,
bool reverse
) const
{
tmp<Field<Type> > tint =
(
pointInterpolate<Type>(mPointMap, sPointMap, tpf(), reverse)
);
tpf.clear();
return tint;
}
//- Interpolate point field
template<class Type>
tmp<Field<Type> >
coupleMap::faceInterpolate
(
const label mStart,
const label sStart,
const Field<Type>& pf,
bool reverse
) const
{
const Map<label>& faceMap = entityMap(coupleMap::FACE);
if (pf.size() != faceMap.size())
{
FatalErrorIn
(
"coupleMap::faceInterpolate"
"(const Field<Type> pf)"
) << "given field does not correspond to patch. Patch size: "
<< faceMap.size() << " field size: " << pf.size()
<< abort(FatalError);
}
tmp<Field<Type> > tresult
(
new Field<Type>
(
faceMap.size(),
pTraits<Type>::zero
)
);
Field<Type>& result = tresult();
if (reverse)
{
forAllConstIter(Map<label>, faceMap, fIter)
{
result[fIter.key() - mStart] = pf[fIter() - sStart];
}
}
else
{
forAllConstIter(Map<label>, faceMap, fIter)
{
result[fIter() - sStart] = pf[fIter.key() - mStart];
}
}
return tresult;
}
template<class Type>
tmp<Field<Type> >
coupleMap::faceInterpolate
(
const label mStart,
const label sStart,
const tmp<Field<Type> >& tpf,
bool reverse
) const
{
tmp<Field<Type> > tint =
(
faceInterpolate<Type>(mStart, sStart, tpf(), reverse)
);
tpf.clear();
return tint;
}
inline label& coupleMap::patchIndex() const
{
return patchIndex_;
}
inline label coupleMap::masterIndex() const
{
return masterIndex_;
}
inline label coupleMap::slaveIndex() const
{
return slaveIndex_;
}
inline bool coupleMap::isTwoDMesh() const
{
return twoDMesh_;
}
inline bool coupleMap::isLocal() const
{
return isLocal_;
}
inline bool coupleMap::isProcessor() const
{
return !isLocal_;
}
inline bool coupleMap::isSend() const
{
return isSend_;
}
inline bool coupleMap::isRecv() const
{
return !isSend_;
}
inline FixedList<label,10>& coupleMap::nEntities() const
{
return nEntities_;
}
inline label& coupleMap::nEntities(const label eType) const
{
return nEntities_[eType];
}
inline Map<label>& coupleMap::entityMap(const label eType) const
{
return entityMap_[eType];
}
inline Map<label>& coupleMap::reverseEntityMap(const label eType) const
{
return reverseEntityMap_[eType];
}
inline FixedList<labelList,12>& coupleMap::entityBuffer() const
{
return entityBuffer_;
}
inline labelList& coupleMap::entityBuffer(const label eType) const
{
return entityBuffer_[eType];
}
inline labelList& coupleMap::entityIndices() const
{
return entityIndices_;
}
inline List<coupleMap::opType>& coupleMap::entityOperations() const
{
return entityOperations_;
}
inline pointField& coupleMap::moveNewPoints() const
{
return moveNewPoints_;
}
inline pointField& coupleMap::moveOldPoints() const
{
return moveOldPoints_;
}
} // End namespace Foam
// ************************************************************************* //

648
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/coupledInfo.C Normal file
View file

@ -0,0 +1,648 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "Time.H"
#include "coupledInfo.H"
#include "dynamicTopoFvMesh.H"
#include "emptyFvPatchFields.H"
#include "emptyFvsPatchFields.H"
#include "fixedValueFvPatchFields.H"
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Constructor for coupledInfo
coupledInfo::coupledInfo
(
const dynamicTopoFvMesh& mesh,
const coupleMap& cMap,
const label mfzIndex,
const label sfzIndex
)
:
mesh_(mesh),
builtMaps_(false),
map_(cMap),
masterFaceZone_(mfzIndex),
slaveFaceZone_(sfzIndex)
{}
coupledInfo::coupledInfo
(
const dynamicTopoFvMesh& mesh,
const bool isTwoDMesh,
const bool isLocal,
const bool isSend,
const label patchIndex,
const label mPatch,
const label sPatch,
const label mfzIndex,
const label sfzIndex
)
:
mesh_(mesh),
builtMaps_(false),
map_
(
IOobject
(
"coupleMap_"
+ Foam::name(mPatch)
+ "_To_"
+ Foam::name(sPatch)
+ word(isLocal ? "_Local" : "_Proc")
+ word(isSend ? "_Send" : "_Recv"),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
true
),
isTwoDMesh,
isLocal,
isSend,
patchIndex,
mPatch,
sPatch
),
masterFaceZone_(mfzIndex),
slaveFaceZone_(sfzIndex)
{}
//- Construct given addressing
coupledInfo::subMeshMapper::subMeshMapper
(
const coupledInfo& cInfo,
const label patchI
)
:
sizeBeforeMapping_(cInfo.baseMesh().boundary()[patchI].size()),
directAddressing_
(
SubList<label>
(
cInfo.map().faceMap(),
cInfo.subMesh().boundary()[patchI].size(),
cInfo.subMesh().boundary()[patchI].patch().start()
)
)
{
// Offset indices
label pStart = cInfo.baseMesh().boundary()[patchI].patch().start();
forAll(directAddressing_, faceI)
{
directAddressing_[faceI] -= pStart;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const dynamicTopoFvMesh& coupledInfo::baseMesh() const
{
return mesh_;
}
void coupledInfo::setMesh
(
label index,
dynamicTopoFvMesh* mesh
)
{
subMesh_.set(mesh);
}
dynamicTopoFvMesh& coupledInfo::subMesh()
{
if (!subMesh_.valid())
{
FatalErrorIn("dynamicTopoFvMesh& coupledInfo::subMesh()")
<< " Sub-mesh pointer has not been set."
<< abort(FatalError);
}
return subMesh_();
}
const dynamicTopoFvMesh& coupledInfo::subMesh() const
{
if (!subMesh_.valid())
{
FatalErrorIn("const dynamicTopoFvMesh& coupledInfo::subMesh() const")
<< " Sub-mesh pointer has not been set."
<< abort(FatalError);
}
return subMesh_();
}
bool coupledInfo::builtMaps() const
{
return builtMaps_;
}
void coupledInfo::setBuiltMaps()
{
builtMaps_ = true;
}
coupleMap& coupledInfo::map()
{
return map_;
}
const coupleMap& coupledInfo::map() const
{
return map_;
}
label coupledInfo::masterFaceZone() const
{
return masterFaceZone_;
}
label coupledInfo::slaveFaceZone() const
{
return slaveFaceZone_;
}
// Set subMesh centres
void coupledInfo::setCentres(PtrList<volVectorField>& centres) const
{
// Fetch reference to subMesh
const dynamicTopoFvMesh& mesh = subMesh();
// Set size
centres.setSize(1);
vectorField Cv(mesh.cellCentres());
vectorField Cf(mesh.faceCentres());
// Create and map the patch field values
label nPatches = mesh.boundary().size();
// Create field parts
PtrList<fvPatchField<vector> > volCentrePatches(nPatches);
// Over-ride and set all patches to fixedValue
for (label patchI = 0; patchI < nPatches; patchI++)
{
volCentrePatches.set
(
patchI,
new fixedValueFvPatchField<vector>
(
mesh.boundary()[patchI],
DimensionedField<vector, volMesh>::null()
)
);
// Slice field to patch (forced assignment)
volCentrePatches[patchI] ==
(
mesh.boundaryMesh()[patchI].patchSlice(Cf)
);
}
// Set the cell-centres pointer.
centres.set
(
0,
new volVectorField
(
IOobject
(
"cellCentres",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimLength,
SubField<vector>(Cv, mesh.nCells()),
volCentrePatches
)
);
}
// Subset volume field
template <class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
coupledInfo::subSetVolField
(
const GeometricField<Type, fvPatchField, volMesh>& fld
) const
{
// Create and map the internal-field values
Field<Type> internalField
(
fld.internalField(),
map().cellMap()
);
// Create and map the patch field values
label nPatches = subMesh().boundary().size();
PtrList<fvPatchField<Type> > patchFields(nPatches);
forAll(patchFields, patchI)
{
if (patchI == (nPatches - 1))
{
// Artificially set last patch
patchFields.set
(
patchI,
new emptyFvPatchField<Type>
(
subMesh().boundary()[patchI],
DimensionedField<Type, volMesh>::null()
)
);
}
else
{
patchFields.set
(
patchI,
fvPatchField<Type>::New
(
fld.boundaryField()[patchI],
subMesh().boundary()[patchI],
DimensionedField<Type, volMesh>::null(),
subMeshMapper(*this, patchI)
)
);
}
}
// Create new field from pieces
tmp<GeometricField<Type, fvPatchField, volMesh> > subFld
(
new GeometricField<Type, fvPatchField, volMesh>
(
IOobject
(
"subField_" + fld.name(),
subMesh().time().timeName(),
subMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
subMesh(),
fld.dimensions(),
internalField,
patchFields
)
);
return subFld;
}
// Subset surface field
template <class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
coupledInfo::subSetSurfaceField
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& fld
) const
{
// Create and map the internal-field values
Field<Type> internalField
(
fld.internalField(),
SubList<label>
(
map().faceMap(),
subMesh().nInternalFaces()
)
);
// Create and map the patch field values
label nPatches = subMesh().boundary().size();
PtrList<fvsPatchField<Type> > patchFields(nPatches);
forAll(patchFields, patchI)
{
if (patchI == (nPatches - 1))
{
// Artificially set last patch
patchFields.set
(
patchI,
new emptyFvsPatchField<Type>
(
subMesh().boundary()[patchI],
DimensionedField<Type, surfaceMesh>::null()
)
);
}
else
{
patchFields.set
(
patchI,
fvsPatchField<Type>::New
(
fld.boundaryField()[patchI],
subMesh().boundary()[patchI],
DimensionedField<Type, surfaceMesh>::null(),
subMeshMapper(*this, patchI)
)
);
}
}
// Create new field from pieces
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > subFld
(
new GeometricField<Type, fvsPatchField, surfaceMesh>
(
IOobject
(
"subField_" + fld.name(),
subMesh().time().timeName(),
subMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
subMesh(),
fld.dimensions(),
internalField,
patchFields
)
);
return subFld;
}
template <class Type>
void coupledInfo::mapVolField
(
const wordList& fieldNames,
const word& fieldType,
OSstream& strStream
) const
{
strStream
<< fieldType << token::NL
<< token::BEGIN_BLOCK << token::NL;
forAll(fieldNames, i)
{
const GeometricField<Type, fvPatchField, volMesh>& fld =
(
mesh_.lookupObject
<
GeometricField<Type, fvPatchField, volMesh>
>(fieldNames[i])
);
tmp<GeometricField<Type, fvPatchField, volMesh> > tsubFld =
(
subSetVolField(fld)
);
// Send field through stream
strStream
<< fieldNames[i]
<< token::NL << token::BEGIN_BLOCK
<< tsubFld
<< token::NL << token::END_BLOCK
<< token::NL;
}
strStream
<< token::END_BLOCK << token::NL;
}
template <class Type>
void coupledInfo::mapSurfaceField
(
const wordList& fieldNames,
const word& fieldType,
OSstream& strStream
) const
{
strStream
<< fieldType << token::NL
<< token::BEGIN_BLOCK << token::NL;
forAll(fieldNames, i)
{
const GeometricField<Type, fvsPatchField, surfaceMesh>& fld =
(
mesh_.lookupObject
<
GeometricField<Type, fvsPatchField, surfaceMesh>
>(fieldNames[i])
);
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsubFld =
(
subSetSurfaceField(fld)
);
// Send field through stream
strStream
<< fieldNames[i]
<< token::NL << token::BEGIN_BLOCK
<< tsubFld
<< token::NL << token::END_BLOCK
<< token::NL;
}
strStream
<< token::END_BLOCK << token::NL;
}
// Set volume field pointer from input dictionary
template <class GeomField>
void coupledInfo::setField
(
const wordList& fieldNames,
const dictionary& fieldDicts,
PtrList<GeomField>& fields
) const
{
// Size up the pointer list
fields.setSize(fieldNames.size());
forAll(fieldNames, i)
{
fields.set
(
i,
new GeomField
(
IOobject
(
fieldNames[i],
subMesh().time().timeName(),
subMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
subMesh(),
fieldDicts.subDict(fieldNames[i])
)
);
}
}
template <class GeomField>
void coupledInfo::resizeMap
(
const label srcIndex,
const subMeshMapper& internalMapper,
const List<labelList>& internalReverseMaps,
const PtrList<subMeshMapper>& boundaryMapper,
const List<labelListList>& boundaryReverseMaps,
const List<PtrList<GeomField> >& srcFields,
GeomField& field
)
{
// autoMap the internal field
field.internalField().autoMap(internalMapper);
// Reverse map for additional cells
forAll(srcFields, pI)
{
// Fetch field for this processor
const GeomField& srcField = srcFields[pI][srcIndex];
field.internalField().rmap
(
srcField.internalField(),
internalReverseMaps[pI]
);
}
// Map physical boundary-fields
forAll(boundaryMapper, patchI)
{
// autoMap the patchField
field.boundaryField()[patchI].autoMap(boundaryMapper[patchI]);
// Reverse map for additional patch faces
forAll(srcFields, pI)
{
// Fetch field for this processor
const GeomField& srcField = srcFields[pI][srcIndex];
field.boundaryField()[patchI].rmap
(
srcField.boundaryField()[patchI],
boundaryReverseMaps[pI][patchI]
);
}
}
}
// Resize all fields in registry
template <class GeomField>
void coupledInfo::resizeMap
(
const wordList& names,
const objectRegistry& mesh,
const subMeshMapper& internalMapper,
const List<labelList>& internalReverseMaps,
const PtrList<subMeshMapper>& boundaryMapper,
const List<labelListList>& boundaryReverseMaps,
const List<PtrList<GeomField> >& srcFields
)
{
forAll(names, indexI)
{
// Fetch field from registry
GeomField& field =
(
const_cast<GeomField&>
(
mesh.lookupObject<GeomField>(names[indexI])
)
);
// Map the field
coupledInfo::resizeMap
(
indexI,
internalMapper,
internalReverseMaps,
boundaryMapper,
boundaryReverseMaps,
srcFields,
field
);
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void coupledInfo::operator=(const coupledInfo& rhs)
{
// Check for assignment to self
if (this == &rhs)
{
FatalErrorIn
(
"void coupledInfo::operator=(const Foam::coupledInfo&)"
)
<< "Attempted assignment to self"
<< abort(FatalError);
}
}
} // End namespace Foam
// ************************************************************************* //

282
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/coupledInfo.H Normal file
View file

@ -0,0 +1,282 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
coupledInfo
Description
An interface class that provides patch coupling functionality
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
coupledInfo.C
\*---------------------------------------------------------------------------*/
#ifndef coupledInfo_H
#define coupledInfo_H
#include "autoPtr.H"
#include "coupleMap.H"
#include "volFieldsFwd.H"
#include "surfaceFieldsFwd.H"
#include "fvPatchFieldMapper.H"
namespace Foam
{
// Class forward declarations
class dynamicTopoFvMesh;
/*---------------------------------------------------------------------------*\
Class coupledInfo Declaration
\*---------------------------------------------------------------------------*/
class coupledInfo
{
// Reference to the parent mesh
const dynamicTopoFvMesh& mesh_;
// Auto pointer to a subMesh
autoPtr<dynamicTopoFvMesh> subMesh_;
// Flag to determine whether maps have been built.
bool builtMaps_;
// For locally coupled patches,
// specify a master / slave index
coupleMap map_;
// Zone IDs for patches associated with faceZones
label masterFaceZone_;
label slaveFaceZone_;
//- Disallow default bitwise assignment
inline void operator=(const coupledInfo&);
public:
// Constructor
inline coupledInfo
(
const dynamicTopoFvMesh& mesh,
const coupleMap& cMap,
const label mfzIndex = -1,
const label sfzIndex = -1
);
inline coupledInfo
(
const dynamicTopoFvMesh& mesh,
const bool isTwoDMesh,
const bool isLocal,
const bool isSend,
const label patchIndex,
const label mPatch,
const label sPatch,
const label mfzIndex = -1,
const label sfzIndex = -1
);
//- Access
// Return a const reference to the parent mesh
inline const dynamicTopoFvMesh& baseMesh() const;
// Set a new subMesh
inline void setMesh(label index, dynamicTopoFvMesh* mesh);
// Return a reference to the subMesh
inline dynamicTopoFvMesh& subMesh();
// Return a const reference to the subMesh
inline const dynamicTopoFvMesh& subMesh() const;
// Have maps been built?
inline bool builtMaps() const;
// Change the flag
inline void setBuiltMaps();
// Return a reference to the coupleMap
inline coupleMap& map();
// Return a const reference to the coupleMap
inline const coupleMap& map() const;
// Return the master / slave face zone IDs
inline label masterFaceZone() const;
inline label slaveFaceZone() const;
//- Interpolation
//- Generic subMesh mapper
class subMeshMapper
:
public fvPatchFieldMapper
{
label sizeBeforeMapping_;
labelField directAddressing_;
public:
// Constructors
//- Construct from components
inline subMeshMapper
(
const label sbm,
const labelList& da
)
:
sizeBeforeMapping_(sbm),
directAddressing_(da)
{}
//- Construct given addressing
inline subMeshMapper
(
const coupledInfo& cInfo,
const label patchI
);
// Destructor
virtual ~subMeshMapper()
{}
// Member Functions
label size() const
{
return directAddressing_.size();
}
label sizeBeforeMapping() const
{
return directAddressing_.size();
}
bool direct() const
{
return true;
}
const unallocLabelList& directAddressing() const
{
return directAddressing_;
}
};
// Set subMesh centres
inline void setCentres(PtrList<volVectorField>& centres) const;
// Subset volume field
template <class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
subSetVolField
(
const GeometricField<Type, fvPatchField, volMesh>& field
) const;
// Subset surface field
template <class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
subSetSurfaceField
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& field
) const;
// Subset the volume fields from registry to output stream
template <class Type>
void mapVolField
(
const wordList& fieldNames,
const word& fieldType,
OSstream& strStream
) const;
// Subset the surface field from registry to output stream
template <class Type>
void mapSurfaceField
(
const wordList& fieldNames,
const word& fieldType,
OSstream& strStream
) const;
// Set volume field pointer from input dictionary
template <class GeomField>
void setField
(
const wordList& fieldNames,
const dictionary& fieldDicts,
PtrList<GeomField>& fields
) const;
// Resize map for individual field
template <class GeomField>
static void resizeMap
(
const label srcIndex,
const subMeshMapper& internalMapper,
const List<labelList>& internalReverseMaps,
const PtrList<subMeshMapper>& boundaryMapper,
const List<labelListList>& boundaryReverseMaps,
const List<PtrList<GeomField> >& srcFields,
GeomField& field
);
// Resize map for all fields in registry
template <class GeomField>
static void resizeMap
(
const wordList& fieldNames,
const objectRegistry& mesh,
const subMeshMapper& internalMapper,
const List<labelList>& internalReverseMaps,
const PtrList<subMeshMapper>& boundaryMapper,
const List<labelListList>& boundaryReverseMaps,
const List<PtrList<GeomField> >& srcFields
);
};
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "coupledInfo.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

1902
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.C
View file

File diff suppressed because it is too large Load diff

386
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.H
View file

@ -37,6 +37,7 @@ SourceFiles
dynamicTopoFvMesh.C
dynamicTopoFvMeshI.H
dynamicTopoFvMeshCheck.C
dynamicTopoFvMeshCoupled.C
dynamicTopoFvMeshReOrder.C
dynamicTopoFvMeshMapping.C
edgeBisect.C
@ -50,6 +51,7 @@ SourceFiles
#include "Switch.H"
#include "tetMetric.H"
#include "topoMapper.H"
#include "DynamicField.H"
#include "threadHandler.H"
#include "dynamicFvMesh.H"
@ -64,8 +66,9 @@ class eMesh;
class Stack;
class changeMap;
class objectMap;
class topoMapper;
class coupledInfo;
class motionSolver;
class convexSetAlgorithm;
class lengthScaleEstimator;
/*---------------------------------------------------------------------------*\
@ -78,9 +81,18 @@ class dynamicTopoFvMesh
{
// Private data
//- Reference to base mesh
const dynamicTopoFvMesh& baseMesh_;
//- Number of boundary patches
label nPatches_;
//- Topology change flag
bool topoChangeFlag_;
//- Is this a subMesh?
bool isSubMesh_;
//- Dynamic mesh dictionary
IOdictionary dict_;
@ -94,9 +106,15 @@ class dynamicTopoFvMesh
//- Edge refinement switch
Switch edgeRefinement_;
//- Load motion solver
Switch loadMotionSolver_;
//- Switch for cell-bandwidth reduction
Switch bandWidthReduction_;
//- Coupled modification switch
mutable Switch coupledModification_;
//- Specify the re-meshing interval
label interval_;
@ -119,27 +137,21 @@ class dynamicTopoFvMesh
// support template typedefs,
// so this is a work-around
template <class T>
class resizableList
class resizable
{
public:
typedef DynamicList<T, 0, 11, 10> Type;
typedef DynamicList<T, 0, 11, 10> ListType;
typedef DynamicField<T, 0, 11, 10> FieldType;
};
template <class T>
class resizableField
{
public:
typedef DynamicField<T> Type;
};
resizableField<point>::Type oldPoints_, points_;
resizableList<face>::Type faces_;
resizableList<label>::Type owner_, neighbour_;
resizableList<cell>::Type cells_;
resizableList<edge>::Type edges_;
resizableList<labelList>::Type pointEdges_, edgePoints_;
resizableList<labelList>::Type edgeFaces_, faceEdges_;
resizableList<scalar>::Type lengthScale_;
resizable<point>::FieldType oldPoints_, points_;
resizable<face>::ListType faces_;
resizable<label>::ListType owner_, neighbour_;
resizable<cell>::ListType cells_;
resizable<edge>::ListType edges_;
resizable<labelList>::ListType pointEdges_;
resizable<labelList>::ListType edgeFaces_, faceEdges_;
resizable<scalar>::ListType lengthScale_;
//- Size information
labelList oldPatchSizes_, patchSizes_;
@ -181,6 +193,7 @@ class dynamicTopoFvMesh
Map<labelList> faceParents_;
List<scalarField> faceWeights_;
List<vectorField> faceCentres_;
Map<labelList> cellParents_;
List<scalarField> cellWeights_;
List<vectorField> cellCentres_;
@ -201,9 +214,8 @@ class dynamicTopoFvMesh
labelHashSet deletedFaces_;
labelHashSet deletedCells_;
//- List of flipped faces / modified old-points
//- List of flipped faces
labelHashSet flipFaces_;
Map<labelList> modPoints_;
//- Run-time statistics
label maxModifications_;
@ -237,11 +249,42 @@ class dynamicTopoFvMesh
// in multi-threaded reOrdering
FixedList<Mutex, 4> entityMutex_;
// Local coupled patch information
PtrList<coupledInfo> patchCoupling_;
// List of processors that this sub-domain talks to
labelList procIndices_;
// Sub-Mesh pointers
PtrList<coupledInfo> sendMeshes_;
PtrList<coupledInfo> recvMeshes_;
// Private Member Functions
//- Disallow default bitwise copy construct
dynamicTopoFvMesh(const dynamicTopoFvMesh&);
//- Construct from components
//- Used for subMeshes only.
dynamicTopoFvMesh
(
const dynamicTopoFvMesh& mesh,
const IOobject& io,
const Xfer<pointField>& points,
const Xfer<pointField>& oldPoints,
const Xfer<edgeList>& edges,
const Xfer<faceList>& faces,
const Xfer<labelListList>& faceEdges,
const Xfer<labelList>& owner,
const Xfer<labelList>& neighbour,
const labelList& faceStarts,
const labelList& faceSizes,
const labelList& edgeStarts,
const labelList& edgeSizes,
const wordList& patchNames,
const dictionary& patchDict
);
//- Disallow default bitwise assignment
void operator=(const dynamicTopoFvMesh&);
@ -256,6 +299,7 @@ class dynamicTopoFvMesh
(
const scalar matchTol,
const bool skipMapping,
const bool mappingOutput,
const label faceStart,
const label faceSize,
const label cellStart,
@ -265,19 +309,13 @@ class dynamicTopoFvMesh
// Static equivalent for multiThreading
static void computeMappingThread(void *argument);
// Compute parents for inverse-distance weighting
void computeParents
(
const label index,
const labelList& mapCandidates,
const labelListList& oldNeighbourList,
const vectorField& oldCentres,
const label dimension,
labelList& parents
) const;
// Routine to invoke threaded mapping
void threadedMapping(scalar matchTol, bool skipMapping);
void threadedMapping
(
scalar matchTol,
bool skipMapping,
bool mappingOutput
);
// Initialize mesh edges and related connectivity lists
void initEdges();
@ -316,11 +354,21 @@ class dynamicTopoFvMesh
// contains a triangular face.
inline label getTriBoundaryEdge(const label fIndex) const;
// Given an input quad-face, determine checkEdges from mesh
static void getCheckEdges
(
const label fIndex,
const dynamicTopoFvMesh& mesh,
changeMap& map,
FixedList<edge,4>& checkEdge,
FixedList<label,4>& checkEdgeIndex
);
// Return length-scale at an face-location in the mesh [2D]
inline scalar faceLengthScale(const label fIndex) const;
scalar faceLengthScale(const label fIndex) const;
// Return length-scale at an edge-location in the mesh [3D]
inline scalar edgeLengthScale(const label eIndex) const;
scalar edgeLengthScale(const label eIndex) const;
// Check for bisection
inline bool checkBisection(const label index) const;
@ -349,6 +397,13 @@ class dynamicTopoFvMesh
// Initialize stacks
inline void initStacks(const labelHashSet& entities);
// Initialize the coupled stack
void initCoupledStack
(
const labelHashSet& entities,
bool useEntities
);
// Method to determine the boundary patch index for a given face
inline label whichPatch(const label index) const;
@ -385,13 +440,6 @@ class dynamicTopoFvMesh
bool forceOp = false
);
// Method for the swapping of an edge in 3D
void swapEdge
(
const label eIndex,
bool forceOp = false
);
// Method for the bisection of an edge in 3D
const changeMap
bisectEdge
@ -420,6 +468,9 @@ class dynamicTopoFvMesh
bool forceOp = false
);
// Slice the mesh at a particular location
void sliceMesh(const labelPair& pointPair);
// Utility method to check for invalid face-bisection.
bool checkBisection
(
@ -429,8 +480,9 @@ class dynamicTopoFvMesh
) const;
// Utility method to check for invalid face-collapse.
bool checkCollapse
static bool checkCollapse
(
const dynamicTopoFvMesh& mesh,
const labelList& triFaces,
const FixedList<label,2>& c0BdyIndex,
const FixedList<label,2>& c1BdyIndex,
@ -440,7 +492,7 @@ class dynamicTopoFvMesh
scalar& collapseQuality,
const bool checkNeighbour,
bool forceOp = false
) const;
);
// Utility method to check for invalid edge-collapse.
bool checkCollapse
@ -449,7 +501,7 @@ class dynamicTopoFvMesh
const point& oldPoint,
const label pointIndex,
const label cellIndex,
labelHashSet& cellsChecked,
DynamicList<label>& cellsChecked,
scalar& collapseQuality,
bool forceOp = false
) const;
@ -491,7 +543,9 @@ class dynamicTopoFvMesh
removeCells
(
const labelList& cList,
const label patch
const label patch,
const word& rcName,
bool checkOnly = false
);
// Remove the specified cell from the mesh
@ -525,13 +579,19 @@ class dynamicTopoFvMesh
const Map<bool>& cellColors
);
// Merge a set of boundary faces into internal
const changeMap
mergeBoundaryFaces
(
const labelList& mergeFaces
);
// Insert the specified edge to the mesh
label insertEdge
(
const label patch,
const edge& newEdge,
const labelList& edgeFaces,
const labelList& edgePoints = labelList::null()
const labelList& edgeFaces
);
// Remove the specified edge from the mesh
@ -555,12 +615,13 @@ class dynamicTopoFvMesh
const label pIndex
);
// Utility method to build edgePoints for an edge [3D]
void buildEdgePoints
// Utility method to build vertexHull for an edge [3D]
void buildVertexHull
(
const label eIndex,
labelList& vertexHull,
const label checkIndex = 0
);
) const;
// Utility to check whether points of an edge lie on a boundary
const FixedList<bool,2> checkEdgeBoundary(const label eIndex) const;
@ -569,7 +630,20 @@ class dynamicTopoFvMesh
inline void setFlip(const label fIndex);
// Utility method to compute the minimum quality of a vertex hull
scalar computeMinQuality(const label eIndex) const;
scalar computeMinQuality
(
const label eIndex,
labelList& hullVertices
) const;
// Compute minQuality given addressing
scalar computeMinQuality
(
const edge& edgeToCheck,
const labelList& hullVertices,
const UList<point>& points,
bool closedRing = false
) const;
// Utility method to compute the quality of a vertex hull
// around an edge after bisection.
@ -580,7 +654,12 @@ class dynamicTopoFvMesh
scalar computeTrisectionQuality(const label fIndex) const;
// Check whether the given edge is on a bounding curve
bool checkBoundingCurve(const label eIndex) const;
bool checkBoundingCurve
(
const label eIndex,
const bool overRidePurityCheck = false,
label* nProcCurves = NULL
) const;
// Allocate dynamic programming tables
void initTables
@ -606,14 +685,28 @@ class dynamicTopoFvMesh
bool fillTables
(
const label eIndex,
const scalar minQuality,
scalar& minQuality,
labelList& m,
labelList& hullVertices,
PtrList<scalarListList>& Q,
PtrList<labelListList>& K,
PtrList<labelListList>& triangulations,
const label checkIndex = 0
) const;
// Fill tables given addressing
void fillTables
(
const edge& edgeToCheck,
const scalar minQuality,
const label m,
const labelList& hullVertices,
const UList<point>& points,
scalarListList& Q,
labelListList& K,
labelListList& triangulations
) const;
// Print out tables for debugging
void printTables
(
@ -629,7 +722,9 @@ class dynamicTopoFvMesh
(
const label eIndex,
const scalar minQuality,
const PtrList<labelListList>& K,
const labelList& hullVertices,
PtrList<scalarListList>& Q,
PtrList<labelListList>& K,
PtrList<labelListList>& triangulations,
const label checkIndex = 0
);
@ -649,7 +744,8 @@ class dynamicTopoFvMesh
(
const label eIndex,
const labelList& hullVertices,
const labelListList& triangulations
const labelListList& triangulations,
bool output = false
) const;
// Check old-volumes for the configuration.
@ -781,6 +877,9 @@ class dynamicTopoFvMesh
// Reset the mesh and generate mapping information
bool resetMesh();
// Write out connectivity for an edge
void writeEdgeConnectivity(const label eIndex) const;
// Output an entity as a VTK file
void writeVTK
(
@ -798,7 +897,9 @@ class dynamicTopoFvMesh
const labelList& cList,
const label primitiveType = 3,
const bool useOldConnectivity = false,
const bool useOldPoints = false
const bool useOldPoints = false,
const UList<scalar>& scalField = UList<scalar>(),
const UList<label>& lablField = UList<label>()
) const;
// Return the status report interval
@ -807,6 +908,18 @@ class dynamicTopoFvMesh
// Check the state of local connectivity lists
void checkConnectivity(const label maxErrors = 0) const;
// Handle mesh slicing events.
void handleMeshSlicing();
// Handle layer addition / removal events
void handleLayerAdditionRemoval();
// Add cell layer above specified patch
const changeMap addCellLayer(const label patchID);
// Remove cell layer above specified patch
const changeMap removeCellLayer(const label patchID);
// Test an edge / face for proximity with other non-neighbouring faces.
// Return the scalar distance to the nearest face.
scalar testProximity
@ -818,9 +931,165 @@ class dynamicTopoFvMesh
// Calculate the edge length-scale for the mesh
void calculateLengthScale(bool dump = false);
// Fill buffers with length-scale info
// and exchange across processors.
void exchangeLengthBuffers();
// Implementing the fillTables operation for coupled edges
bool coupledFillTables
(
const label eIndex,
scalar& minQuality,
labelList& m,
PtrList<scalarListList>& Q,
PtrList<labelListList>& K,
PtrList<labelListList>& triangulations
) const;
// Additional mapping contributions for coupled entities
void computeCoupledWeights
(
const label index,
const label dimension,
labelList& parents,
scalarField& weights,
vectorField& centres,
bool output = false
);
// Fetch length-scale info for processor entities
scalar processorLengthScale(const label index) const;
// Method to determine whether the master entity is locally coupled
bool locallyCoupledEntity
(
const label index,
bool checkSlaves = false,
bool checkProcs = true,
bool checkFace = false
) const;
// Method to determine the locally coupled patch index
label locallyCoupledEdgePatch(const label eIndex) const;
// Method to determine if the entity is on a processor boundary
bool processorCoupledEntity
(
const label index,
bool checkFace = false,
bool checkEdge = false,
bool checkPure = false,
FixedList<label, 2>* patchLabels = NULL,
FixedList<vector, 2>* patchNormals = NULL
) const;
// Read optional dictionary parameters
void readOptionalParameters(bool reRead = false);
// Execute load balancing operations on the mesh
void executeLoadBalancing(const dictionary& balanceDict);
// Identify coupled patches.
bool identifyCoupledPatches();
// Read coupled patch information from dictionary.
void readCoupledPatches();
// Initialize coupled patch connectivity for topology modifications.
static void initCoupledConnectivity(void *argument);
// Move coupled subMesh points
void moveCoupledSubMeshes();
// Build a list of entities that need to be avoided
// by regular topo-changes.
void buildEntitiesToAvoid
(
labelHashSet& entities,
bool checkSubMesh
);
// Check whether the specified edge is a coupled master edge.
bool isCoupledMaster(const label eIndex) const;
// Set coupled modification
void setCoupledModification() const;
// Unset coupled modification
void unsetCoupledModification() const;
// Insert the cells around the coupled master entity to the mesh
const changeMap insertCells(const label mIndex);
// Handle topology changes for coupled patches
void handleCoupledPatches(labelHashSet& entities);
// Synchronize topology operations across processors
void syncCoupledPatches(labelHashSet& entities);
// Check the state of coupled boundaries
bool checkCoupledBoundaries(bool report = true) const;
// Build patch sub-meshes for processor patches
void buildProcessorPatchMeshes();
// Build patch sub-mesh for a specified processor
void buildProcessorPatchMesh
(
coupledInfo& subMesh,
Map<labelList>& commonCells
);
// Build coupled maps for locally coupled patches
void buildLocalCoupledMaps();
// Build coupled maps for coupled processor patches
void buildProcessorCoupledMaps();
// Introduce a new processor patch to the mesh
label createProcessorPatch(const label proc);
// If a patch is of processor type,
// get the neighbouring processor ID
label getNeighbourProcessor(const label patch) const;
// If the number of patches have changed during run-time,
// reset boundaries with new processor patches
void resetBoundaries();
// Initialize subMesh field transfers for mapping
void initFieldTransfers
(
wordList& fieldTypes,
List<wordList>& fieldNames,
List<List<char> >& sendBuffer,
List<List<char> >& recvBuffer
);
// Synchronize field transfers for mapping
void syncFieldTransfers
(
wordList& fieldTypes,
List<wordList>& fieldNames,
List<List<char> >& recvBuffer
);
// Initialize coupled boundary ordering
void initCoupledBoundaryOrdering
(
List<pointField>& centres,
List<pointField>& anchors
) const;
// Synchronize coupled boundary ordering
bool syncCoupledBoundaryOrdering
(
List<pointField>& centres,
List<pointField>& anchors,
labelListList& faceMaps,
labelListList& rotations
) const;
// Dump cell-quality statistics
bool meshQuality(bool outputOption);
@ -840,8 +1109,11 @@ public:
// Member Functions
//- Update mesh corresponding to the given map
virtual void updateMesh(const mapPolyMesh& mpm);
// Map all fields in time using a customized mapper
virtual void mapFields(const mapPolyMesh& meshMap) const;
virtual void mapFields(const mapPolyMesh& mpm) const;
// Update the mesh for motion / topology changes
// Return true if topology changes have occurred

1284
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C
View file

File diff suppressed because it is too large Load diff

10088
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshCoupled.C Normal file
View file

File diff suppressed because it is too large Load diff

354
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshI.H
View file

@ -199,272 +199,6 @@ inline label dynamicTopoFvMesh::getTriBoundaryEdge
}
// Return length-scale at an face-location in the mesh [2D]
inline scalar dynamicTopoFvMesh::faceLengthScale
(
const label fIndex
) const
{
// Reset the scale first
scalar scale = 0.0;
label facePatch = whichPatch(fIndex);
// Determine whether the face is internal
if (facePatch < 0)
{
# ifdef FULLDEBUG
// Check whether neighbour is valid
if (neighbour_[fIndex] == -1)
{
FatalErrorIn
(
"inline scalar dynamicTopoFvMesh::faceLengthScale"
"(const label fIndex) const"
)
<< nl << "Face: " << fIndex
<< ": " << faces_[fIndex]
<< " is not internal."
<< abort(FatalError);
}
# endif
scale =
(
0.5 *
(
lengthScale_[owner_[fIndex]]
+ lengthScale_[neighbour_[fIndex]]
)
);
}
else
{
// Fetch the fixed-length scale
scale = lengthEstimator().fixedLengthScale(fIndex, facePatch);
// If this is a floating face, pick the owner length-scale
if (lengthEstimator().isFreePatch(facePatch))
{
scale = lengthScale_[owner_[fIndex]];
}
// If proximity-based refinement is requested,
// test the proximity to the nearest non-neighbour.
if (lengthEstimator().isProximityPatch(facePatch))
{
label proximityFace = -1;
// Perform a proximity-check.
scalar distance = testProximity(fIndex, proximityFace);
if (debug > 3 && self() == 0)
{
if
(
(proximityFace > -1) &&
((distance / 5.0) < scale)
)
{
Info << " Closest opposing face detected for face: " << nl
<< '\t' << fIndex
<< " :: " << faces_[fIndex]
<< " was face:\n"
<< '\t' << proximityFace
<< " :: " << polyMesh::faces()[proximityFace] << nl
<< " with distance: " << distance
<< endl;
}
}
scale =
(
Foam::min
(
scale,
((distance / 3.0) - SMALL)/lengthEstimator().ratioMax()
)
);
}
// Limit scales if necessary
lengthEstimator().limitScale(scale);
}
return scale;
}
// Compute length-scale at an edge-location in the mesh [3D]
inline scalar dynamicTopoFvMesh::edgeLengthScale
(
const label eIndex
) const
{
// Reset the scale first
scalar scale = 0.0;
const labelList& eFaces = edgeFaces_[eIndex];
label edgePatch = whichEdgePatch(eIndex);
// Determine whether the edge is internal
if (edgePatch < 0)
{
forAll(eFaces, faceI)
{
# ifdef FULLDEBUG
// Check whether neighbour is valid
if (neighbour_[eFaces[faceI]] == -1)
{
FatalErrorIn
(
"inline scalar dynamicTopoFvMesh::edgeLengthScale"
"(const label eIndex) const"
)
<< nl << "Face: " << eFaces[faceI]
<< ": " << faces_[eFaces[faceI]]
<< " is not internal, while edge: "
<< eIndex << ": " << edges_[eIndex] << " is."
<< abort(FatalError);
}
# endif
scale += lengthScale_[owner_[eFaces[faceI]]];
scale += lengthScale_[neighbour_[eFaces[faceI]]];
}
scale /= (2.0*eFaces.size());
}
else
{
// Search for boundary faces, and average their scale
forAll(eFaces, faceI)
{
if (neighbour_[eFaces[faceI]] == -1)
{
scale +=
(
lengthEstimator().fixedLengthScale
(
eFaces[faceI],
edgePatch
)
);
}
}
scale *= 0.5;
// If proximity-based refinement is requested,
// test the proximity to the nearest non-neighbour.
if (lengthEstimator().isProximityPatch(edgePatch))
{
label proximityFace = -1;
// Perform a proximity-check.
scalar distance = testProximity(eIndex, proximityFace);
if (debug > 3 && self() == 0)
{
if
(
(proximityFace > -1) &&
((distance / 5.0) < scale)
)
{
Info << " Closest opposing face detected for edge: " << nl
<< '\t' << eIndex
<< " :: " << edges_[eIndex]
<< " was face:\n"
<< '\t' << proximityFace
<< " :: " << polyMesh::faces()[proximityFace] << nl
<< " with distance: " << distance
<< endl;
}
}
scale =
(
Foam::min
(
scale,
((distance / 3.0) - SMALL)/lengthEstimator().ratioMax()
)
);
}
// If curvature-based refinement is requested,
// test the variation in face-normal directions.
if (lengthEstimator().isCurvaturePatch(edgePatch))
{
// Obtain face-normals for both faces.
label count = 0;
FixedList<vector, 2> fNorm;
forAll(eFaces, faceI)
{
if (neighbour_[eFaces[faceI]] == -1)
{
// Obtain the normal.
fNorm[count] = faces_[eFaces[faceI]].normal(points_);
// Normalize it.
fNorm[count] /= mag(fNorm[count]);
count++;
}
}
scalar deviation = (fNorm[0] & fNorm[1]);
scalar refDeviation = lengthEstimator().curvatureDeviation();
if (mag(deviation) < refDeviation)
{
// Fetch the edge
const edge& edgeToCheck = edges_[eIndex];
// Get the edge-length.
scalar length =
(
linePointRef
(
points_[edgeToCheck.start()],
points_[edgeToCheck.end()]
).mag()
);
if (debug > 3 && self() == 0)
{
Info << "Deviation: " << deviation << nl
<< "curvatureDeviation: " << refDeviation
<< ", Edge: " << eIndex << ", Length: " << length
<< ", Scale: " << scale << nl
<< " Half-length: " << (0.5*length) << nl
<< " MinRatio: "
<< (lengthEstimator().ratioMin()*scale)
<< endl;
}
scale =
(
Foam::min
(
scale,
((length - SMALL)/lengthEstimator().ratioMax())
)
);
}
}
// Limit scales if necessary
lengthEstimator().limitScale(scale);
}
return scale;
}
// Check for edge bisection
inline bool dynamicTopoFvMesh::checkBisection
(
@ -475,6 +209,12 @@ inline bool dynamicTopoFvMesh::checkBisection
if (twoDMesh_)
{
// If this entity was deleted, skip it.
if (faces_[index].empty())
{
return false;
}
// Fetch the edge
const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
@ -493,6 +233,12 @@ inline bool dynamicTopoFvMesh::checkBisection
}
else
{
// If this entity was deleted, skip it.
if (edgeFaces_[index].empty())
{
return false;
}
// Fetch the edge
const edge& edgeToCheck = edges_[index];
@ -531,6 +277,12 @@ inline bool dynamicTopoFvMesh::checkCollapse
if (twoDMesh_)
{
// If this entity was deleted, skip it.
if (faces_[index].empty())
{
return false;
}
// Fetch the edge
const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
@ -549,6 +301,12 @@ inline bool dynamicTopoFvMesh::checkCollapse
}
else
{
// If this entity was deleted, skip it.
if (edgeFaces_[index].empty())
{
return false;
}
// Fetch the edge
const edge& edgeToCheck = edges_[index];
@ -637,6 +395,15 @@ inline void dynamicTopoFvMesh::initStacks
{
forAll(faces_, faceI)
{
// For coupled meshes, avoid certain faces.
if (patchCoupling_.size() || procIndices_.size())
{
if (entities.found(faceI))
{
continue;
}
}
if (faces_[faceI].size() == 4)
{
stack(tID[tIndex]).insert(faceI);
@ -649,6 +416,15 @@ inline void dynamicTopoFvMesh::initStacks
{
forAll(edges_, edgeI)
{
// For coupled meshes, avoid certain edges.
if (patchCoupling_.size() || procIndices_.size())
{
if (entities.found(edgeI))
{
continue;
}
}
if (edgeFaces_[edgeI].size())
{
stack(tID[tIndex]).insert(edgeI);
@ -657,6 +433,32 @@ inline void dynamicTopoFvMesh::initStacks
}
}
}
if (debug > 3 && Pstream::parRun())
{
Pout<< nl << "Stack size: " << stack(0).size() << endl;
if (debug > 4)
{
// Write out stack entities
labelList stackElements(stack(0).size(), -1);
forAll(stackElements, elemI)
{
stackElements[elemI] = stack(0)[elemI];
}
label elemType = twoDMesh_ ? 2 : 1;
writeVTK
(
"Stack_"
+ Foam::name(Pstream::myProcNo()),
stackElements,
elemType
);
}
}
}
@ -686,9 +488,11 @@ inline label dynamicTopoFvMesh::whichPatch
// If not in any of the above, it's possible that the face was added
// at the end of the list. Check addedFacePatches_ for the patch info
if (addedFacePatches_.found(index))
Map<label>::const_iterator it = addedFacePatches_.find(index);
if (it != addedFacePatches_.end())
{
return addedFacePatches_[index];
return it();
}
else
{
@ -737,9 +541,11 @@ inline label dynamicTopoFvMesh::whichEdgePatch
// If not in any of the above, it's possible that the edge was added
// at the end of the list. Check addedEdgePatches_ for the patch info
if (addedEdgePatches_.found(index))
Map<label>::const_iterator it = addedEdgePatches_.find(index);
if (it != addedEdgePatches_.end())
{
return addedEdgePatches_[index];
return it();
}
else
{
@ -799,13 +605,15 @@ inline void dynamicTopoFvMesh::setFlip(const label fIndex)
{
if (fIndex < nOldFaces_)
{
if (flipFaces_.found(fIndex))
labelHashSet::iterator it = flipFaces_.find(fIndex);
if (it == flipFaces_.end())
{
flipFaces_.erase(fIndex);
flipFaces_.insert(fIndex);
}
else
{
flipFaces_.insert(fIndex);
flipFaces_.erase(it);
}
}
}

620
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshMapping.C
View file

@ -37,14 +37,21 @@ Author
#include "dynamicTopoFvMesh.H"
#include "Time.H"
#include "meshOps.H"
#include "IOmanip.H"
#include "triFace.H"
#include "objectMap.H"
#include "faceSetAlgorithm.H"
#include "cellSetAlgorithm.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTemplateTypeNameAndDebugWithName(IOList<objectMap>, "objectMapIOList", 0);
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Compute mapping weights for modified entities
@ -52,77 +59,397 @@ void dynamicTopoFvMesh::computeMapping
(
const scalar matchTol,
const bool skipMapping,
const bool mappingOutput,
const label faceStart,
const label faceSize,
const label cellStart,
const label cellSize
)
{
// Convex-set algorithm for cells
cellSetAlgorithm cellAlgorithm
(
(*this),
oldPoints_,
edges_,
faces_,
cells_,
owner_,
neighbour_
);
label nInconsistencies = 0;
scalar maxFaceError = 0.0, maxCellError = 0.0;
DynamicList<scalar> cellErrors(10), faceErrors(10);
DynamicList<objectMap> failedCells(10), failedFaces(10);
// Compute cell mapping
for (label cellI = cellStart; cellI < (cellStart + cellSize); cellI++)
{
label cIndex = cellsFromCells_[cellI].index();
labelList& masterObjects = cellsFromCells_[cellI].masterObjects();
if (skipMapping)
{
// Set empty mapping parameters
const labelList& mo = cellParents_[cIndex];
cellsFromCells_[cellI].masterObjects() = mo;
cellWeights_[cellI].setSize(mo.size(), (1.0/(mo.size() + VSMALL)));
cellCentres_[cellI].setSize(mo.size(), vector::zero);
// Dummy map from cell[0]
masterObjects = labelList(1, 0);
cellWeights_[cellI].setSize(1, 1.0);
cellCentres_[cellI].setSize(1, vector::zero);
}
else
{
// Compute master objects for inverse-distance weighting
computeParents
// Obtain weighting factors for this cell.
cellAlgorithm.computeWeights
(
cIndex,
0,
cellParents_[cIndex],
polyMesh::cellCells(),
polyMesh::cellCentres(),
3,
cellsFromCells_[cellI].masterObjects()
masterObjects,
cellWeights_[cellI],
cellCentres_[cellI]
);
// Add contributions from subMeshes, if any.
computeCoupledWeights
(
cIndex,
cellAlgorithm.dimension(),
masterObjects,
cellWeights_[cellI],
cellCentres_[cellI]
);
// Compute error
scalar error = mag(1.0 - sum(cellWeights_[cellI]));
if (error > matchTol)
{
bool consistent = false;
// Check whether any edges lie on boundary patches.
// These cells can have relaxed weights to account
// for mild convexity.
const cell& cellToCheck = cells_[cIndex];
if (twoDMesh_)
{
const labelList& parents = cellParents_[cIndex];
forAll(parents, cI)
{
const cell& pCell = polyMesh::cells()[parents[cI]];
forAll(pCell, fI)
{
const face& pFace = polyMesh::faces()[pCell[fI]];
if (pFace.size() == 3)
{
continue;
}
label fP = boundaryMesh().whichPatch(pCell[fI]);
// Disregard processor patches
if (getNeighbourProcessor(fP) > -1)
{
continue;
}
if (fP > -1)
{
consistent = true;
break;
}
}
if (consistent)
{
break;
}
}
}
else
{
forAll(cellToCheck, fI)
{
const labelList& fE = faceEdges_[cellToCheck[fI]];
forAll(fE, eI)
{
label eP = whichEdgePatch(fE[eI]);
// Disregard processor patches
if (getNeighbourProcessor(eP) > -1)
{
continue;
}
if (eP > -1)
{
consistent = true;
break;
}
}
if (consistent)
{
break;
}
}
}
if (!consistent)
{
nInconsistencies++;
// Add to list
cellErrors.append(error);
// Accumulate error stats
maxCellError = Foam::max(maxCellError, error);
failedCells.append
(
objectMap
(
cIndex,
cellParents_[cIndex]
)
);
}
}
}
}
// Convex-set algorithm for faces
faceSetAlgorithm faceAlgorithm
(
(*this),
oldPoints_,
edges_,
faces_,
cells_,
owner_,
neighbour_
);
// Compute face mapping
for (label faceI = faceStart; faceI < (faceStart + faceSize); faceI++)
{
label fIndex = facesFromFaces_[faceI].index();
label patchIndex = whichPatch(fIndex);
labelList& masterObjects = facesFromFaces_[faceI].masterObjects();
// Skip mapping for internal faces.
if (patchIndex == -1)
label patchIndex = whichPatch(fIndex);
label neiProc = getNeighbourProcessor(patchIndex);
// Skip mapping for internal / processor faces.
if (patchIndex == -1 || neiProc > -1)
{
// Set dummy masters, so that the conventional
// faceMapper doesn't incur a seg-fault.
facesFromFaces_[faceI].masterObjects() = labelList(1, 0);
// faceMapper doesn't crash-and-burn
masterObjects = labelList(1, 0);
continue;
}
if (skipMapping)
{
// Set empty mapping parameters
const labelList& mo = faceParents_[fIndex];
facesFromFaces_[faceI].masterObjects() = mo;
faceWeights_[faceI].setSize(mo.size(), (1.0/(mo.size() + VSMALL)));
faceCentres_[faceI].setSize(mo.size(), vector::zero);
// Dummy map from patch[0]
masterObjects = labelList(1, 0);
faceWeights_[faceI].setSize(1, 1.0);
faceCentres_[faceI].setSize(1, vector::zero);
}
else
{
// Compute master objects for inverse-distance weighting
computeParents
// Obtain weighting factors for this face.
faceAlgorithm.computeWeights
(
fIndex,
boundaryMesh()[patchIndex].start(),
faceParents_[fIndex],
boundaryMesh()[patchIndex].faceFaces(),
boundaryMesh()[patchIndex].faceCentres(),
2,
facesFromFaces_[faceI].masterObjects()
masterObjects,
faceWeights_[faceI],
faceCentres_[faceI]
);
// Add contributions from subMeshes, if any.
computeCoupledWeights
(
fIndex,
faceAlgorithm.dimension(),
masterObjects,
faceWeights_[faceI],
faceCentres_[faceI]
);
// Compute error
scalar error = mag(1.0 - sum(faceWeights_[faceI]));
if (error > matchTol)
{
bool consistent = false;
// Check whether any edges lie on bounding curves.
// These faces can have relaxed weights to account
// for addressing into patches on the other side
// of the curve.
const labelList& fEdges = faceEdges_[fIndex];
forAll(fEdges, eI)
{
if (checkBoundingCurve(fEdges[eI]))
{
consistent = true;
}
}
if (!consistent)
{
nInconsistencies++;
// Add to list
faceErrors.append(error);
// Accumulate error stats
maxFaceError = Foam::max(maxFaceError, error);
failedFaces.append
(
objectMap
(
fIndex,
faceParents_[fIndex]
)
);
}
}
}
}
if (nInconsistencies)
{
Pout<< " Mapping errors: "
<< " max cell error: " << maxCellError
<< " max face error: " << maxFaceError
<< endl;
if (debug || mappingOutput)
{
if (failedCells.size())
{
Pout<< " failedCells: " << endl;
forAll(failedCells, cellI)
{
label index = failedCells[cellI].index();
Pout<< " Cell: " << index
<< " Error: " << cellErrors[cellI]
<< endl;
}
}
if (failedFaces.size())
{
Pout<< " failedFaces: " << endl;
forAll(failedFaces, faceI)
{
label index = failedFaces[faceI].index();
label patchIndex = whichPatch(index);
const polyBoundaryMesh& boundary = boundaryMesh();
Pout<< " Face: " << index
<< " Patch: " << boundary[patchIndex].name()
<< " Error: " << faceErrors[faceI]
<< endl;
}
}
if (debug > 3 || mappingOutput)
{
// Prepare lists
labelList objects;
scalarField weights;
vectorField centres;
if (failedCells.size())
{
forAll(failedCells, cellI)
{
label cIndex = failedCells[cellI].index();
cellAlgorithm.computeWeights
(
cIndex,
0,
cellParents_[cIndex],
polyMesh::cellCells(),
objects,
weights,
centres,
true
);
computeCoupledWeights
(
cIndex,
cellAlgorithm.dimension(),
objects,
weights,
centres,
true
);
// Clear lists
objects.clear();
weights.clear();
centres.clear();
}
}
if (failedFaces.size())
{
forAll(failedFaces, faceI)
{
label fIndex = failedFaces[faceI].index();
label patchIndex = whichPatch(fIndex);
const polyBoundaryMesh& boundary = boundaryMesh();
faceAlgorithm.computeWeights
(
fIndex,
boundary[patchIndex].start(),
faceParents_[fIndex],
boundary[patchIndex].faceFaces(),
objects,
weights,
centres,
true
);
computeCoupledWeights
(
fIndex,
faceAlgorithm.dimension(),
objects,
weights,
centres,
true
);
// Clear lists
objects.clear();
weights.clear();
centres.clear();
}
}
}
}
}
}
@ -144,16 +471,18 @@ void dynamicTopoFvMesh::computeMappingThread(void *argument)
// Recast the pointers for the argument
scalar& matchTol = *(static_cast<scalar*>(thread->operator()(0)));
bool& skipMapping = *(static_cast<bool*>(thread->operator()(1)));
label& faceStart = *(static_cast<label*>(thread->operator()(2)));
label& faceSize = *(static_cast<label*>(thread->operator()(3)));
label& cellStart = *(static_cast<label*>(thread->operator()(4)));
label& cellSize = *(static_cast<label*>(thread->operator()(5)));
bool& mappingOutput = *(static_cast<bool*>(thread->operator()(2)));
label& faceStart = *(static_cast<label*>(thread->operator()(3)));
label& faceSize = *(static_cast<label*>(thread->operator()(4)));
label& cellStart = *(static_cast<label*>(thread->operator()(5)));
label& cellSize = *(static_cast<label*>(thread->operator()(6)));
// Now calculate addressing
mesh.computeMapping
(
matchTol,
skipMapping,
mappingOutput,
faceStart, faceSize,
cellStart, cellSize
);
@ -169,7 +498,8 @@ void dynamicTopoFvMesh::computeMappingThread(void *argument)
void dynamicTopoFvMesh::threadedMapping
(
scalar matchTol,
bool skipMapping
bool skipMapping,
bool mappingOutput
)
{
label nThreads = threader_->getNumThreads();
@ -177,7 +507,7 @@ void dynamicTopoFvMesh::threadedMapping
// If mapping is being skipped, issue a warning.
if (skipMapping)
{
Info << " *** Mapping is being skipped *** " << endl;
Info<< " *** Mapping is being skipped *** " << endl;
}
// Check if single-threaded
@ -187,6 +517,7 @@ void dynamicTopoFvMesh::threadedMapping
(
matchTol,
skipMapping,
mappingOutput,
0, facesFromFaces_.size(),
0, cellsFromCells_.size()
);
@ -212,8 +543,8 @@ void dynamicTopoFvMesh::threadedMapping
if (debug > 2)
{
Info << " Mapping Faces: " << index[0] << endl;
Info << " Mapping Cells: " << index[1] << endl;
Pout<< " Mapping Faces: " << index[0] << nl
<< " Mapping Cells: " << index[1] << endl;
}
forAll(index, indexI)
@ -234,8 +565,8 @@ void dynamicTopoFvMesh::threadedMapping
if (debug > 2)
{
Info << " Load starts: " << tStarts[indexI] << endl;
Info << " Load sizes: " << tSizes[indexI] << endl;
Pout<< " Load starts: " << tStarts[indexI] << nl
<< " Load sizes: " << tSizes[indexI] << endl;
}
}
@ -243,7 +574,7 @@ void dynamicTopoFvMesh::threadedMapping
forAll(hdl, i)
{
// Size up the argument list
hdl[i].setSize(6);
hdl[i].setSize(7);
// Set match tolerance
hdl[i].set(0, &matchTol);
@ -251,25 +582,26 @@ void dynamicTopoFvMesh::threadedMapping
// Set the skipMapping flag
hdl[i].set(1, &skipMapping);
// Set the mappingOutput flag
hdl[i].set(2, &mappingOutput);
// Set the start/size indices
hdl[i].set(2, &(tStarts[0][i]));
hdl[i].set(3, &(tSizes[0][i]));
hdl[i].set(4, &(tStarts[1][i]));
hdl[i].set(5, &(tSizes[1][i]));
hdl[i].set(3, &(tStarts[0][i]));
hdl[i].set(4, &(tSizes[0][i]));
hdl[i].set(5, &(tStarts[1][i]));
hdl[i].set(6, &(tSizes[1][i]));
}
// Prior to multi-threaded operation,
// force calculation of demand-driven data.
polyMesh::cells();
primitiveMesh::cellCells();
primitiveMesh::cellCentres();
const polyBoundaryMesh& boundary = boundaryMesh();
forAll(boundary, patchI)
{
boundary[patchI].faceFaces();
boundary[patchI].faceCentres();
}
// Execute threads in linear sequence
@ -277,130 +609,6 @@ void dynamicTopoFvMesh::threadedMapping
}
// Compute parents for inverse-distance weighting
void dynamicTopoFvMesh::computeParents
(
const label index,
const labelList& mapCandidates,
const labelListList& oldNeighbourList,
const vectorField& oldCentres,
const label dimension,
labelList& parents
) const
{
if (parents.size())
{
FatalErrorIn
(
"\n\n"
"void dynamicTopoFvMesh::computeParents\n"
"(\n"
" const label index,\n"
" const labelList& mapCandidates,\n"
" const labelListList& oldNeighbourList,\n"
" const vectorField& oldCentres,\n"
" const label dimension,\n"
" labelList& parents\n"
") const\n"
)
<< " Addressing has already been calculated." << nl
<< " Index: " << index << nl
<< " Type: " << (dimension == 2 ? "Face" : "Cell") << nl
<< " mapCandidates: " << mapCandidates << nl
<< " Parents: " << parents << nl
<< abort(FatalError);
}
// Figure out the patch offset and centre
label offset = -1;
vector centre = vector::zero;
if (dimension == 2)
{
offset = boundaryMesh()[whichPatch(index)].start();
centre = faces_[index].centre(oldPoints_);
}
else
if (dimension == 3)
{
offset = 0;
scalar dummyVol = 0.0;
meshOps::cellCentreAndVolume
(
index,
oldPoints_,
faces_,
cells_,
owner_,
centre,
dummyVol
);
}
// Maintain a check-list
labelHashSet checked;
// Insert candidates first
forAll(mapCandidates, indexI)
{
checked.insert(mapCandidates[indexI] - offset);
}
// Loop for three outward levels from candidates
for (label i = 0; i < 3; i++)
{
// Fetch the set of candidates
const labelList checkList = checked.toc();
forAll(checkList, indexI)
{
const labelList& oldNei = oldNeighbourList[checkList[indexI]];
forAll(oldNei, entityI)
{
if (!checked.found(oldNei[entityI]))
{
checked.insert(oldNei[entityI]);
}
}
}
}
// Loop through accumulated candidates and fetch the nearest one.
scalar minDist = GREAT;
label minLocation = -1;
const labelList checkList = checked.toc();
forAll(checkList, indexI)
{
scalar dist = magSqr(oldCentres[checkList[indexI]] - centre);
if (dist < minDist)
{
minDist = dist;
minLocation = checkList[indexI];
}
}
// Set the final list
const labelList& neiList = oldNeighbourList[minLocation];
parents.setSize(neiList.size() + 1, -1);
// Fill indices
forAll(neiList, indexI)
{
parents[indexI] = neiList[indexI] + offset;
}
// Set final index
parents[neiList.size()] = minLocation + offset;
}
// Set fill-in mapping information for a particular cell
void dynamicTopoFvMesh::setCellMapping
(
@ -413,9 +621,9 @@ void dynamicTopoFvMesh::setCellMapping
{
if (debug > 3)
{
Info << "Inserting mapping cell: " << cIndex << nl
<< " mapCells: " << mapCells
<< endl;
Pout<< "Inserting mapping cell: " << cIndex
<< " mapCells: " << mapCells
<< endl;
}
// Insert index into the list, and overwrite if necessary
@ -445,7 +653,7 @@ void dynamicTopoFvMesh::setCellMapping
}
// Update cell-parents information
labelHashSet masterCells;
DynamicList<label> masterCells(5);
forAll(mapCells, cellI)
{
@ -456,7 +664,10 @@ void dynamicTopoFvMesh::setCellMapping
if (mapCells[cellI] < nOldCells_)
{
masterCells.insert(mapCells[cellI]);
if (findIndex(masterCells, mapCells[cellI]) == -1)
{
masterCells.append(mapCells[cellI]);
}
}
else
if (cellParents_.found(mapCells[cellI]))
@ -465,12 +676,15 @@ void dynamicTopoFvMesh::setCellMapping
forAll(nParents, cI)
{
masterCells.insert(nParents[cI]);
if (findIndex(masterCells, nParents[cI]) == -1)
{
masterCells.append(nParents[cI]);
}
}
}
}
cellParents_.set(cIndex, masterCells.toc());
cellParents_.set(cIndex, masterCells);
}
@ -482,13 +696,33 @@ void dynamicTopoFvMesh::setFaceMapping
)
{
label patch = whichPatch(fIndex);
label neiProc = getNeighbourProcessor(patch);
if (debug > 3)
{
Info << "Inserting mapping face: " << fIndex << nl
<< " patch: " << patch << nl
<< " mapFaces: " << mapFaces
<< endl;
const polyBoundaryMesh& boundary = boundaryMesh();
word pName;
if (patch == -1)
{
pName = "Internal";
}
else
if (patch < boundary.size())
{
pName = boundaryMesh()[patch].name();
}
else
{
pName = "New patch: " + Foam::name(patch);
}
Pout<< "Inserting mapping face: " << fIndex
<< " patch: " << pName
<< " mapFaces: " << mapFaces
<< " neiProc: " << neiProc
<< endl;
}
bool foundError = false;
@ -500,15 +734,16 @@ void dynamicTopoFvMesh::setFaceMapping
foundError = true;
}
// Check to ensure that boundary faces map
// only from other faces on the same patch
if (patch > -1 && mapFaces.empty())
// Check to ensure that mapFaces is non-empty for physical patches
if (patch > -1 && mapFaces.empty() && neiProc == -1)
{
foundError = true;
}
if (foundError)
{
writeVTK("mapFace_" + Foam::name(fIndex), fIndex, 2);
FatalErrorIn
(
"\n"
@ -524,7 +759,8 @@ void dynamicTopoFvMesh::setFaceMapping
<< " 2. Mapping specified for an internal face, " << nl
<< " when none was expected." << nl << nl
<< " Face: " << fIndex << nl
<< " Patch: " << patch << nl
<< " Patch: "
<< (patch > -1 ? boundaryMesh()[patch].name() : "Internal") << nl
<< " Owner: " << owner_[fIndex] << nl
<< " Neighbour: " << neighbour_[fIndex] << nl
<< " mapFaces: " << mapFaces << nl
@ -556,14 +792,14 @@ void dynamicTopoFvMesh::setFaceMapping
facesFromFaces_[index].masterObjects() = labelList(0);
}
// For internal faces, set dummy maps / weights, and bail out
if (patch == -1)
// For internal / processor faces, bail out
if (patch == -1 || neiProc > -1)
{
return;
}
// Update face-parents information
labelHashSet masterFaces;
DynamicList<label> masterFaces(5);
forAll(mapFaces, faceI)
{
@ -574,7 +810,10 @@ void dynamicTopoFvMesh::setFaceMapping
if (mapFaces[faceI] < nOldFaces_)
{
masterFaces.insert(mapFaces[faceI]);
if (findIndex(masterFaces, mapFaces[faceI]) == -1)
{
masterFaces.append(mapFaces[faceI]);
}
}
else
if (faceParents_.found(mapFaces[faceI]))
@ -583,12 +822,15 @@ void dynamicTopoFvMesh::setFaceMapping
forAll(nParents, fI)
{
masterFaces.insert(nParents[fI]);
if (findIndex(masterFaces, nParents[fI]) == -1)
{
masterFaces.append(nParents[fI]);
}
}
}
}
faceParents_.set(fIndex, masterFaces.toc());
faceParents_.set(fIndex, masterFaces);
}

491
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshReOrder.C
View file

@ -25,6 +25,7 @@ License
\*---------------------------------------------------------------------------*/
#include "objectMap.H"
#include "coupledInfo.H"
#include "dynamicTopoFvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -51,10 +52,10 @@ void dynamicTopoFvMesh::reOrderPoints
{
if (threaded)
{
Info << "Thread: " << self() << ": ";
Info<< "Thread: " << self() << ": ";
}
Info << "ReOrdering points..." << flush;
Info<< "ReOrdering points..." << flush;
}
// Allocate for the mapping information
@ -204,6 +205,52 @@ void dynamicTopoFvMesh::reOrderPoints
// Reset all zones
pointZones.updateMesh();
// Loop through all local coupling maps, and renumber points.
forAll(patchCoupling_, patchI)
{
if (!patchCoupling_(patchI))
{
continue;
}
const coupleMap& cMap = patchCoupling_[patchI].map();
// Obtain references
Map<label>& mtsMap = cMap.entityMap(coupleMap::POINT);
Map<label> newMtsMap, newStmMap;
forAllIter(Map<label>, mtsMap, pIter)
{
label newMaster = -1, newSlave = -1;
if (pIter.key() < nOldPoints_)
{
newMaster = reversePointMap_[pIter.key()];
}
else
{
newMaster = addedPointRenumbering_[pIter.key()];
}
if (pIter() < nOldPoints_)
{
newSlave = reversePointMap_[pIter()];
}
else
{
newSlave = addedPointRenumbering_[pIter()];
}
// Update the map.
newMtsMap.insert(newMaster, newSlave);
newStmMap.insert(newSlave, newMaster);
}
// Overwrite the old maps.
cMap.transferMaps(coupleMap::POINT, newMtsMap, newStmMap);
}
// Clear local point copies
points_.clear();
oldPoints_.clear();
@ -214,7 +261,7 @@ void dynamicTopoFvMesh::reOrderPoints
if (debug)
{
Info << "Done." << endl;
Info<< "Done." << endl;
}
}
@ -274,10 +321,10 @@ void dynamicTopoFvMesh::reOrderEdges
{
if (threaded)
{
Info << "Thread: " << self() << ": ";
Info<< "Thread: " << self() << ": ";
}
Info << "ReOrdering edges..." << flush;
Info<< "ReOrdering edges..." << flush;
}
// Allocate for mapping information
@ -286,10 +333,8 @@ void dynamicTopoFvMesh::reOrderEdges
label edgeInOrder = 0, allEdges = edges_.size();
edgeList oldEdges(allEdges);
labelListList oldEdgeFaces(allEdges);
labelListList oldEdgePoints(allEdges);
addedEdgeRenumbering_.clear();
Map<label> addedEdgeReverseRenumbering;
// Transfer old edge-based lists, and clear them
forAll(edges_, edgeI)
@ -300,20 +345,10 @@ void dynamicTopoFvMesh::reOrderEdges
edges_.setSize(nEdges_); edgeFaces_.setSize(nEdges_);
if (!twoDMesh_)
{
forAll(edgePoints_, edgeI)
{
oldEdgePoints[edgeI].transfer(edgePoints_[edgeI]);
}
edgePoints_.setSize(nEdges_);
}
// Keep track of inserted boundary edge indices
labelList boundaryPatchIndices(edgePatchStarts_);
// Loop through all edges and add internal ones first
// Loop through all edges and add them
forAll(oldEdges, edgeI)
{
// Ensure that we're adding valid edges
@ -325,8 +360,10 @@ void dynamicTopoFvMesh::reOrderEdges
// Determine which patch this edge belongs to
label patch = whichEdgePatch(edgeI);
// Update maps for boundary edges. Edge insertion for
// boundaries will be done after internal edges.
// Obtain references
edge& thisEdge = oldEdges[edgeI];
labelList& thisEF = oldEdgeFaces[edgeI];
if (patch >= 0)
{
label bEdgeIndex = boundaryPatchIndices[patch]++;
@ -339,18 +376,20 @@ void dynamicTopoFvMesh::reOrderEdges
}
else
{
addedEdgeRenumbering_.insert(edgeI, bEdgeIndex);
addedEdgeReverseRenumbering.insert(bEdgeIndex, edgeI);
edgeMap_[bEdgeIndex] = -1;
addedEdgeRenumbering_.insert(edgeI, bEdgeIndex);
}
// Insert entities into local lists...
edges_[bEdgeIndex] = thisEdge;
edgeFaces_[bEdgeIndex] = thisEF;
// Insert entities into mesh-reset lists
edges[bEdgeIndex] = thisEdge;
edgeFaces[bEdgeIndex].transfer(thisEF);
}
else
{
// Obtain references
edge& thisEdge = oldEdges[edgeI];
labelList& thisEF = oldEdgeFaces[edgeI];
// Renumber internal edges and add normally.
if (edgeI < nOldEdges_)
{
@ -370,46 +409,10 @@ void dynamicTopoFvMesh::reOrderEdges
edges[edgeInOrder] = thisEdge;
edgeFaces[edgeInOrder].transfer(thisEF);
if (!twoDMesh_)
{
edgePoints_[edgeInOrder].transfer(oldEdgePoints[edgeI]);
}
edgeInOrder++;
}
}
// All internal edges have been inserted. Now insert boundary edges.
label oldIndex;
for(label i = nInternalEdges_; i < nEdges_; i++)
{
if (edgeMap_[i] == -1)
{
// This boundary edge was added during the topology change
oldIndex = addedEdgeReverseRenumbering[i];
}
else
{
oldIndex = edgeMap_[i];
}
// Insert entities into local Lists...
edges_[edgeInOrder] = oldEdges[oldIndex];
edgeFaces_[edgeInOrder] = oldEdgeFaces[oldIndex];
// Insert entities into mesh-reset lists
edges[edgeInOrder] = oldEdges[oldIndex];
edgeFaces[edgeInOrder].transfer(oldEdgeFaces[oldIndex]);
if (!twoDMesh_)
{
edgePoints_[edgeInOrder].transfer(oldEdgePoints[oldIndex]);
}
edgeInOrder++;
}
// Now that we're done with edges, unlock it
if (threaded)
{
@ -417,15 +420,15 @@ void dynamicTopoFvMesh::reOrderEdges
}
// Final check to ensure everything went okay
if (edgeInOrder != nEdges_)
if (edgeInOrder != nInternalEdges_)
{
FatalErrorIn("dynamicTopoFvMesh::reOrderEdges()") << nl
<< " Algorithm did not visit every edge in the mesh."
<< " Something's messed up." << nl
<< abort(FatalError);
<< " Algorithm did not visit every internal edge in the mesh."
<< " Something's messed up." << nl
<< abort(FatalError);
}
// Renumber all edges / edgePoints with updated point information
// Renumber all edges with updated point information
label pIndex = -1;
if (threaded)
@ -463,24 +466,6 @@ void dynamicTopoFvMesh::reOrderEdges
thisEdge[1] = pIndex;
thisREdge[1] = pIndex;
// Renumber edgePoints
if (!twoDMesh_)
{
labelList& ePoints = edgePoints_[edgeI];
forAll(ePoints, pointI)
{
if (ePoints[pointI] < nOldPoints_)
{
ePoints[pointI] = reversePointMap_[ePoints[pointI]];
}
else
{
ePoints[pointI] = addedPointRenumbering_[ePoints[pointI]];
}
}
}
}
if (threaded)
@ -553,11 +538,56 @@ void dynamicTopoFvMesh::reOrderEdges
{
// Invert edges to obtain pointEdges
pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_);
// Loop through all local coupling maps, and renumber edges.
forAll(patchCoupling_, patchI)
{
if (!patchCoupling_(patchI))
{
continue;
}
// Obtain references
const coupleMap& cMap = patchCoupling_[patchI].map();
const Map<label>& mtsMap = cMap.entityMap(coupleMap::EDGE);
Map<label> newMtsMap, newStmMap;
forAllConstIter(Map<label>, mtsMap, mIter)
{
label newMaster = -1, newSlave = -1;
if (mIter.key() < nOldEdges_)
{
newMaster = reverseEdgeMap_[mIter.key()];
}
else
{
newMaster = addedEdgeRenumbering_[mIter.key()];
}
if (mIter() < nOldEdges_)
{
newSlave = reverseEdgeMap_[mIter()];
}
else
{
newSlave = addedEdgeRenumbering_[mIter()];
}
// Update the map.
newMtsMap.insert(newMaster, newSlave);
newStmMap.insert(newSlave, newMaster);
}
// Overwrite the old maps.
cMap.transferMaps(coupleMap::EDGE, newMtsMap, newStmMap);
}
}
if (debug)
{
Info << "Done." << endl;
Info<< "Done." << endl;
}
}
@ -624,10 +654,10 @@ void dynamicTopoFvMesh::reOrderFaces
{
if (threaded)
{
Info << "Thread: " << self() << ": ";
Info<< "Thread: " << self() << ": ";
}
Info << "ReOrdering faces..." << flush;
Info<< "ReOrdering faces..." << flush;
}
// Allocate for mapping information
@ -639,6 +669,9 @@ void dynamicTopoFvMesh::reOrderFaces
labelListList oldFaceEdges(allFaces);
addedFaceRenumbering_.clear();
// Track reverse renumbering for added faces.
// - Required during coupled patch re-ordering.
Map<label> addedFaceReverseRenumbering;
// Make a copy of the old face-based lists, and clear them
@ -709,7 +742,6 @@ void dynamicTopoFvMesh::reOrderFaces
// Handle boundaries first. If any coupled interfaces need to be
// updated, they can be reshuffled after interior faces are done.
// Update maps for boundaries now.
for (label faceI = nOldInternalFaces_; faceI < allFaces; faceI++)
{
if (visited[faceI] == -1)
@ -725,18 +757,46 @@ void dynamicTopoFvMesh::reOrderFaces
}
else
{
addedFaceRenumbering_.insert(faceI, bFaceIndex);
addedFaceReverseRenumbering.insert(bFaceIndex, faceI);
faceMap_[bFaceIndex] = -1;
addedFaceRenumbering_.insert(faceI, bFaceIndex);
addedFaceReverseRenumbering.insert(bFaceIndex, faceI);
}
// Renumber owner
label ownerRenumber =
(
oldOwner[faceI] < nOldCells_
? reverseCellMap_[oldOwner[faceI]]
: addedCellRenumbering_[oldOwner[faceI]]
);
// Insert entities into local lists...
owner_[bFaceIndex] = ownerRenumber;
neighbour_[bFaceIndex] = -1;
faces_[bFaceIndex] = oldFaces[faceI];
faceEdges_[bFaceIndex] = oldFaceEdges[faceI];
// Insert entities into mesh-reset lists...
owner[bFaceIndex] = ownerRenumber;
faces[bFaceIndex].transfer(oldFaces[faceI]);
faceEdges[bFaceIndex].transfer(oldFaceEdges[faceI]);
// Mark this face as visited
visited[faceI] = 0;
}
}
// Prepare centres and anchors for coupled interfaces
List<pointField> centres(nPatches_), anchors(nPatches_);
// Now that points and boundary faces are up-to-date,
// send and receive centres and anchor points for coupled patches
initCoupledBoundaryOrdering
(
centres,
anchors
);
// Upper-triangular ordering of internal faces:
// Insertion cannot be done in one go as the faces need to be
@ -885,43 +945,105 @@ void dynamicTopoFvMesh::reOrderFaces
}
}
// All internal faces have been inserted. Now insert boundary faces.
label oldIndex;
// Prepare faceMaps and rotations for coupled interfaces
labelListList patchMaps(nPatches_), rotations(nPatches_);
for (label i = nInternalFaces_; i < nFaces_; i++)
{
if (faceMap_[i] == -1)
{
// This boundary face was added during the topology change
oldIndex = addedFaceReverseRenumbering[i];
}
else
{
oldIndex = faceMap_[i];
}
// Renumber owner/neighbour
label ownerRenumber =
// Now compute patchMaps for coupled interfaces
bool anyChange =
(
syncCoupledBoundaryOrdering
(
oldOwner[oldIndex] < nOldCells_
? reverseCellMap_[oldOwner[oldIndex]]
: addedCellRenumbering_[oldOwner[oldIndex]]
);
centres,
anchors,
patchMaps,
rotations
)
);
// Insert entities into local listsLists...
faces_[faceInOrder] = oldFaces[oldIndex];
owner_[faceInOrder] = ownerRenumber;
neighbour_[faceInOrder] = -1;
faceEdges_[faceInOrder] = oldFaceEdges[oldIndex];
if (anyChange)
{
forAll(patchMaps, pI)
{
const labelList& patchMap = patchMaps[pI];
// Insert entities into mesh-reset lists
// NOTE: From OF-1.5 onwards, neighbour array
// does not store -1 for boundary faces
faces[faceInOrder].transfer(oldFaces[oldIndex]);
owner[faceInOrder] = ownerRenumber;
faceEdges[faceInOrder].transfer(oldFaceEdges[oldIndex]);
if (patchMap.size())
{
// Faces in this patch need to be shuffled
label pStart = patchStarts_[pI];
faceInOrder++;
forAll(patchMap, fI)
{
label oldPos = fI + pStart;
label newPos = patchMap[fI] + pStart;
// Rotate the face in-place
const face& oldFace = faces_[oldPos];
// Copy the old face, and rotate if necessary
face newFace(oldFace);
label nPos = rotations[pI][fI];
if (nPos != 0)
{
forAll(oldFace, fpI)
{
label nfpI = (fpI + nPos) % oldFace.size();
if (nfpI < 0)
{
nfpI += oldFace.size();
}
newFace[nfpI] = oldFace[fpI];
}
}
// Map to new locations, using the
// 'old' lists as temporaries.
// - All boundary neighbours are '-1',
// so use that as a faceMap temporary
oldNeighbour[newPos] = faceMap_[oldPos];
// Check if this face was added
if (faceMap_[oldPos] == -1)
{
// Fetch the index for the new face
label addedIndex = addedFaceReverseRenumbering[oldPos];
// Renumber for the added index
addedFaceRenumbering_[addedIndex] = newPos;
}
else
{
// Fetch the index for the existing face
label existIndex = faceMap_[oldPos];
// Renumber the reverse map
reverseFaceMap_[existIndex] = newPos;
}
oldOwner[newPos] = owner_[oldPos];
oldFaces[newPos].transfer(newFace);
oldFaceEdges[newPos].transfer(faceEdges_[oldPos]);
}
// Now copy / transfer back to original lists
forAll(patchMap, fI)
{
label index = (fI + pStart);
faceMap_[index] = oldNeighbour[index];
owner_[index] = oldOwner[index];
faces_[index] = oldFaces[index];
faceEdges_[index] = oldFaceEdges[index];
owner[index] = oldOwner[index];
faces[index].transfer(oldFaces[index]);
faceEdges[index].transfer(oldFaceEdges[index]);
}
}
}
}
// Now that we're done with faces, unlock it
@ -1104,9 +1226,54 @@ void dynamicTopoFvMesh::reOrderFaces
// Reset all zones
faceZones.updateMesh();
// Loop through all local coupling maps, and renumber faces.
forAll(patchCoupling_, patchI)
{
if (!patchCoupling_(patchI))
{
continue;
}
// Obtain references
const coupleMap& cMap = patchCoupling_[patchI].map();
const Map<label>& mtsMap = cMap.entityMap(coupleMap::FACE);
Map<label> newMtsMap, newStmMap;
forAllConstIter(Map<label>, mtsMap, mIter)
{
label newMaster = -1, newSlave = -1;
if (mIter.key() < nOldFaces_)
{
newMaster = reverseFaceMap_[mIter.key()];
}
else
{
newMaster = addedFaceRenumbering_[mIter.key()];
}
if (mIter() < nOldFaces_)
{
newSlave = reverseFaceMap_[mIter()];
}
else
{
newSlave = addedFaceRenumbering_[mIter()];
}
// Update the map.
newMtsMap.insert(newMaster, newSlave);
newStmMap.insert(newSlave, newMaster);
}
// Overwrite the old maps.
cMap.transferMaps(coupleMap::FACE, newMtsMap, newStmMap);
}
if (debug)
{
Info << "Done." << endl;
Info<< "Done." << endl;
}
}
@ -1184,10 +1351,10 @@ void dynamicTopoFvMesh::reOrderCells
{
if (threaded)
{
Info << "Thread: " << self() << ": ";
Info<< "Thread: " << self() << ": ";
}
Info << "ReOrdering cells..." << flush;
Info<< "ReOrdering cells..." << flush;
}
// Allocate for mapping information
@ -1507,7 +1674,7 @@ void dynamicTopoFvMesh::reOrderCells
if (debug)
{
Info << "Done." << endl;
Info<< "Done." << endl;
}
}
@ -1558,47 +1725,53 @@ void dynamicTopoFvMesh::reOrderMesh
{
if (debug)
{
Info << endl;
Info << "=================" << endl;
Info << " Mesh reOrdering " << endl;
Info << "=================" << endl;
Info << "Mesh Info [n]:" << endl;
Info << "Points: " << nOldPoints_ << endl;
Info << "Edges: " << nOldEdges_ << endl;
Info << "Faces: " << nOldFaces_ << endl;
Info << "Cells: " << nOldCells_ << endl;
Info << "Internal Edges: " << nOldInternalEdges_ << endl;
Info << "Internal Faces: " << nOldInternalFaces_ << endl;
Info<< nl
<< "=================" << nl
<< " Mesh reOrdering " << nl
<< "=================" << nl
<< " Mesh Info [n]:" << nl
<< " Points: " << nOldPoints_ << nl
<< " Edges: " << nOldEdges_ << nl
<< " Faces: " << nOldFaces_ << nl
<< " Cells: " << nOldCells_ << nl
<< " Internal Edges: " << nOldInternalEdges_ << nl
<< " Internal Faces: " << nOldInternalFaces_ << nl
<< " nPatches: " << nPatches_ << nl;
if (debug > 1)
{
Info << "Patch Starts [Face]: " << oldPatchStarts_ << endl;
Info << "Patch Sizes [Face]: " << oldPatchSizes_ << endl;
Info << "Patch Starts [Edge]: " << oldEdgePatchStarts_ << endl;
Info << "Patch Sizes [Edge]: " << oldEdgePatchSizes_ << endl;
Info<< " Patch Starts [Face]: " << oldPatchStarts_ << nl
<< " Patch Sizes [Face]: " << oldPatchSizes_ << nl
<< " Patch Starts [Edge]: " << oldEdgePatchStarts_ << nl
<< " Patch Sizes [Edge]: " << oldEdgePatchSizes_ << nl;
}
Info << "=================" << endl;
Info << "Mesh Info [n+1]:" << endl;
Info << "Points: " << nPoints_ << endl;
Info << "Edges: " << nEdges_ << endl;
Info << "Faces: " << nFaces_ << endl;
Info << "Cells: " << nCells_ << endl;
Info << "Internal Edges: " << nInternalEdges_ << endl;
Info << "Internal Faces: " << nInternalFaces_ << endl;
Info<< "=================" << nl
<< " Mesh Info [n+1]:" << nl
<< " Points: " << nPoints_ << nl
<< " Edges: " << nEdges_ << nl
<< " Faces: " << nFaces_ << nl
<< " Cells: " << nCells_ << nl
<< " Internal Edges: " << nInternalEdges_ << nl
<< " Internal Faces: " << nInternalFaces_ << nl
<< " nPatches: " << nPatches_ << nl;
if (debug > 1)
{
Info << "Patch Starts [Face]: " << patchStarts_ << endl;
Info << "Patch Sizes: [Face]: " << patchSizes_ << endl;
Info << "Patch Starts [Edge]: " << edgePatchStarts_ << endl;
Info << "Patch Sizes: [Edge]: " << edgePatchSizes_ << endl;
Info<< " Patch Starts [Face]: " << patchStarts_ << nl
<< " Patch Sizes: [Face]: " << patchSizes_ << nl
<< " Patch Starts [Edge]: " << edgePatchStarts_ << nl
<< " Patch Sizes: [Edge]: " << edgePatchSizes_ << nl;
}
Info << "=================" << endl;
Info<< "=================" << endl;
// Check connectivity structures for consistency
// before entering the reOrdering phase.
checkConnectivity();
}
if (threader_->multiThreaded())
if (threader_->multiThreaded() && !Pstream::parRun())
{
// Initialize multi-threaded reOrdering
threadedMeshReOrdering

27
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/eMesh/eMesh.C
View file

@ -42,7 +42,6 @@ defineTypeNameAndDebug(eMesh, 0);
word eMesh::meshSubDir = "eMesh";
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void eMesh::clearGeom() const
@ -65,8 +64,6 @@ void eMesh::clearAddressing() const
edges_.clear();
deleteDemandDrivenData(pePtr_);
deleteDemandDrivenData(epPtr_);
deleteDemandDrivenData(fePtr_);
deleteDemandDrivenData(efPtr_);
}
@ -149,8 +146,6 @@ eMesh::eMesh(const polyMesh& pMesh, const word& subDir)
),
*this
),
pePtr_(NULL),
epPtr_(NULL),
fePtr_(NULL),
efPtr_(NULL)
{
@ -393,28 +388,6 @@ void eMesh::setInstance(const fileName& inst)
}
const labelListList& eMesh::pointEdges() const
{
if (!pePtr_)
{
calcPointEdges();
}
return *pePtr_;
}
const labelListList& eMesh::edgePoints() const
{
if (!epPtr_)
{
calcEdgePoints();
}
return *epPtr_;
}
const labelListList& eMesh::faceEdges() const
{
if (!fePtr_)

19
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/eMesh/eMesh.H
View file

@ -80,12 +80,6 @@ class eMesh
//- Mapping between ordered edge-data and regular edges.
labelList reverseEdgeMap_;
//- Point-edges
mutable labelListList* pePtr_;
//- Edge-points
mutable labelListList* epPtr_;
//- Face-edges
mutable labelListIOList* fePtr_;
@ -111,12 +105,6 @@ class eMesh
//- Calculate ordered edges (and edgeFaces)
void calcOrderedEdgeList();
//- Calculate ordered edgePoints
void calcEdgePoints() const;
//- Calculate point-edges
void calcPointEdges() const;
//- Calculate face-edges
void calcFaceEdges() const;
@ -197,13 +185,6 @@ public:
// Demand-driven data
//- Return constant reference to the pointEdges list
const labelListList& pointEdges() const;
//- Return constant reference to the edgePoints list
//- Valid for 3D tetrahedral meshes only.
const labelListList& edgePoints() const;
//- Return constant reference to the faceEdges list
const labelListList& faceEdges() const;

191
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/eMesh/eMeshDemandDrivenData.C
View file

@ -24,6 +24,8 @@ License
Description
Demand-driven edge mesh data
\*---------------------------------------------------------------------------*/
#include "eMesh.H"
@ -144,195 +146,6 @@ void eMesh::calcOrderedEdgeList()
calcFaceEdges();
}
void eMesh::calcEdgePoints() const
{
if (debug)
{
Info<< "void eMesh::calcEdgePoints() const : "
<< "Calculating ordered edgePoints" << endl;
}
if (epPtr_)
{
FatalErrorIn
(
"void eMesh::calcEdgePoints() const"
) << "edgePoints already allocated."
<< abort(FatalError);
}
if (!efPtr_->size())
{
FatalErrorIn
(
"void eMesh::calcEdgePoints() const"
) << "edgeFaces doesn't exist."
<< abort(FatalError);
}
// Size the list
epPtr_ = new labelListList(nEdges_);
labelListList& edgePoints = *epPtr_;
// EdgePoints are ordered such that points appear in
// counter-clockwise direction around point[0] of the edge.
// If a boundary point is present, that must begin the list.
// NOTE: Will work only on tetrahedral meshes!
bool found;
label faceIndex = -1, cellIndex = -1;
const labelList& owner = mesh_.faceOwner();
const labelList& neighbour = mesh_.faceNeighbour();
const cellList& cells = mesh_.cells();
const faceList& faces = mesh_.faces();
const labelListList& eFaces = this->edgeFaces();
for (label eIndex=0; eIndex < nEdges_; eIndex++)
{
const edge& e = edges_[eIndex];
const labelList& eFace = eFaces[eIndex];
// Size the list
edgePoints[eIndex].setSize(eFace.size(), -1);
if (eIndex < nInternalEdges_)
{
// Pick the first face and start with that
faceIndex = eFace[0];
}
else
{
// Need to find a properly oriented start-face
forAll(eFace, faceI)
{
if
(
(!mesh_.isInternalFace(eFace[faceI]))
&& (edgeDirection(faces[eFace[faceI]], e) == 1)
)
{
faceIndex = eFace[faceI];
break;
}
}
}
// Shuffle vertices to appear in CCW order
for (label j=0; j < eFace.size(); j++)
{
const face& f = faces[faceIndex];
// Add the isolated point
edgePoints[eIndex][j] = findIsolatedPoint(f, e);
// Figure out how this edge is oriented.
if (edgeDirection(f, e) == 1)
{
// Counter-clockwise. Pick the owner.
cellIndex = owner[faceIndex];
}
else
if (mesh_.isInternalFace(faceIndex))
{
// Clockwise. Pick the neighbour.
cellIndex = neighbour[faceIndex];
}
else
{
// Looks like we've hit a boundary face. Break out.
break;
}
const cell& cellToCheck = cells[cellIndex];
found = false;
// Assuming tet-cells,
// Loop through edgeFaces and get the next face
forAll(eFace, faceI)
{
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[0]
)
{
faceIndex = cellToCheck[0];
found = true; break;
}
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[1]
)
{
faceIndex = cellToCheck[1];
found = true; break;
}
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[2]
)
{
faceIndex = cellToCheck[2];
found = true; break;
}
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[3]
)
{
faceIndex = cellToCheck[3];
found = true; break;
}
}
# ifdef FULLDEBUG
if (!found)
{
// Something's terribly wrong
FatalErrorIn
(
"void eMesh::calcEdgePoints() const"
)
<< " Failed to determine a vertex ring. " << nl
<< " edgeFaces connectivity is inconsistent. " << nl
<< "Edge: " << eIndex << ":: " << e << nl
<< "edgeFaces: " << eFace
<< abort(FatalError);
}
# endif
}
}
}
void eMesh::calcPointEdges() const
{
if (debug)
{
Info<< "void eMesh::calcPointEdges() const : "
<< "Calculating PointEdges" << endl;
}
if (pePtr_)
{
FatalErrorIn
(
"void eMesh::calcPointEdges() const"
) << "pePtr_ already allocated"
<< abort(FatalError);
}
pePtr_ = new labelListList(mesh_.nPoints());
invertManyToMany(mesh_.nPoints(), edges(), *pePtr_);
}
void eMesh::calcFaceEdges() const
{

10
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/eMesh/ePatches/ePatch/ePatch.C
View file

@ -80,16 +80,6 @@ ePatch::ePatch
size_(readLabel(dict.lookup("size")))
{}
ePatch::ePatch(const ePatch& p)
:
patchIdentifier(p, p.index()),
boundaryMesh_(p.boundaryMesh_),
start_(p.start()),
size_(p.size())
{}
ePatch::ePatch(const ePatch& p, const eBoundaryMesh& bm)
:
patchIdentifier(p, p.index()),

12
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/eMesh/ePatches/ePatch/ePatch.H
View file

@ -69,6 +69,13 @@ private:
//- Size of the patch
label size_;
// Demand-driven private data
// Private Member Functions
//- Disallow construct as copy
ePatch(const ePatch&);
protected:
@ -163,9 +170,6 @@ public:
const eBoundaryMesh& bm
);
//- Construct as copy
ePatch(const ePatch&);
//- Construct as copy, resetting the boundary mesh
ePatch(const ePatch&, const eBoundaryMesh&);
@ -185,7 +189,7 @@ public:
);
//- Return a pointer to a new patch created
//- Return a pointer to a new patch created
// on freestore from dictionary
static autoPtr<ePatch> New
(

4154
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeBisect.C
View file

File diff suppressed because it is too large Load diff

5479
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeCollapse.C
View file

File diff suppressed because it is too large Load diff

1695
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeSwap.C
View file

File diff suppressed because it is too large Load diff

150
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/conservativeMapFields.H Normal file
View file

@ -0,0 +1,150 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "topoMapper.H"
#include "topoCellMapper.H"
#include "topoSurfaceMapper.H"
#include "topoBoundaryMeshMapper.H"
namespace Foam
{
// Conservatively map all volFields in the registry
template <class Type>
void conservativeMapVolFields(const topoMapper& mapper)
{
// Define a few typedefs for convenience
typedef typename outerProduct<vector, Type>::type gCmptType;
typedef GeometricField<Type, fvPatchField, volMesh> volType;
typedef GeometricField<gCmptType, fvPatchField, volMesh> gradVolType;
HashTable<const volType*> fields(mapper.mesh().lookupClass<volType>());
// Store old-times before mapping
forAllIter(typename HashTable<const volType*>, fields, fIter)
{
fIter()->storeOldTimes();
}
// Fetch internal/boundary mappers
const topoCellMapper& fMap = mapper.volMap();
const topoBoundaryMeshMapper& bMap = mapper.boundaryMap();
// Now map all fields
forAllIter(typename HashTable<const volType*>, fields, fIter)
{
volType& field = const_cast<volType&>(*fIter());
if (fvMesh::debug)
{
Info<< "Conservatively mapping "
<< field.typeName
<< ' ' << field.name()
<< endl;
}
// Map the internal field
fMap.mapInternalField
(
field.name(),
mapper.gradient<gradVolType>(field.name()).internalField(),
field.internalField()
);
// Map patch fields
forAll(bMap, patchI)
{
bMap[patchI].mapPatchField
(
field.name(),
field.boundaryField()[patchI]
);
}
// Set the field instance
field.instance() = field.mesh().thisDb().time().timeName();
}
}
// Conservatively map all surfaceFields in the registry
template <class Type>
void conservativeMapSurfaceFields(const topoMapper& mapper)
{
// Define a few typedefs for convenience
typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfType;
HashTable<const surfType*> fields(mapper.mesh().lookupClass<surfType>());
// Store old-times before mapping
forAllIter(typename HashTable<const surfType*>, fields, fIter)
{
fIter()->storeOldTimes();
}
// Fetch internal/boundary mappers
const topoSurfaceMapper& fMap = mapper.surfaceMap();
const topoBoundaryMeshMapper& bMap = mapper.boundaryMap();
// Now map all fields
forAllIter(typename HashTable<const surfType*>, fields, fIter)
{
surfType& field = const_cast<surfType&>(*fIter());
if (fvMesh::debug)
{
Info<< "Conservatively mapping "
<< field.typeName
<< ' ' << field.name()
<< endl;
}
// Map the internal field
fMap.mapInternalField
(
field.name(),
field.internalField()
);
// Map patch fields
forAll(bMap, patchI)
{
bMap[patchI].mapPatchField
(
field.name(),
field.boundaryField()[patchI]
);
}
// Set the field instance
field.instance() = field.mesh().thisDb().time().timeName();
}
}
} // End namespace Foam
// ************************************************************************* //

2
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/fluxCorrector.C
View file

@ -33,7 +33,7 @@ Author
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
\*---------------------------------------------------------------------------*/
#include "fluxCorrector.H"
#include "dlLibraryTable.H"

6
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/fluxCorrector.H
View file

@ -49,8 +49,6 @@ SourceFiles
namespace Foam
{
// Forward class declarations
/*---------------------------------------------------------------------------*\
Class fluxCorrector Declaration
\*---------------------------------------------------------------------------*/
@ -142,10 +140,6 @@ public:
} // End namespace Foam
#ifdef NoRepository
#include "fluxCorrector.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

112
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoCellMapper.C
View file

@ -33,9 +33,8 @@ Author
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
\*---------------------------------------------------------------------------*/
#include "IOmanip.H"
#include "topoMapper.H"
#include "mapPolyMesh.H"
#include "topoCellMapper.H"
@ -306,6 +305,7 @@ const List<vectorField>& topoCellMapper::intersectionCentres() const
return *centresPtr_;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
@ -319,6 +319,7 @@ topoCellMapper::topoCellMapper
mpm_(mpm),
tMapper_(mapper),
direct_(false),
sizeBeforeMapping_(mpm.nOldCells()),
directAddrPtr_(NULL),
interpolationAddrPtr_(NULL),
weightsPtr_(NULL),
@ -326,6 +327,18 @@ topoCellMapper::topoCellMapper
volumesPtr_(NULL),
centresPtr_(NULL)
{
// Fetch offset sizes from topoMapper
const labelList& sizes = tMapper_.cellSizes();
// Add offset sizes
if (sizes.size())
{
forAll(sizes, pI)
{
sizeBeforeMapping_ += sizes[pI];
}
}
// Check for possibility of direct mapping
if
(
@ -363,7 +376,7 @@ label topoCellMapper::size() const
//- Return size before mapping
label topoCellMapper::sizeBeforeMapping() const
{
return mpm_.nOldCells();
return sizeBeforeMapping_;
}
@ -459,99 +472,6 @@ const labelList& topoCellMapper::insertedObjectLabels() const
}
//- Conservatively map the internal field
template <class Type, class gradType>
void topoCellMapper::mapInternalField
(
const word& fieldName,
const Field<gradType>& gF,
Field<Type>& iF
) const
{
if (iF.size() != sizeBeforeMapping() || gF.size() != sizeBeforeMapping())
{
FatalErrorIn
(
"\n\n"
"void topoCellMapper::mapInternalField<Type>\n"
"(\n"
" const word& fieldName,\n"
" const Field<gradType>& gF,\n"
" Field<Type>& iF\n"
") const\n"
) << "Incompatible size before mapping." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << iF.size() << nl
<< " Gradient Field size: " << gF.size() << nl
<< " map size: " << sizeBeforeMapping() << nl
<< abort(FatalError);
}
// Fetch addressing
const labelListList& cAddressing = addressing();
const List<scalarField>& wC = intersectionWeights();
const List<vectorField>& xC = intersectionCentres();
// Fetch geometry
const vectorField& centres = tMapper_.internalCentres();
// Compute the integral of the source field
Type intSource = sum(iF * tMapper_.cellVolumes());
// Copy the original field
Field<Type> fieldCpy(iF);
// Resize to current dimensions
iF.setSize(size());
// Map the internal field
forAll(iF, cellI)
{
const labelList& addr = cAddressing[cellI];
iF[cellI] = pTraits<Type>::zero;
forAll(addr, cellJ)
{
// Accumulate volume-weighted Taylor-series interpolate
iF[cellI] +=
(
wC[cellI][cellJ] *
(
fieldCpy[addr[cellJ]]
+ (
gF[addr[cellJ]] &
(
xC[cellI][cellJ]
- centres[addr[cellJ]]
)
)
)
);
}
}
// Compute the integral of the target field
Type intTarget = sum(iF * mesh_.cellVolumes());
if (polyMesh::debug)
{
int oldP = Info().precision();
// Compare the global integral
Info << " Field : " << fieldName
<< " integral errors : "
<< setprecision(15)
<< " source : " << mag(intSource)
<< " target : " << mag(intTarget)
<< " norm : "
<< (mag(intTarget - intSource) / (mag(intSource) + VSMALL))
<< setprecision(oldP)
<< endl;
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void topoCellMapper::operator=(const topoCellMapper& rhs)

13
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoCellMapper.H
View file

@ -43,6 +43,7 @@ SourceFiles
#ifndef topoCellMapper_H
#define topoCellMapper_H
#include "topoMapper.H"
#include "morphFieldMapper.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -50,10 +51,6 @@ SourceFiles
namespace Foam
{
// Class forward declarations
class topoMapper;
class mapPolyMesh;
/*---------------------------------------------------------------------------*\
Class topoCellMapper Declaration
\*---------------------------------------------------------------------------*/
@ -76,6 +73,9 @@ class topoCellMapper
//- Is the mapping direct
bool direct_;
//- Size before mapping
mutable label sizeBeforeMapping_;
// Demand-driven private data
//- Direct addressing
@ -178,9 +178,12 @@ public:
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "topoCellMapper.C"
# include "topoCellMapperTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

112
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoCellMapperTemplates.C Normal file
View file

@ -0,0 +1,112 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Conservatively map the internal field
template <class Type, class gradType>
void topoCellMapper::mapInternalField
(
const word& fieldName,
const Field<gradType>& gF,
Field<Type>& iF
) const
{
if (iF.size() != sizeBeforeMapping() || gF.size() != sizeBeforeMapping())
{
FatalErrorIn
(
"\n\n"
"void topoCellMapper::mapInternalField<Type>\n"
"(\n"
" const word& fieldName,\n"
" const Field<gradType>& gF,\n"
" Field<Type>& iF\n"
") const\n"
) << "Incompatible size before mapping." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << iF.size() << nl
<< " Gradient Field size: " << gF.size() << nl
<< " map size: " << sizeBeforeMapping() << nl
<< abort(FatalError);
}
// If we have direct addressing, map and bail out
if (direct())
{
iF.autoMap(*this);
return;
}
// Fetch addressing
const labelListList& cAddressing = addressing();
const List<scalarField>& wC = intersectionWeights();
const List<vectorField>& xC = intersectionCentres();
// Fetch geometry
const vectorField& centres = tMapper_.internalCentres();
// Copy the original field
Field<Type> fieldCpy(iF);
// Resize to current dimensions
iF.setSize(size());
// Map the internal field
forAll(iF, cellI)
{
const labelList& addr = cAddressing[cellI];
iF[cellI] = pTraits<Type>::zero;
forAll(addr, cellJ)
{
// Accumulate volume-weighted Taylor-series interpolate
iF[cellI] +=
(
wC[cellI][cellJ] *
(
fieldCpy[addr[cellJ]]
+ (
gF[addr[cellJ]] &
(
xC[cellI][cellJ]
- centres[addr[cellJ]]
)
)
)
);
}
}
}
} // End namespace Foam
// ************************************************************************* //

473
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoMapper.C
View file

@ -35,120 +35,18 @@ Author
\*----------------------------------------------------------------------------*/
#include "fvc.H"
#include "topoMapper.H"
#include "fluxCorrector.H"
#include "topoCellMapper.H"
#include "leastSquaresGrad.H"
#include "topoSurfaceMapper.H"
#include "topoBoundaryMeshMapper.H"
#include "fixedValueFvPatchFields.H"
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Store gradients of fields on the mesh prior to topology changes
template <class Type, class gradType>
void topoMapper::storeGradients
(
HashTable<autoPtr<gradType> >& gradTable
) const
{
// Define a few typedefs for convenience
typedef typename outerProduct<vector, Type>::type gCmptType;
typedef GeometricField<Type, fvPatchField, volMesh> volType;
typedef GeometricField<gCmptType, fvPatchField, volMesh> gVolType;
// Fetch all fields from registry
HashTable<const volType*> fields
(
mesh_.objectRegistry::lookupClass<volType>()
);
forAllConstIter(typename HashTable<const volType*>, fields, fIter)
{
const volType& field = *fIter();
// Compute the gradient.
tmp<gVolType> tGrad;
// If the fvSolution dictionary contains an entry,
// use that, otherwise, default to leastSquares
word gradName("grad(" + field.name() + ')');
if (mesh_.schemesDict().subDict("gradSchemes").found(gradName))
{
tGrad = fvc::grad(field);
}
else
{
tGrad = fv::leastSquaresGrad<Type>(mesh_).grad(field);
}
// Make a new entry, but don't register the field.
gradTable.insert
(
field.name(),
autoPtr<gradType>
(
new gradType
(
IOobject
(
tGrad().name(),
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
tGrad()
)
)
);
}
}
//- Fetch the gradient field (template specialisation)
template <>
const volVectorField& topoMapper::gradient(const word& name) const
{
if (!sGrads_.found(name))
{
FatalErrorIn
(
"const volVectorField& "
"topoMapper::gradient(const word& name) const"
) << nl << " Gradient for: " << name
<< " has not been stored."
<< abort(FatalError);
}
return sGrads_[name]();
}
//- Fetch the gradient field (template specialisation)
template <>
const volTensorField& topoMapper::gradient(const word& name) const
{
if (!vGrads_.found(name))
{
FatalErrorIn
(
"const volTensorField& "
"topoMapper::gradient(const word& name) const"
) << nl << " Gradient for: " << name
<< " has not been stored."
<< abort(FatalError);
}
return vGrads_[name]();
}
//- Store gradients prior to mesh reset
void topoMapper::storeGradients() const
{
@ -157,11 +55,8 @@ void topoMapper::storeGradients() const
if (fvMesh::debug)
{
Info << "Registered volScalarFields: " << endl;
Info << sGrads_.toc() << endl;
Info << "Registered volVectorFields: " << endl;
Info << vGrads_.toc() << endl;
Info<< "Registered volScalarFields: " << scalarGrads() << endl;
Info<< "Registered volVectorFields: " << vectorGrads() << endl;
}
}
@ -169,52 +64,78 @@ void topoMapper::storeGradients() const
//- Store geometric information
void topoMapper::storeGeometry() const
{
if (cellCentresPtr_)
// Wipe out existing information
deleteDemandDrivenData(cellCentresPtr_);
vectorField Cv(mesh_.cellCentres());
vectorField Cf(mesh_.faceCentres());
// Create and map the patch field values
label nPatches = mesh_.boundary().size();
// Create field parts
PtrList<fvPatchField<vector> > volCentrePatches(nPatches);
// Over-ride and set all patches to fixedValue
for (label patchI = 0; patchI < nPatches; patchI++)
{
deleteDemandDrivenData(cellVolumesPtr_);
deleteDemandDrivenData(cellCentresPtr_);
volCentrePatches.set
(
patchI,
new fixedValueFvPatchField<vector>
(
mesh_.boundary()[patchI],
DimensionedField<vector, volMesh>::null()
)
);
patchAreasPtr_.clear();
patchCentresPtr_.clear();
// Slice field to patch (forced assignment)
volCentrePatches[patchI] ==
(
mesh_.boundaryMesh()[patchI].patchSlice(Cf)
);
}
// Set the cell-volumes pointer.
cellVolumesPtr_ = new scalarField(mesh_.cellVolumes());
// Set the cell-centres pointer.
cellCentresPtr_ = new vectorField(mesh_.cellCentres());
// Set patch-areas
patchAreasPtr_.setSize(mesh_.boundaryMesh().size());
forAll(mesh_.boundaryMesh(), patchI)
{
patchAreasPtr_.set
cellCentresPtr_ =
(
new volVectorField
(
patchI,
new scalarField
IOobject
(
mag(mesh_.boundaryMesh()[patchI].faceAreas())
)
);
}
// Set patch-centres.
patchCentresPtr_.setSize(mesh_.boundaryMesh().size());
forAll(mesh_.boundaryMesh(), patchI)
{
patchCentresPtr_.set
(
patchI,
new vectorField
(
mesh_.boundaryMesh()[patchI].faceCentres()
)
);
}
"cellCentres",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimLength,
SubField<vector>(Cv, mesh_.nCells()),
volCentrePatches
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
//- Construct from mesh and dictionary
topoMapper::topoMapper
(
const fvMesh& mesh,
const dictionary& dict
)
:
mesh_(mesh),
dict_(dict),
cellMap_(NULL),
surfaceMap_(NULL),
boundaryMap_(NULL),
fluxCorrector_(fluxCorrector::New(mesh, dict)),
cellCentresPtr_(NULL)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
@ -287,6 +208,26 @@ void topoMapper::setCellWeights
}
//- Set cell / patch offset information
void topoMapper::setOffsets
(
const labelList& cellSizes,
const labelList& cellStarts,
const labelList& faceSizes,
const labelList& faceStarts,
const labelListList& patchSizes,
const labelListList& patchStarts
) const
{
cellSizes_ = cellSizes;
cellStarts_ = cellStarts;
faceSizes_ = faceSizes;
faceStarts_ = faceStarts;
patchSizes_ = patchSizes;
patchStarts_ = patchStarts;
}
//- Fetch face weights
const List<scalarField>& topoMapper::faceWeights() const
{
@ -315,6 +256,48 @@ const List<vectorField>& topoMapper::cellCentres() const
}
//- Fetch cell sizes
const labelList& topoMapper::cellSizes() const
{
return cellSizes_;
}
//- Fetch face sizes
const labelList& topoMapper::faceSizes() const
{
return faceSizes_;
}
//- Fetch patch sizes
const labelListList& topoMapper::patchSizes() const
{
return patchSizes_;
}
//- Fetch cell starts
const labelList& topoMapper::cellStarts() const
{
return cellStarts_;
}
//- Fetch face starts
const labelList& topoMapper::faceStarts() const
{
return faceStarts_;
}
//- Fetch patch starts
const labelListList& topoMapper::patchStarts() const
{
return patchStarts_;
}
//- Store mesh information for the mapping stage
void topoMapper::storeMeshInformation() const
{
@ -325,20 +308,19 @@ void topoMapper::storeMeshInformation() const
storeGeometry();
}
//- Return stored cell volume information
const scalarField& topoMapper::cellVolumes() const
//- Return non-const access to cell centres
volVectorField& topoMapper::volCentres() const
{
if (!cellVolumesPtr_)
if (!cellCentresPtr_)
{
FatalErrorIn
(
"const scalarField& topoMapper::cellVolumes()"
"const vectorField& topoMapper::volCentres() const"
) << nl << " Pointer has not been set. "
<< abort(FatalError);
}
return *cellVolumesPtr_;
return *cellCentresPtr_;
}
@ -349,7 +331,7 @@ const vectorField& topoMapper::internalCentres() const
{
FatalErrorIn
(
"const vectorField& topoMapper::internalCentres()"
"const vectorField& topoMapper::internalCentres() const"
) << nl << " Pointer has not been set. "
<< abort(FatalError);
}
@ -358,154 +340,70 @@ const vectorField& topoMapper::internalCentres() const
}
//- Return stored patch areas information
const scalarField& topoMapper::patchAreas(const label i) const
{
if (!patchAreasPtr_.set(i))
{
FatalErrorIn
(
"const scalarField& topoMapper::patchAreas"
"(const label i) const"
) << nl << " Pointer has not been set at index: " << i
<< abort(FatalError);
}
return patchAreasPtr_[i];
}
//- Return stored patch centre information
const vectorField& topoMapper::patchCentres(const label i) const
{
if (!patchCentresPtr_.set(i))
if (!cellCentresPtr_)
{
FatalErrorIn
(
"const vectorField& topoMapper::patchCentres"
"(const label i) const"
) << nl << " Pointer has not been set at index: " << i
) << nl << " Pointer has not been set. index: " << i
<< abort(FatalError);
}
return patchCentresPtr_[i];
return (*cellCentresPtr_).boundaryField()[i];
}
// Conservatively map all volFields in the registry
template <class Type>
void topoMapper::conservativeMapVolFields() const
//- Return names of stored scalar gradients
const wordList topoMapper::scalarGrads() const
{
// Define a few typedefs for convenience
typedef typename outerProduct<vector, Type>::type gCmptType;
typedef GeometricField<Type, fvPatchField, volMesh> volType;
typedef GeometricField<gCmptType, fvPatchField, volMesh> gradVolType;
HashTable<const volType*> fields(mesh_.lookupClass<volType>());
// Store old-times before mapping
forAllIter(typename HashTable<const volType*>, fields, fIter)
{
volType& field = const_cast<volType&>(*fIter());
field.storeOldTimes();
}
// Fetch internal/boundary mappers
const topoCellMapper& fMap = volMap();
const topoBoundaryMeshMapper& bMap = boundaryMap();
// Now map all fields
forAllIter(typename HashTable<const volType*>, fields, fIter)
{
volType& field = const_cast<volType&>(*fIter());
if (fvMesh::debug)
{
Info << "Conservatively mapping "
<< field.typeName
<< ' ' << field.name()
<< endl;
}
// Map the internal field
fMap.mapInternalField
(
field.name(),
gradient<gradVolType>(field.name()).internalField(),
field.internalField()
);
// Map patch fields
forAll(bMap, patchI)
{
bMap[patchI].mapPatchField
(
field.name(),
field.boundaryField()[patchI]
);
}
// Set the field instance
field.instance() = field.mesh().thisDb().time().timeName();
}
return sGrads_.toc();
}
// Conservatively map all surfaceFields in the registry
template <class Type>
void topoMapper::conservativeMapSurfaceFields() const
//- Return names of stored vector gradients
const wordList topoMapper::vectorGrads() const
{
// Define a few typedefs for convenience
typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfType;
return vGrads_.toc();
}
HashTable<const surfType*> fields(mesh_.lookupClass<surfType>());
// Store old-times before mapping
forAllIter(typename HashTable<const surfType*>, fields, fIter)
//- Fetch the gradient field (template specialisation)
template <>
volVectorField& topoMapper::gradient(const word& name) const
{
if (!sGrads_.found(name))
{
surfType& field = const_cast<surfType&>(*fIter());
field.storeOldTimes();
}
// Fetch internal/boundary mappers
const topoSurfaceMapper& fMap = surfaceMap();
const topoBoundaryMeshMapper& bMap = boundaryMap();
// Now map all fields
forAllIter(typename HashTable<const surfType*>, fields, fIter)
{
surfType& field = const_cast<surfType&>(*fIter());
if (fvMesh::debug)
{
Info << "Conservatively mapping "
<< field.typeName
<< ' ' << field.name()
<< endl;
}
// Map the internal field
fMap.mapInternalField
FatalErrorIn
(
field.name(),
field.internalField()
);
// Map patch fields
forAll(bMap, patchI)
{
bMap[patchI].mapPatchField
(
field.name(),
field.boundaryField()[patchI]
);
}
// Set the field instance
field.instance() = field.mesh().thisDb().time().timeName();
"volVectorField& topoMapper::gradient(const word& name) const"
) << nl << " Gradient for: " << name
<< " has not been stored."
<< abort(FatalError);
}
return sGrads_[name]();
}
//- Fetch the gradient field (template specialisation)
template <>
volTensorField& topoMapper::gradient(const word& name) const
{
if (!vGrads_.found(name))
{
FatalErrorIn
(
"volTensorField& topoMapper::gradient(const word& name) const"
) << nl << " Gradient for: " << name
<< " has not been stored."
<< abort(FatalError);
}
return vGrads_[name]();
}
@ -533,7 +431,7 @@ const topoCellMapper& topoMapper::volMap() const
{
FatalErrorIn
(
"const topoCellMapper& topoMapper::volMap()"
"const topoCellMapper& topoMapper::volMap() const"
) << nl << " Volume mapper has not been set. "
<< abort(FatalError);
}
@ -549,7 +447,7 @@ const topoSurfaceMapper& topoMapper::surfaceMap() const
{
FatalErrorIn
(
"const topoSurfaceMapper& topoMapper::surfaceMap()"
"const topoSurfaceMapper& topoMapper::surfaceMap() const"
) << nl << " Surface mapper has not been set. "
<< abort(FatalError);
}
@ -565,7 +463,7 @@ const topoBoundaryMeshMapper& topoMapper::boundaryMap() const
{
FatalErrorIn
(
"const topoBoundaryMeshMapper& topoMapper::boundaryMap()"
"const topoBoundaryMeshMapper& topoMapper::boundaryMap() const"
) << nl << " Boundary mapper has not been set. "
<< abort(FatalError);
}
@ -581,7 +479,7 @@ const fluxCorrector& topoMapper::surfaceFluxCorrector() const
{
FatalErrorIn
(
"const fluxCorrector& topoMapper::surfaceFluxCorrector()"
"const fluxCorrector& topoMapper::surfaceFluxCorrector() const"
) << nl << " fluxCorrector has not been set. "
<< abort(FatalError);
}
@ -591,7 +489,7 @@ const fluxCorrector& topoMapper::surfaceFluxCorrector() const
//- Clear out member data
void topoMapper::clear()
void topoMapper::clear() const
{
// Clear out mappers
cellMap_.clear();
@ -603,10 +501,7 @@ void topoMapper::clear()
vGrads_.clear();
// Wipe out geomtry information
deleteDemandDrivenData(cellVolumesPtr_);
deleteDemandDrivenData(cellCentresPtr_);
patchAreasPtr_.clear();
patchCentresPtr_.clear();
// Clear maps
faceWeights_.clear();
@ -614,6 +509,16 @@ void topoMapper::clear()
faceCentres_.clear();
cellCentres_.clear();
// Clear sizes / offsets
cellSizes_.clear();
cellStarts_.clear();
faceSizes_.clear();
faceStarts_.clear();
patchSizes_.clear();
patchStarts_.clear();
}

95
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoMapper.H
View file

@ -42,7 +42,9 @@ SourceFiles
#ifndef topoMapper_H
#define topoMapper_H
#include "fluxCorrector.H"
#include "autoPtr.H"
#include "IOmanip.H"
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -50,10 +52,10 @@ namespace Foam
{
// Class forward declarations
class fluxCorrector;
class topoCellMapper;
class topoSurfaceMapper;
class topoBoundaryMeshMapper;
class mapPolyMesh;
/*---------------------------------------------------------------------------*\
Class topoMapper Declaration
@ -88,19 +90,24 @@ class topoMapper
mutable HashTable<autoPtr<volTensorField> > vGrads_;
//- Geometric information on the old mesh
mutable scalarField* cellVolumesPtr_;
mutable vectorField* cellCentresPtr_;
mutable PtrList<scalarField> patchAreasPtr_;
mutable PtrList<vectorField> patchCentresPtr_;
mutable volVectorField* cellCentresPtr_;
//- Intersection volume weights
//- Intersection weights
mutable List<scalarField> faceWeights_;
mutable List<scalarField> cellWeights_;
//- Intersection centre weights
//- Intersection centres
mutable List<vectorField> faceCentres_;
mutable List<vectorField> cellCentres_;
//- Sizes / starts for mapping
mutable labelList cellSizes_;
mutable labelList cellStarts_;
mutable labelList faceSizes_;
mutable labelList faceStarts_;
mutable labelListList patchSizes_;
mutable labelListList patchStarts_;
// Private Member Functions
//- Disallow default bitwise copy construct
@ -117,10 +124,6 @@ class topoMapper
HashTable<autoPtr<gradType> >& gradTable
) const;
//- Fetch the gradient field
template <class Type>
const Type& gradient(const word& name) const;
//- Store gradients prior to mesh reset
void storeGradients() const;
@ -132,19 +135,7 @@ public:
// Constructors
//- Construct from mesh and dictionary
topoMapper(const fvMesh& mesh, const dictionary& dict)
:
mesh_(mesh),
dict_(dict),
cellMap_(NULL),
surfaceMap_(NULL),
boundaryMap_(NULL),
fluxCorrector_(fluxCorrector::New(mesh, dict)),
cellVolumesPtr_(NULL),
cellCentresPtr_(NULL),
patchAreasPtr_(mesh.boundary().size()),
patchCentresPtr_(mesh.boundary().size())
{}
topoMapper(const fvMesh& mesh, const dictionary& dict);
// Destructor
@ -175,6 +166,17 @@ public:
const Xfer<List<vectorField> >& centres
) const;
//- Set cell / patch offset information
void setOffsets
(
const labelList& cellSizes,
const labelList& cellStarts,
const labelList& faceSizes,
const labelList& faceStarts,
const labelListList& patchSizes,
const labelListList& patchStarts
) const;
//- Fetch face weights
const List<scalarField>& faceWeights() const;
@ -187,28 +189,45 @@ public:
//- Fetch cell centres
const List<vectorField>& cellCentres() const;
//- Fetch cell sizes
const labelList& cellSizes() const;
//- Fetch face sizes
const labelList& faceSizes() const;
//- Fetch patch sizes
const labelListList& patchSizes() const;
//- Fetch cell starts
const labelList& cellStarts() const;
//- Fetch face starts
const labelList& faceStarts() const;
//- Fetch patch starts
const labelListList& patchStarts() const;
//- Store mesh information for the mapping stage
void storeMeshInformation() const;
//- Return stored cell volume information
const scalarField& cellVolumes() const;
//- Return non-const access to cell centres
volVectorField& volCentres() const;
//- Return stored cell centre information
const vectorField& internalCentres() const;
//- Return stored patch area information
const scalarField& patchAreas(const label i) const;
//- Return stored patch centre information
const vectorField& patchCentres(const label i) const;
// Conservatively map all volFields in the registry
template <class Type>
void conservativeMapVolFields() const;
//- Return names of stored scalar gradients
const wordList scalarGrads() const;
// Conservatively map all surfaceFields in the registry
//- Return names of stored vector gradients
const wordList vectorGrads() const;
//- Fetch the gradient field
template <class Type>
void conservativeMapSurfaceFields() const;
Type& gradient(const word& name) const;
//- Correct fluxes after topology change
void correctFluxes() const;
@ -226,7 +245,7 @@ public:
const fluxCorrector& surfaceFluxCorrector() const;
//- Clear out member data
void clear();
void clear() const;
};
@ -234,8 +253,10 @@ public:
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "topoMapper.C"
# include "topoMapperTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

106
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoMapperTemplates.C Normal file
View file

@ -0,0 +1,106 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "fvc.H"
#include "leastSquaresGrad.H"
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Store gradients of fields on the mesh prior to topology changes
template <class Type, class gradType>
void topoMapper::storeGradients
(
HashTable<autoPtr<gradType> >& gradTable
) const
{
// Define a few typedefs for convenience
typedef typename outerProduct<vector, Type>::type gCmptType;
typedef GeometricField<Type, fvPatchField, volMesh> volType;
typedef GeometricField<gCmptType, fvPatchField, volMesh> gVolType;
// Fetch all fields from registry
HashTable<const volType*> fields
(
mesh_.objectRegistry::lookupClass<volType>()
);
// Store old-times before gradient computation
forAllIter(typename HashTable<const volType*>, fields, fIter)
{
fIter()->storeOldTimes();
}
forAllConstIter(typename HashTable<const volType*>, fields, fIter)
{
const volType& field = *fIter();
// Compute the gradient.
tmp<gVolType> tGrad;
// If the fvSolution dictionary contains an entry,
// use that, otherwise, default to leastSquares
word gradName("grad(" + field.name() + ')');
if (mesh_.schemesDict().subDict("gradSchemes").found(gradName))
{
tGrad = fvc::grad(field);
}
else
{
tGrad = fv::leastSquaresGrad<Type>(mesh_).grad(field);
}
// Make a new entry, but don't register the field.
gradTable.insert
(
field.name(),
autoPtr<gradType>
(
new gradType
(
IOobject
(
tGrad().name(),
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
tGrad()
)
)
);
}
}
} // End namespace Foam
// ************************************************************************* //

301
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoPatchMapper.C
View file

@ -33,9 +33,10 @@ Author
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
\*---------------------------------------------------------------------------*/
#include "IOmanip.H"
#include "meshOps.H"
#include "topoMapper.H"
#include "mapPolyMesh.H"
#include "topoPatchMapper.H"
@ -73,9 +74,17 @@ void topoPatchMapper::calcInsertedFaceAddressing() const
}
// Information from the old patch
const label oldPatchSize = mpm_.oldPatchSizes()[patch_.index()];
const label oldPatchSize = sizeBeforeMapping();
const label oldPatchStart = mpm_.oldPatchStarts()[patch_.index()];
const label oldPatchEnd = oldPatchStart + oldPatchSize;
if (oldPatchSize == 0)
{
// Nothing to map from. Return empty.
insertedFaceLabelsPtr_ = new labelList(0);
insertedFaceIndexMapPtr_ = new labelList(0);
insertedFaceAddressingPtr_ = new labelListList(0);
return;
}
// Allocate for inserted face labels and addressing
label nInsertedFaces = 0;
@ -108,6 +117,19 @@ void topoPatchMapper::calcInsertedFaceAddressing() const
{
if (fffI.masterObjects().empty())
{
// Write out for post-processing
meshOps::writeVTK
(
mpm_.mesh(),
"patchFaceInsError_"
+ Foam::name(fffI.index()),
labelList(1, fffI.index()),
2,
mpm_.mesh().points(),
List<edge>(0),
mpm_.mesh().faces()
);
FatalErrorIn
(
"void topoPatchMapper::"
@ -130,23 +152,26 @@ void topoPatchMapper::calcInsertedFaceAddressing() const
// Make an entry for addressing
labelList& addr = insertedAddressing[nInsertedFaces];
// Renumber addressing to patch.
// Also, check mapping for hits into
// other patches / internal faces.
// Check for illegal addressing
addr = fffI.masterObjects();
forAll(addr, faceI)
{
if
(
addr[faceI] >= oldPatchStart
&& addr[faceI] < oldPatchEnd
)
{
addr[faceI] -= oldPatchStart;
}
else
if (addr[faceI] < 0 || addr[faceI] >= oldPatchSize)
{
// Write out for post-processing
meshOps::writeVTK
(
mpm_.mesh(),
"patchFacePatchError_"
+ Foam::name(fffI.index()),
labelList(1, fffI.index()),
2,
mpm_.mesh().points(),
List<edge>(0),
mpm_.mesh().faces()
);
FatalErrorIn
(
"void topoPatchMapper::"
@ -154,10 +179,13 @@ void topoPatchMapper::calcInsertedFaceAddressing() const
)
<< "Addressing into another patch is not allowed."
<< nl << " Patch face index: " << faceI
<< nl << " Patch: " << patch_.name()
<< nl << " fffI.index: " << fffI.index()
<< nl << " pStart: " << pStart
<< nl << " addr: " << addr
<< nl << " addr[faceI]: " << addr[faceI]
<< nl << " oldPatchStart: " << oldPatchStart
<< nl << " oldPatchSize: " << oldPatchSize
<< nl << " oldPatchEnd: " << oldPatchEnd
<< abort(FatalError);
}
}
@ -185,21 +213,31 @@ void topoPatchMapper::calcAddressing() const
}
// Information from the old patch
const label oldPatchSize = mpm_.oldPatchSizes()[patch_.index()];
const label oldPatchSize = sizeBeforeMapping();
const label oldPatchStart = mpm_.oldPatchStarts()[patch_.index()];
const label oldPatchEnd = oldPatchStart + oldPatchSize;
// Assemble the maps: slice to patch
if (direct())
{
if (oldPatchSize == 0)
{
// Nothing to map from. Return empty.
directAddrPtr_ = new labelList(0);
return;
}
// Direct mapping - slice to size
directAddrPtr_ = new labelList(patch_.patch().patchSlice(mpm_.faceMap()));
directAddrPtr_ =
(
new labelList(patch_.patch().patchSlice(mpm_.faceMap()))
);
labelList& addr = *directAddrPtr_;
// Shift to local patch indices.
// Also, check mapping for hits into other patches / internal faces.
forAll (addr, faceI)
forAll(addr, faceI)
{
if
(
@ -211,22 +249,55 @@ void topoPatchMapper::calcAddressing() const
}
else
{
// Relax addressing requirement for
// processor patch faces. These require
// cell-to-face interpolation anyway.
if (isA<processorPolyPatch>(patch_.patch()))
{
// Artificially map from face[0] of this patch.
addr[faceI] = 0;
continue;
}
// Write out for post-processing
meshOps::writeVTK
(
mpm_.mesh(),
"patchFacePatchError_"
+ Foam::name(faceI),
labelList(1, faceI),
2,
mpm_.mesh().points(),
List<edge>(0),
mpm_.mesh().faces()
);
FatalErrorIn
(
"void topoPatchMapper::calcAddressing() const"
)
<< "Addressing into another patch is not allowed."
<< nl << " Patch: " << patch_.name()
<< nl << " Patch index: " << patch_.index()
<< nl << " Patch face index: " << faceI
<< nl << " addr[faceI]: " << addr[faceI]
<< nl << " oldPatchStart: " << oldPatchStart
<< nl << " oldPatchSize: " << oldPatchSize
<< nl << " oldPatchEnd: " << oldPatchEnd
<< nl << " nInserted: " << insertedObjectLabels().size()
<< abort(FatalError);
}
}
}
else
{
if (oldPatchSize == 0)
{
// Nothing to map from. Return empty.
interpolationAddrPtr_ = new labelListList(0);
return;
}
// Interpolative addressing
interpolationAddrPtr_ = new labelListList(patchSize(), labelList(0));
labelListList& addr = *interpolationAddrPtr_;
@ -267,6 +338,8 @@ void topoPatchMapper::calcAddressing() const
"void topoPatchMapper::calcAddressing() const"
)
<< "Addressing into another patch is not allowed."
<< nl << " Patch: " << patch_.name()
<< nl << " Patch index: " << patch_.index()
<< nl << " Patch face index: " << faceI
<< nl << " faceMap[faceI]: " << oldFace
<< nl << " oldPatchStart: " << oldPatchStart
@ -284,6 +357,17 @@ void topoPatchMapper::calcAddressing() const
{
if (addr[faceI].empty())
{
// Relax addressing requirement for
// processor patch faces. These require
// cell-to-face interpolation anyway.
if (isA<processorPolyPatch>(patch_.patch()))
{
// Artificially map from face[0] of this patch.
addr[faceI] = labelList(1, 0);
continue;
}
FatalErrorIn
(
"void topoPatchMapper::calcAddressing() const"
@ -316,6 +400,13 @@ void topoPatchMapper::calcInverseDistanceWeights() const
// Fetch interpolative addressing
const labelListList& addr = addressing();
if (addr.empty())
{
// Nothing to map from. Return empty.
weightsPtr_ = new scalarListList(0);
return;
}
// Allocate memory
weightsPtr_ = new scalarListList(patchSize());
scalarListList& w = *weightsPtr_;
@ -386,12 +477,20 @@ void topoPatchMapper::calcIntersectionWeightsAndCentres() const
// Fetch interpolative addressing
const labelListList& addr = addressing();
// Allocate memory
areasPtr_ = new List<scalarField>(patchSize(), scalarField(0));
List<scalarField>& a = *areasPtr_;
if (addr.empty())
{
// Nothing to map from. Return empty.
areasPtr_ = new List<scalarList>(0);
centresPtr_ = new List<vectorList>(0);
return;
}
centresPtr_ = new List<vectorField>(patchSize(), vectorField(0));
List<vectorField>& x = *centresPtr_;
// Allocate memory
areasPtr_ = new List<scalarList>(patchSize(), scalarList(0));
List<scalarList>& a = *areasPtr_;
centresPtr_ = new List<vectorList>(patchSize(), vectorList(0));
List<vectorList>& x = *centresPtr_;
// Obtain stored face-centres
const vectorField& faceCentres = tMapper_.patchCentres(patch_.index());
@ -418,8 +517,8 @@ void topoPatchMapper::calcIntersectionWeightsAndCentres() const
// Check if this is indeed a mapped face
if (mo.size() == 1 && x[faceI].empty() && a[faceI].empty())
{
x[faceI] = vectorField(1, faceCentres[mo[0]]);
a[faceI] = scalarField(1, 1.0);
x[faceI] = vectorList(1, faceCentres[mo[0]]);
a[faceI] = scalarList(1, 1.0);
}
}
@ -441,13 +540,13 @@ void topoPatchMapper::calcIntersectionWeightsAndCentres() const
}
const List<scalarField>& topoPatchMapper::intersectionWeights() const
const List<scalarList>& topoPatchMapper::intersectionWeights() const
{
if (direct())
{
FatalErrorIn
(
"const List<scalarField>& "
"const List<scalarList>& "
"topoPatchMapper::intersectionWeights() const"
) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError);
@ -462,13 +561,13 @@ const List<scalarField>& topoPatchMapper::intersectionWeights() const
}
const List<vectorField>& topoPatchMapper::intersectionCentres() const
const List<vectorList>& topoPatchMapper::intersectionCentres() const
{
if (direct())
{
FatalErrorIn
(
"const List<vectorField>& "
"const List<vectorList>& "
"topoPatchMapper::intersectionCentres() const"
) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError);
@ -497,6 +596,8 @@ topoPatchMapper::topoPatchMapper
mpm_(mpm),
tMapper_(mapper),
direct_(false),
sizeBeforeMapping_(0),
conservative_(false),
directAddrPtr_(NULL),
interpolationAddrPtr_(NULL),
weightsPtr_(NULL),
@ -506,6 +607,39 @@ topoPatchMapper::topoPatchMapper
areasPtr_(NULL),
centresPtr_(NULL)
{
// Compute sizeBeforeMapping.
// - This needs to be done before insertedObjects
// is computed to determine direct mapping
if (isA<emptyPolyPatch>(patch_.patch()))
{
sizeBeforeMapping_ = 0;
}
else
{
label patchIndex = patch_.index();
label totalSize = mpm_.oldPatchSizes()[patchIndex];
// Fetch offset sizes from topoMapper
const labelListList& sizes = tMapper_.patchSizes();
// Add offset sizes
if (sizes.size())
{
// Fetch number of physical patches
label nPhysical = sizes[0].size();
if (patchIndex < nPhysical)
{
forAll(sizes, pI)
{
totalSize += sizes[pI][patchIndex];
}
}
}
sizeBeforeMapping_ = totalSize;
}
// Check for the possibility of direct mapping
if (insertedObjects())
{
@ -543,12 +677,7 @@ label topoPatchMapper::size() const
//- Return size before mapping
label topoPatchMapper::sizeBeforeMapping() const
{
if (patch_.type() == "empty")
{
return 0;
}
return mpm_.oldPatchSizes()[patch_.index()];
return sizeBeforeMapping_;
}
@ -616,6 +745,16 @@ const scalarListList& topoPatchMapper::weights() const
<< abort(FatalError);
}
if (conservative_)
{
if (!areasPtr_ && !centresPtr_)
{
calcIntersectionWeightsAndCentres();
}
return *areasPtr_;
}
if (!weightsPtr_)
{
calcInverseDistanceWeights();
@ -668,96 +807,6 @@ const labelListList& topoPatchMapper::insertedFaceAddressing() const
}
//- Map the patch field
template <class Type>
void topoPatchMapper::mapPatchField
(
const word& fieldName,
Field<Type>& pF
) const
{
// To invoke inverse-distance weighting, use this:
// pF.autoMap(*this);
// Check for possibility of direct mapping
if (direct())
{
pF.autoMap(*this);
}
else
{
if (pF.size() != sizeBeforeMapping())
{
FatalErrorIn
(
"\n\n"
"void topoCellMapper::mapPatchField<Type>\n"
"(\n"
" const word& fieldName,\n"
" Field<Type>& iF\n"
") const\n"
) << "Incompatible size before mapping." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << pF.size() << nl
<< " map size: " << sizeBeforeMapping() << nl
<< abort(FatalError);
}
// Fetch addressing
const labelListList& pAddressing = addressing();
const List<scalarField>& wF = intersectionWeights();
// Compute the integral of the source field
Type intSource = sum(pF * tMapper_.patchAreas(patch_.index()));
// Copy the original field
Field<Type> fieldCpy(pF);
// Resize to current dimensions
pF.setSize(size());
// Map the patch field
forAll(pF, faceI)
{
const labelList& addr = pAddressing[faceI];
pF[faceI] = pTraits<Type>::zero;
// Accumulate area-weighted interpolate
forAll(addr, faceJ)
{
pF[faceI] +=
(
wF[faceI][faceJ] * fieldCpy[addr[faceJ]]
);
}
}
// Compute the integral of the target field
const polyPatch& ppI = mpm_.mesh().boundaryMesh()[patch_.index()];
Type intTarget = sum(pF * mag(ppI.faceAreas()));
if (polyMesh::debug)
{
int oldP = Info().precision();
// Compare the global integral
Info << " Field : " << fieldName
<< " Patch : " << ppI.name()
<< " integral errors : "
<< setprecision(15)
<< " source : " << mag(intSource)
<< " target : " << mag(intTarget)
<< " norm : "
<< (mag(intTarget - intSource) / (mag(intSource) + VSMALL))
<< setprecision(oldP)
<< endl;
}
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void topoPatchMapper::operator=(const topoPatchMapper& rhs)

26
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoPatchMapper.H
View file

@ -43,17 +43,15 @@ SourceFiles
#ifndef topoPatchMapper_H
#define topoPatchMapper_H
#include "morphFieldMapper.H"
#include "vectorList.H"
#include "topoMapper.H"
#include "fvPatchFieldMapper.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
class topoMapper;
class mapPolyMesh;
/*---------------------------------------------------------------------------*\
Class topoPatchMapper Declaration
\*---------------------------------------------------------------------------*/
@ -76,6 +74,12 @@ class topoPatchMapper
//- Is the mapping direct
bool direct_;
//- Size before mapping
mutable label sizeBeforeMapping_;
//- Is the mapping conservative
mutable bool conservative_;
// Demand-driven private data
//- Direct addressing
@ -97,10 +101,10 @@ class topoPatchMapper
mutable labelListList* insertedFaceAddressingPtr_;
//- Interpolation areas
mutable List<scalarField>* areasPtr_;
mutable List<scalarList>* areasPtr_;
//- Interpolation centres
mutable List<vectorField>* centresPtr_;
mutable List<vectorList>* centresPtr_;
// Private Member Functions
@ -123,10 +127,10 @@ class topoPatchMapper
void calcIntersectionWeightsAndCentres() const;
//- Return intersection area weights
const List<scalarField>& intersectionWeights() const;
const List<scalarList>& intersectionWeights() const;
//- Return intersection area centres
const List<vectorField>& intersectionCentres() const;
const List<vectorList>& intersectionCentres() const;
//- Clear out local storage
void clearOut();
@ -196,8 +200,10 @@ public:
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "topoPatchMapper.C"
# include "topoPatchMapperTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

52
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoPatchMapperTemplates.C Normal file
View file

@ -0,0 +1,52 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "mapPolyMesh.H"
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Map the patch field
template <class Type>
void topoPatchMapper::mapPatchField
(
const word& fieldName,
Field<Type>& pF
) const
{
// Specify that mapping is conservative
conservative_ = true;
// Map patchField onto itself
pF.autoMap(*this);
}
} // End namespace Foam
// ************************************************************************* //

72
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoSurfaceMapper.C
View file

@ -33,7 +33,7 @@ Author
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
\*---------------------------------------------------------------------------*/
#include "topoMapper.H"
#include "mapPolyMesh.H"
@ -157,11 +157,24 @@ topoSurfaceMapper::topoSurfaceMapper
mpm_(mpm),
tMapper_(mapper),
direct_(false),
sizeBeforeMapping_(mpm.nOldInternalFaces()),
directAddrPtr_(NULL),
interpolationAddrPtr_(NULL),
weightsPtr_(NULL),
insertedFaceLabelsPtr_(NULL)
{
// Fetch offset sizes from topoMapper
const labelList& sizes = tMapper_.faceSizes();
// Add offset sizes
if (sizes.size())
{
forAll(sizes, pI)
{
sizeBeforeMapping_ += sizes[pI];
}
}
// Calculate the insertedFaces list
calcInsertedFaceLabels();
@ -188,7 +201,7 @@ label topoSurfaceMapper::size() const
//- Return size before mapping
label topoSurfaceMapper::sizeBeforeMapping() const
{
return mpm_.nOldInternalFaces();
return sizeBeforeMapping_;
}
@ -291,61 +304,6 @@ const labelHashSet& topoSurfaceMapper::flipFaceFlux() const
}
//- Map the internal field
template <class Type>
void topoSurfaceMapper::mapInternalField
(
const word& fieldName,
Field<Type>& iF
) const
{
if (iF.size() != sizeBeforeMapping())
{
FatalErrorIn
(
"\n\n"
"void topoSurfaceMapper::mapInternalField<Type>\n"
"(\n"
" Field<Type>& iF\n"
") const\n"
) << "Incompatible size before mapping." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << iF.size() << nl
<< " map size: " << sizeBeforeMapping() << nl
<< abort(FatalError);
}
// Map the internal field
iF.autoMap(*this);
// Flip the flux
const labelList flipFaces = flipFaceFlux().toc();
forAll (flipFaces, i)
{
if (flipFaces[i] < iF.size())
{
iF[flipFaces[i]] *= -1.0;
}
else
{
FatalErrorIn
(
"\n\n"
"void topoSurfaceMapper::mapInternalField<Type>\n"
"(\n"
" Field<Type>& iF\n"
") const\n"
) << "Cannot flip boundary face fluxes." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << iF.size() << nl
<< " Face flip index: " << flipFaces[i] << nl
<< abort(FatalError);
}
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void topoSurfaceMapper::operator=(const topoSurfaceMapper& rhs)

12
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoSurfaceMapper.H
View file

@ -43,6 +43,7 @@ SourceFiles
#ifndef topoSurfaceMapper_H
#define topoSurfaceMapper_H
#include "topoMapper.H"
#include "morphFieldMapper.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -50,10 +51,6 @@ SourceFiles
namespace Foam
{
// Class forward declarations
class topoMapper;
class mapPolyMesh;
/*---------------------------------------------------------------------------*\
Class topoSurfaceMapper Declaration
\*---------------------------------------------------------------------------*/
@ -76,6 +73,9 @@ class topoSurfaceMapper
//- Is the mapping direct
bool direct_;
//- Size before mapping
mutable label sizeBeforeMapping_;
// Demand-driven private data
//- Direct addressing
@ -165,8 +165,10 @@ public:
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "topoSurfaceMapper.C"
# include "topoSurfaceMapperTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

89
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/fieldMapping/topoSurfaceMapperTemplates.C Normal file
View file

@ -0,0 +1,89 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Map the internal field
template <class Type>
void topoSurfaceMapper::mapInternalField
(
const word& fieldName,
Field<Type>& iF
) const
{
if (iF.size() != sizeBeforeMapping())
{
FatalErrorIn
(
"\n\n"
"void topoSurfaceMapper::mapInternalField<Type>\n"
"(\n"
" Field<Type>& iF\n"
") const\n"
) << "Incompatible size before mapping." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << iF.size() << nl
<< " map size: " << sizeBeforeMapping() << nl
<< abort(FatalError);
}
// Map the internal field
iF.autoMap(*this);
// Flip the flux
const labelList flipFaces = flipFaceFlux().toc();
forAll (flipFaces, i)
{
if (flipFaces[i] < iF.size())
{
iF[flipFaces[i]] *= -1.0;
}
else
{
FatalErrorIn
(
"\n\n"
"void topoSurfaceMapper::mapInternalField<Type>\n"
"(\n"
" Field<Type>& iF\n"
") const\n"
) << "Cannot flip boundary face fluxes." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << iF.size() << nl
<< " Face flip index: " << flipFaces[i] << nl
<< abort(FatalError);
}
}
}
} // End namespace Foam
// ************************************************************************* //

773
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOps.C
View file

@ -33,7 +33,7 @@ Author
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
\*---------------------------------------------------------------------------*/
#include "Time.H"
#include "meshOps.H"
@ -41,8 +41,8 @@ Author
#include "Pstream.H"
#include "triFace.H"
#include "IOmanip.H"
#include "HashSet.H"
#include "polyMesh.H"
#include "OFstream.H"
#include "triPointRef.H"
#include "tetPointRef.H"
@ -56,7 +56,7 @@ namespace meshOps
// Utility method to build a hull of cells
// connected to the edge (for 2D simplical meshes)
void constructPrismHull
inline void constructPrismHull
(
const label eIndex,
const UList<face>& faces,
@ -126,7 +126,7 @@ void constructPrismHull
// Utility method to build a hull of cells (and faces)
// around an edge (for 3D simplical meshes)
void constructHull
inline void constructHull
(
const label eIndex,
const UList<face>& faces,
@ -136,7 +136,7 @@ void constructHull
const UList<label>& neighbour,
const UList<labelList>& faceEdges,
const UList<labelList>& edgeFaces,
const UList<labelList>& edgePoints,
const labelList& vertexHull,
labelList& hullEdges,
labelList& hullFaces,
labelList& hullCells,
@ -145,9 +145,9 @@ void constructHull
{
// [1] hullEdges is an ordered list of edge-labels around eIndex,
// but not connected to it.
// - Ordering is in the same manner as edgePoints.
// - Ordering is in the same manner as vertexHull.
// [2] hullFaces is an ordered list of face-labels connected to eIndex.
// - Ordering is in the same manner as edgePoints.
// - Ordering is in the same manner as vertexHull.
// [3] hullCells is an ordered list of cell-labels connected to eIndex.
// - For boundary hulls, the last cell label is -1
// [4] ringEntities are edges and faces connected to eIndex[0] and eIndex[1]
@ -162,7 +162,6 @@ void constructHull
// Obtain a reference to this edge, and its edgeFaces
const edge& edgeToCheck = edges[eIndex];
const labelList& eFaces = edgeFaces[eIndex];
const labelList& hullVertices = edgePoints[eIndex];
// Loop through all faces of this edge and add them to hullFaces
forAll(eFaces, faceI)
@ -170,7 +169,7 @@ void constructHull
const face& faceToCheck = faces[eFaces[faceI]];
// Find the isolated point on this face,
// and compare it with hullVertices
// and compare it with vertexHull
meshOps::findIsolatedPoint
(
faceToCheck,
@ -181,9 +180,9 @@ void constructHull
found = false;
forAll(hullVertices, indexI)
forAll(vertexHull, indexI)
{
if (hullVertices[indexI] == otherPoint)
if (vertexHull[indexI] == otherPoint)
{
// Fill in the position of this face on the hull
hullFaces[indexI] = eFaces[faceI];
@ -235,8 +234,8 @@ void constructHull
if (hullCells[indexI] != -1)
{
label nextI = hullVertices.fcIndex(indexI);
label nextHullPoint = hullVertices[nextI];
label nextI = vertexHull.fcIndex(indexI);
label nextHullPoint = vertexHull[nextI];
const cell& currCell = cells[hullCells[indexI]];
// Look for the ring-faces
@ -316,7 +315,7 @@ void constructHull
<< " edgeFaces connectivity is inconsistent. " << nl
<< " Edge: " << eIndex << ":: " << edgeToCheck << nl
<< " edgeFaces: " << eFaces << nl
<< " edgePoints: " << hullVertices
<< " vertexHull: " << vertexHull
<< abort(FatalError);
}
}
@ -334,7 +333,7 @@ void constructHull
// Renaud Waldura
// Dijkstra's Shortest Path Algorithm in Java
// http://renaud.waldura.com/
bool Dijkstra
inline bool Dijkstra
(
const Map<point>& points,
const Map<edge>& edges,
@ -464,13 +463,47 @@ bool Dijkstra
}
}
/*
// Write out the path
if (debug > 3)
{
if (foundEndPoint)
{
DynamicList<label> pathNodes(50);
label currentPoint = endPoint;
while (currentPoint != startPoint)
{
pathNodes.append(currentPoint);
currentPoint = pi[currentPoint];
}
pathNodes.append(startPoint);
pathNodes.shrink();
writeVTK
(
"DijkstraPath_"
+ Foam::name(startPoint)
+ '_'
+ Foam::name(endPoint),
pathNodes,
0
);
}
}
*/
return foundEndPoint;
}
// Select a list of elements from connectivity,
// and output to a VTK format
void writeVTK
inline void writeVTK
(
const polyMesh& mesh,
const word& name,
@ -480,7 +513,10 @@ void writeVTK
const UList<edge>& edges,
const UList<face>& faces,
const UList<cell>& cells,
const UList<label>& owner
const UList<label>& owner,
const UList<scalar>& scalField,
const UList<label>& lablField,
const UList<vector>& vectField
)
{
label nTotalCells = 0;
@ -502,236 +538,320 @@ void writeVTK
continue;
}
// Are we looking at points?
if (primitiveType == 0)
bool isPolyhedron = false;
switch (primitiveType)
{
// Size the list
cpList[nCells].setSize(1);
cpList[nCells] = cList[cellI];
nTotalCells++;
}
// Are we looking at edges?
if (primitiveType == 1)
{
// Size the list
cpList[nCells].setSize(2);
const edge& tEdge = edges[cList[cellI]];
cpList[nCells][0] = tEdge[0];
cpList[nCells][1] = tEdge[1];
nTotalCells += 2;
}
// Are we looking at faces?
if (primitiveType == 2)
{
const face& tFace = faces[cList[cellI]];
if (tFace.size() == 3)
// Are we looking at points?
case 0:
{
// Size the list
cpList[nCells].setSize(3);
cpList[nCells].setSize(1);
// Write out in order
cpList[nCells][0] = tFace[0];
cpList[nCells][1] = tFace[1];
cpList[nCells][2] = tFace[2];
cpList[nCells] = cList[cellI];
nTotalCells += 3;
nTotalCells++;
break;
}
else
if (tFace.size() == 4)
// Are we looking at edges?
case 1:
{
// Size the list
cpList[nCells].setSize(4);
cpList[nCells].setSize(2);
// Write out in order
cpList[nCells][0] = tFace[0];
cpList[nCells][1] = tFace[1];
cpList[nCells][2] = tFace[2];
cpList[nCells][3] = tFace[3];
const edge& tEdge = edges[cList[cellI]];
nTotalCells += 4;
cpList[nCells][0] = tEdge[0];
cpList[nCells][1] = tEdge[1];
nTotalCells += 2;
break;
}
}
// Are we looking at cells?
if (primitiveType == 3)
{
const cell& tCell = cells[cList[cellI]];
if (tCell.size() == 4)
// Are we looking at faces?
case 2:
{
// Point-ordering for tetrahedra
const face& baseFace = faces[tCell[0]];
const face& checkFace = faces[tCell[1]];
const face& tFace = faces[cList[cellI]];
// Size the list
cpList[nCells].setSize(4);
// Size the list and assign
cpList[nCells] = tFace;
// Get the fourth point
label apexPoint =
(
meshOps::findIsolatedPoint(baseFace, checkFace)
);
nTotalCells += tFace.size();
// Something's wrong with connectivity.
if (apexPoint == -1)
{
FatalErrorIn
(
"void writeVTK\n"
"(\n"
" const polyMesh& mesh,\n"
" const word& name,\n"
" const labelList& cList,\n"
" const label primitiveType,\n"
" const UList<point>& points,\n"
" const UList<edge>& edges,\n"
" const UList<face>& faces,\n"
" const UList<cell>& cells,\n"
" const UList<label>& owner\n"
") const\n"
)
<< "Cell: " << cList[cellI]
<< ":: " << tCell
<< " has inconsistent connectivity."
<< abort(FatalError);
}
// Write-out in order
label ownCell = owner[tCell[0]];
if (ownCell == cList[cellI])
{
cpList[nCells][0] = baseFace[2];
cpList[nCells][1] = baseFace[1];
cpList[nCells][2] = baseFace[0];
cpList[nCells][3] = apexPoint;
}
else
{
cpList[nCells][0] = baseFace[0];
cpList[nCells][1] = baseFace[1];
cpList[nCells][2] = baseFace[2];
cpList[nCells][3] = apexPoint;
}
nTotalCells += 4;
break;
}
else
if (tCell.size() == 5)
// Are we looking at cells?
case 3:
{
// Point-ordering for wedge cells
label firstTriFace = -1;
const cell& tCell = cells[cList[cellI]];
// Size the list
cpList[nCells].setSize(6);
// Figure out triangle faces
forAll(tCell, faceI)
// Is face information provided?
if (faces.size())
{
const face& cFace = faces[tCell[faceI]];
if (cFace.size() == 3)
if (tCell.size() == 4)
{
if (firstTriFace == -1)
{
firstTriFace = tCell[faceI];
// Point-ordering for tetrahedra
const face& baseFace = faces[tCell[0]];
const face& checkFace = faces[tCell[1]];
// Right-handedness is assumed here.
// Tri-faces are always on the boundary.
cpList[nCells][0] = cFace[0];
cpList[nCells][1] = cFace[1];
cpList[nCells][2] = cFace[2];
// Size the list
cpList[nCells].setSize(4);
// Get the fourth point
label apexPoint =
(
meshOps::findIsolatedPoint(baseFace, checkFace)
);
// Something's wrong with connectivity.
if (apexPoint == -1)
{
FatalErrorIn
(
"void meshOps::writeVTK\n"
"(\n"
" const polyMesh& mesh,\n"
" const word& name,\n"
" const labelList& cList,\n"
" const label primitiveType,\n"
" const UList<point>& meshPoints,\n"
" const UList<edge>& edges,\n"
" const UList<face>& faces,\n"
" const UList<cell>& cells,\n"
" const UList<label>& owner\n"
" const UList<scalar>& scalField,\n"
" const UList<label>& lablField,\n"
" const UList<vector>& vectField\n"
")\n"
)
<< "Cell: " << cList[cellI]
<< ":: " << tCell
<< " has inconsistent connectivity."
<< abort(FatalError);
}
// Write-out in order
label ownCell = owner[tCell[0]];
if (ownCell == cList[cellI])
{
cpList[nCells][0] = baseFace[2];
cpList[nCells][1] = baseFace[1];
cpList[nCells][2] = baseFace[0];
cpList[nCells][3] = apexPoint;
}
else
{
// Detect the three other points.
forAll(tCell, faceJ)
cpList[nCells][0] = baseFace[0];
cpList[nCells][1] = baseFace[1];
cpList[nCells][2] = baseFace[2];
cpList[nCells][3] = apexPoint;
}
nTotalCells += 4;
}
else
if (tCell.size() > 4)
{
// Point ordering for polyhedra
isPolyhedron = true;
// First obtain the face count
label npF = 0;
forAll(tCell, faceI)
{
npF += faces[tCell[faceI]].size();
}
// Size the list
cpList[nCells].setSize(tCell.size() + npF + 1);
// Fill in facePoints
label nP = 0;
// Fill in the number of faces
cpList[nCells][nP++] = tCell.size();
forAll(tCell, faceI)
{
const face& checkFace = faces[tCell[faceI]];
if (checkFace.empty())
{
const face& nFace = faces[tCell[faceJ]];
if (nFace.size() == 4)
{
// Search for vertices on cFace
// in this face.
forAll(cFace, I)
{
label i = nFace.which(cFace[I]);
if (i != -1)
{
label p = nFace.prevLabel(i);
label n = nFace.nextLabel(i);
if (p == cpList[nCells][0])
{
cpList[nCells][3] = cFace[I];
}
if (p == cpList[nCells][1])
{
cpList[nCells][4] = cFace[I];
}
if (p == cpList[nCells][2])
{
cpList[nCells][5] = cFace[I];
}
if (n == cpList[nCells][0])
{
cpList[nCells][3] = cFace[I];
}
if (n == cpList[nCells][1])
{
cpList[nCells][4] = cFace[I];
}
if (n == cpList[nCells][2])
{
cpList[nCells][5] = cFace[I];
}
}
}
}
FatalErrorIn("void meshOps::writeVTK()")
<< " Empty face: " << tCell[faceI]
<< abort(FatalError);
}
break;
// First fill in face size
cpList[nCells][nP++] = checkFace.size();
// Next fill in points in order
if (owner[tCell[faceI]] == cList[cellI])
{
forAll(checkFace, pI)
{
cpList[nCells][nP++] = checkFace[pI];
}
}
else
{
forAllReverse(checkFace, pI)
{
cpList[nCells][nP++] = checkFace[pI];
}
}
}
nTotalCells += (tCell.size() + npF + 1);
}
else
{
FatalErrorIn("void meshOps::writeVTK()")
<< " Wrong cell size: " << tCell << nl
<< " Index: " << cList[cellI] << nl
<< abort(FatalError);
}
}
else
{
// No face information.
// Assume cell-to-node information.
if (tCell.size() == 4)
{
// Build a face out of first three points.
triPointRef tpr
(
meshPoints[tCell[0]],
meshPoints[tCell[1]],
meshPoints[tCell[2]]
);
// Fetch fourth point
const point& d = meshPoints[tCell[3]];
scalar dDotN = ((d - tpr.a()) & tpr.normal());
// Size the list
cpList[nCells].setSize(4);
if (dDotN > 0.0)
{
// Correct orientation
cpList[nCells][0] = tCell[0];
cpList[nCells][1] = tCell[1];
cpList[nCells][2] = tCell[2];
cpList[nCells][3] = tCell[3];
}
else
{
// Flip triangle
cpList[nCells][0] = tCell[2];
cpList[nCells][1] = tCell[1];
cpList[nCells][2] = tCell[0];
cpList[nCells][3] = tCell[3];
}
nTotalCells += 4;
}
else
{
FatalErrorIn("void meshOps::writeVTK()")
<< " Wrong cell size: " << tCell << nl
<< " Index: " << cList[cellI] << nl
<< abort(FatalError);
}
}
nTotalCells += 6;
break;
}
default:
{
FatalErrorIn("void meshOps::writeVTK() const")
<< " Incorrect primitiveType: "
<< primitiveType
<< abort(FatalError);
}
}
// Renumber to local ordering
forAll(cpList[nCells], pointI)
if (isPolyhedron)
{
// Check if this point was added to the map
if (!pointMap.found(cpList[nCells][pointI]))
// Polyhedra also contain face sizes, so use a different scheme
label pCtr = 0;
// Fetch number of faces
label npF = cpList[nCells][pCtr++];
for (label i = 0; i < npF; i++)
{
// Point was not found, so add it
points[nPoints] = meshPoints[cpList[nCells][pointI]];
// Fetch number of points
label npP = cpList[nCells][pCtr++];
// Update the map
pointMap.insert(cpList[nCells][pointI], nPoints);
reversePointMap.insert(nPoints, cpList[nCells][pointI]);
// Read points in order
for (label j = 0; j < npP; j++)
{
// Check if this point was added to the map
if (!pointMap.found(cpList[nCells][pCtr]))
{
// Resize pointField if necessary
if (nPoints >= points.size())
{
points.setSize(2 * nPoints);
}
// Increment the number of points
nPoints++;
// Point was not found, so add it
points[nPoints] = meshPoints[cpList[nCells][pCtr]];
// Update the map
pointMap.insert(cpList[nCells][pCtr], nPoints);
reversePointMap.insert(nPoints, cpList[nCells][pCtr]);
// Increment the number of points
nPoints++;
}
// Renumber it.
cpList[nCells][pCtr] = pointMap[cpList[nCells][pCtr]];
// Increment point counter
pCtr++;
}
}
}
else
{
forAll(cpList[nCells], pointI)
{
// Check if this point was added to the map
if (!pointMap.found(cpList[nCells][pointI]))
{
// Resize pointField if necessary
if (nPoints >= points.size())
{
points.setSize(2 * nPoints);
}
// Renumber it.
cpList[nCells][pointI] = pointMap[cpList[nCells][pointI]];
// Point was not found, so add it
points[nPoints] = meshPoints[cpList[nCells][pointI]];
// Update the map
pointMap.insert(cpList[nCells][pointI], nPoints);
reversePointMap.insert(nPoints, cpList[nCells][pointI]);
// Increment the number of points
nPoints++;
}
// Renumber it.
cpList[nCells][pointI] = pointMap[cpList[nCells][pointI]];
}
}
// Update the cell map.
@ -752,13 +872,16 @@ void writeVTK
cpList,
primitiveType,
reversePointMap,
reverseCellMap
reverseCellMap,
scalField,
lablField,
vectField
);
}
// Actual routine to write out the VTK file
void writeVTK
inline void writeVTK
(
const polyMesh& mesh,
const word& name,
@ -769,7 +892,10 @@ void writeVTK
const labelListList& cpList,
const label primitiveType,
const Map<label>& reversePointMap,
const Map<label>& reverseCellMap
const Map<label>& reverseCellMap,
const UList<scalar>& scalField,
const UList<label>& lablField,
const UList<vector>& vectField
)
{
// Make the directory
@ -830,48 +956,54 @@ void writeVTK
{
forAll(cpList, i)
{
if (cpList[i].size() == 1)
switch (primitiveType)
{
// Vertex
file << "1" << nl;
}
case 0:
{
file << "1" << nl;
break;
}
if (cpList[i].size() == 2)
{
// Edge
file << "3" << nl;
}
case 1:
{
file << "3" << nl;
break;
}
if (cpList[i].size() == 3)
{
// Triangle face
file << "5" << nl;
}
// General Polygonal Face
case 2:
{
file << "7" << nl;
break;
}
if
(
(cpList[i].size() == 4) &&
(primitiveType == 2)
)
{
// Quad face
file << "9" << nl;
}
// Cell
case 3:
{
if (cpList[i].size() == 4)
{
// Tetrahedron
file << "10" << nl;
}
else
if (cpList[i].size() > 4)
{
// Polyhedron
file << "42" << nl;
}
if
(
(cpList[i].size() == 4) &&
(primitiveType == 3)
)
{
// Tetrahedron
file << "10" << nl;
}
break;
}
if (cpList[i].size() == 6)
{
// Wedge
file << "13" << nl;
default:
{
FatalErrorIn("void meshOps::writeVTK() const")
<< " Incorrect primitiveType: "
<< primitiveType
<< abort(FatalError);
}
}
}
}
@ -885,44 +1017,169 @@ void writeVTK
}
}
// Write out indices for visualization.
label nCellFields = 0, nPointFields = 0;
if (reverseCellMap.size())
{
file << "CELL_DATA " << nCells << endl;
nCellFields++;
}
file << "FIELD CellFields 1" << endl;
if (scalField.size() == nCells && scalField.size())
{
nCellFields++;
}
file << "CellIds 1 " << nCells << " int" << endl;
if (lablField.size() == nCells && lablField.size())
{
nCellFields++;
}
for (label i = 0; i < nCells; i++)
{
file << reverseCellMap[i] << ' ';
}
file << endl;
if (vectField.size() == nCells && vectField.size())
{
nCellFields++;
}
// Write out indices for visualization.
if (nCellFields)
{
file << "CELL_DATA " << nCells << endl;
file << "FIELD CellFields " << nCellFields << endl;
if (reverseCellMap.size())
{
file << "CellIds 1 " << nCells << " int" << endl;
for (label i = 0; i < nCells; i++)
{
file << reverseCellMap[i] << ' ';
}
file << endl;
}
if (scalField.size() == nCells)
{
file << "CellScalars 1 " << nCells << " double" << endl;
forAll(scalField, i)
{
file << scalField[i] << ' ';
}
file << endl;
}
if (lablField.size() == nCells)
{
file << "CellLabels 1 " << nCells << " int" << endl;
forAll(lablField, i)
{
file << lablField[i] << ' ';
}
file << endl;
}
if (vectField.size() == nCells)
{
file << "CellVectors 3 " << nCells << " double" << endl;
forAll(vectField, i)
{
file << vectField[i].x() << ' '
<< vectField[i].y() << ' '
<< vectField[i].z() << ' ';
}
file << endl;
}
}
if (reversePointMap.size())
{
nPointFields++;
}
if (scalField.size() == nPoints && scalField.size())
{
nPointFields++;
}
if (lablField.size() == nPoints && lablField.size())
{
nPointFields++;
}
if (vectField.size() == nPoints && vectField.size())
{
nPointFields++;
}
// Write out indices for visualization.
if (nPointFields)
{
file << "POINT_DATA " << nPoints << endl;
file << "FIELD PointFields 1" << endl;
file << "FIELD PointFields " << nPointFields << endl;
file << "PointIds 1 " << nPoints << " int" << endl;
for (label i = 0; i < nPoints; i++)
if (reversePointMap.size())
{
file << reversePointMap[i] << ' ';
file << "PointIds 1 " << nPoints << " int" << endl;
for (label i = 0; i < nPoints; i++)
{
file << reversePointMap[i] << ' ';
}
file << endl;
}
file << endl;
if (scalField.size() == nPoints)
{
file << "PointScalars 1 " << nPoints << " double" << endl;
forAll(scalField, i)
{
file << scalField[i] << ' ';
}
file << endl;
}
if (lablField.size() == nPoints)
{
file << "PointLabels 1 " << nPoints << " int" << endl;
forAll(lablField, i)
{
file << lablField[i] << ' ';
}
file << endl;
}
if (vectField.size() == nPoints)
{
file << "PointVectors 3 " << nPoints << " double" << endl;
forAll(vectField, i)
{
file << vectField[i].x() << ' '
<< vectField[i].y() << ' '
<< vectField[i].z() << ' ';
}
file << endl;
}
}
}
} // End namespace meshOps
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

120
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOps.H
View file

@ -43,7 +43,6 @@ SourceFiles
#define meshOps_H
#include "Map.H"
#include "line.H"
#include "point.H"
#include "label.H"
#include "scalar.H"
@ -51,8 +50,10 @@ SourceFiles
#include "cellList.H"
#include "edgeList.H"
#include "faceList.H"
#include "triangle.H"
#include "triPointRef.H"
#include "tetPointRef.H"
#include "vectorField.H"
#include "linePointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -93,7 +94,7 @@ namespace meshOps
// Utility method to build a hull of cells
// connected to the edge (for 2D simplical meshes)
void constructPrismHull
inline void constructPrismHull
(
const label eIndex,
const UList<face>& faces,
@ -107,7 +108,7 @@ namespace meshOps
// Utility method to build a hull of cells (and faces)
// around an edge (for 3D simplical meshes)
void constructHull
inline void constructHull
(
const label eIndex,
const UList<face>& faces,
@ -117,7 +118,7 @@ namespace meshOps
const UList<label>& neighbour,
const UList<labelList>& faceEdges,
const UList<labelList>& edgeFaces,
const UList<labelList>& edgePoints,
const labelList& hullVertices,
labelList& hullEdges,
labelList& hullFaces,
labelList& hullCells,
@ -147,12 +148,12 @@ namespace meshOps
(
const label cIndex,
const label fIndex,
const UList<face> faces,
const UList<cell> cells
const UList<face>& faces,
const UList<cell>& cells
);
// Given a cell index, find the centroid / volume of a cell.
inline bool cellCentreAndVolume
inline void cellCentreAndVolume
(
const label cIndex,
const vectorField& points,
@ -160,33 +161,28 @@ namespace meshOps
const UList<cell>& cells,
const UList<label>& owner,
vector& centre,
scalar& volume,
bool checkClosed = false
scalar& volume
);
// Determine whether a segment intersects a triangular face
template <class T>
inline bool segmentTriFaceIntersection
(
const triangle<Vector<T>, const Vector<T>&>& faceToCheck,
const line<Vector<T>, const Vector<T>&>& edgeToCheck,
const T& matchTol,
Vector<T>& intPoint
const triPointRef& faceToCheck,
const linePointRef& edgeToCheck,
vector& intPoint
);
// Determine whether the particular point lies
// inside the given triangular face
template <class T>
inline bool pointInTriFace
(
const triangle<Vector<T>, const Vector<T>&>& faceToCheck,
const Vector<T>& checkPoint,
const T& matchTol,
const triPointRef& faceToCheck,
const vector& checkPoint,
bool testCoplanarity
);
// Dijkstra's algorithm for the shortest path problem
bool Dijkstra
inline bool Dijkstra
(
const Map<point>& points,
const Map<edge>& edges,
@ -195,17 +191,6 @@ namespace meshOps
Map<label>& pi
);
// Method to insert labels in a face, so that
// right-handedness is preserved.
template <class T>
inline void insertPointLabels
(
const Vector<T>& refNorm,
const Field<Vector<T> >& points,
const labelHashSet& pLabels,
face& modFace
);
// Method to insert a label between two labels in a list
inline void insertLabel
(
@ -247,23 +232,75 @@ namespace meshOps
List<Type>& list
);
// Parallel send
inline void pWrite
(
const label toID,
const label& data
);
// Parallel send (for fixed lists)
template <class Type, unsigned Size>
inline void pWrite
(
const label toID,
const FixedList<Type,Size>& data
);
// Parallel send (for lists)
template <class Type>
inline void pWrite
(
const label toID,
const UList<Type>& data
);
// Parallel receive
inline void pRead
(
const label fromID,
label& data
);
// Parallel receive (for fixed lists)
template <class Type, unsigned Size>
inline void pRead
(
const label fromID,
FixedList<Type,Size>& data
);
// Parallel receive (for lists)
template <class Type>
inline void pRead
(
const label fromID,
UList<Type>& data
);
// Wait for buffer transfer completion.
inline void waitForBuffers();
// Select a list of elements from connectivity,
// and output to a VTK format
void writeVTK
inline void writeVTK
(
const polyMesh& mesh,
const word& name,
const labelList& cList,
const label primitiveType,
const UList<point>& meshPoints,
const UList<edge>& edges,
const UList<face>& faces,
const UList<cell>& cells,
const UList<label>& owner
const UList<edge>& edges = List<edge>(0),
const UList<face>& faces = List<face>(0),
const UList<cell>& cells = List<cell>(0),
const UList<label>& owner = List<label>(0),
const UList<scalar>& scalField = UList<scalar>(),
const UList<label>& lablField = UList<label>(),
const UList<vector>& vectField = UList<vector>()
);
// Actual routine to write out the VTK file
void writeVTK
inline void writeVTK
(
const polyMesh& mesh,
const word& name,
@ -274,7 +311,10 @@ namespace meshOps
const labelListList& cpList = labelListList(0),
const label primitiveType = 3,
const Map<label>& reversePointMap = Map<label>(),
const Map<label>& reverseCellMap = Map<label>()
const Map<label>& reverseCellMap = Map<label>(),
const UList<scalar>& scalField = UList<scalar>(),
const UList<label>& lablField = UList<label>(),
const UList<vector>& vectField = UList<vector>()
);
} // End namespace meshOps
@ -287,6 +327,10 @@ namespace meshOps
#include "meshOpsI.H"
#ifdef NoRepository
# include "meshOps.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

328
src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOpsI.H
View file

@ -216,8 +216,8 @@ inline label tetApexPoint
(
const label cIndex,
const label fIndex,
const UList<face> faces,
const UList<cell> cells
const UList<face>& faces,
const UList<cell>& cells
)
{
label apexPoint = -1;
@ -241,13 +241,13 @@ inline label tetApexPoint
if (!foundApex)
{
Info << "fIndex: " << fIndex << ":: " << faces[fIndex] << endl;
Info << "cIndex: " << cIndex << ":: " << cellToCheck << endl;
Pout<< "fIndex: " << fIndex << ":: " << faces[fIndex] << endl;
Pout<< "cIndex: " << cIndex << ":: " << cellToCheck << endl;
forAll(cellToCheck, faceI)
{
Info << '\t' << cellToCheck[faceI] << ":: "
<< faces[cellToCheck[faceI]] << endl;
Pout<< '\t' << cellToCheck[faceI] << ":: "
<< faces[cellToCheck[faceI]] << endl;
}
FatalErrorIn
@ -257,8 +257,8 @@ inline label tetApexPoint
"(\n"
" const label cIndex,\n"
" const label fIndex,\n"
" const UList<face> faces,\n"
" const UList<cell> cells\n"
" const UList<face>& faces,\n"
" const UList<cell>& cells\n"
")\n"
)
<< "Could not find an apex point in cell: " << cIndex
@ -272,7 +272,7 @@ inline label tetApexPoint
// Given a cell index, find the centroid / volume of a cell.
// - If checking is enabled, return whether the cell is closed
inline bool cellCentreAndVolume
inline void cellCentreAndVolume
(
const label cIndex,
const vectorField& points,
@ -280,8 +280,7 @@ inline bool cellCentreAndVolume
const UList<cell>& cells,
const UList<label>& owner,
vector& centre,
scalar& volume,
bool checkClosed
scalar& volume
)
{
// Reset inputs
@ -291,18 +290,12 @@ inline bool cellCentreAndVolume
const cell& cellToCheck = cells[cIndex];
// Average face-centres to get an estimate centroid
vector cEst(vector::zero), fCentre(vector::zero), fArea(vector::zero);
vector sumClosed(vector::zero), sumMagClosed(vector::zero);
vector cEst(vector::zero), fCentre(vector::zero);
forAll(cellToCheck, faceI)
{
const face& checkFace = faces[cellToCheck[faceI]];
if (checkFace.empty())
{
continue;
}
cEst += checkFace.centre(points);
}
@ -312,93 +305,82 @@ inline bool cellCentreAndVolume
{
const face& checkFace = faces[cellToCheck[faceI]];
if (checkFace.empty())
if (checkFace.size() == 3)
{
continue;
}
fArea = checkFace.normal(points);
fCentre = checkFace.centre(points);
// Flip if necessary
if (owner[cellToCheck[faceI]] != cIndex)
{
fArea *= -1.0;
}
// Calculate 3*face-pyramid volume
scalar pyr3Vol = fArea & (fCentre - cEst);
// Calculate face-pyramid centre
vector pc = ((3.0 / 4.0) * fCentre) + ((1.0 / 4.0) * cEst);
// Accumulate volume-weighted face-pyramid centre
centre += pyr3Vol*pc;
// Accumulate face-pyramid volume
volume += pyr3Vol;
// Accumulate areas, if necessary
if (checkClosed)
{
sumClosed += fArea;
sumMagClosed += cmptMag(fArea);
}
}
centre /= volume + VSMALL;
volume *= (1.0 / 3.0);
if (checkClosed)
{
// Check the sum across components
scalar closed = 0.0;
for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
{
closed =
tetPointRef tpr
(
mag(sumClosed[cmpt]) / (sumMagClosed[cmpt] + VSMALL)
points[checkFace[2]],
points[checkFace[1]],
points[checkFace[0]],
cEst
);
if (closed > 1e-06)
scalar tetVol = tpr.mag();
// Flip if necessary
if (owner[cellToCheck[faceI]] != cIndex)
{
// Volume is invalid
return false;
tetVol *= -1.0;
}
// Accumulate volume-weighted tet centre
centre += tetVol*tpr.centre();
// Accumulate tet volume
volume += tetVol;
}
else
{
forAll(checkFace, pI)
{
tetPointRef tpr
(
points[checkFace[pI]],
points[checkFace.prevLabel(pI)],
checkFace.centre(points),
cEst
);
scalar tetVol = tpr.mag();
// Flip if necessary
if (owner[cellToCheck[faceI]] != cIndex)
{
tetVol *= -1.0;
}
// Accumulate volume-weighted tet centre
centre += tetVol*tpr.centre();
// Accumulate tet volume
volume += tetVol;
}
}
}
// Volume is valid
return true;
centre /= volume + VSMALL;
}
// Determine whether a segment intersects a triangular face
template <class T>
inline bool segmentTriFaceIntersection
(
const triangle<Vector<T>, const Vector<T>&>& faceToCheck,
const line<Vector<T>, const Vector<T>&>& edgeToCheck,
const T& matchTol,
Vector<T>& intPoint
const triPointRef& faceToCheck,
const linePointRef& edgeToCheck,
vector& intPoint
)
{
// Fetch references
const Vector<T>& p1 = edgeToCheck.start();
const Vector<T>& p2 = edgeToCheck.end();
const vector& p1 = edgeToCheck.start();
const vector& p2 = edgeToCheck.end();
const Vector<T>& a = faceToCheck.a();
const Vector<T>& b = faceToCheck.b();
const Vector<T>& c = faceToCheck.c();
// Define custom face-normal
Vector<T> n = (T(0.5) * ((b - a)^(c - a)));
// Define face-normal
vector n = faceToCheck.normal();
n /= mag(n) + VSMALL;
// Compute uValue
T numerator = n & (a - p1);
T denominator = n & (p2 - p1);
scalar numerator = n & (faceToCheck.a() - p1);
scalar denominator = n & (p2 - p1);
// Check if the edge is parallel to the face
if (mag(denominator) < VSMALL)
@ -406,17 +388,16 @@ inline bool segmentTriFaceIntersection
return false;
}
T u = (numerator / denominator);
T tolerance = (matchTol * mag(p2 - p1));
scalar u = (numerator / denominator);
// Check for intersection along line.
if ((u > tolerance) && (u < (pTraits<T>::one - tolerance)))
if ((u >= 0.0) && (u <= 1.0))
{
// Compute point of intersection
intPoint = p1 + u*(p2 - p1);
// Also make sure that intPoint lies within the face
if (pointInTriFace(faceToCheck, intPoint, matchTol, false))
if (pointInTriFace(faceToCheck, intPoint, false))
{
return true;
}
@ -429,49 +410,45 @@ inline bool segmentTriFaceIntersection
// Determine whether the particular point lies
// inside the given triangular face
template <class T>
inline bool pointInTriFace
(
const triangle<Vector<T>, const Vector<T>&>& faceToCheck,
const Vector<T>& cP,
const T& matchTol,
const triPointRef& faceToCheck,
const vector& cP,
bool testCoplanarity
)
{
// Copy inputs
const Vector<T>& a = faceToCheck.a();
const Vector<T>& b = faceToCheck.b();
const Vector<T>& c = faceToCheck.c();
const T pointFive(0.5);
const vector& a = faceToCheck.a();
const vector& b = faceToCheck.b();
const vector& c = faceToCheck.c();
// Compute the normal
Vector<T> nf = (pointFive * ((b - a)^(c - a)));
vector nf = faceToCheck.normal();
if ( ((pointFive * ((b - a)^(cP - a))) & nf) < pTraits<T>::zero)
if ( ((0.5 * ((b - a)^(cP - a))) & nf) < 0.0)
{
return false;
}
if ( ((pointFive * ((c - b)^(cP - b))) & nf) < pTraits<T>::zero)
if ( ((0.5 * ((c - b)^(cP - b))) & nf) < 0.0)
{
return false;
}
if ( ((pointFive * ((a - c)^(cP - c))) & nf) < pTraits<T>::zero)
if ( ((0.5 * ((a - c)^(cP - c))) & nf) < 0.0)
{
return false;
}
if (testCoplanarity)
{
Vector<T> ftp(a - cP);
vector ftp(a - cP);
// Normalize vectors
nf /= mag(nf) + VSMALL;
ftp /= mag(ftp) + VSMALL;
if (mag(ftp & nf) > matchTol)
if (mag(ftp & nf) > VSMALL)
{
return false;
}
@ -482,53 +459,120 @@ inline bool pointInTriFace
}
// Method to insert labels in a face, so that
// right-handedness is preserved.
template <class T>
inline void insertPointLabels
// Parallel blocking send
inline void pWrite
(
const Vector<T>& refNorm,
const Field<Vector<T> >& points,
const labelHashSet& pLabels,
face& modFace
const label toID,
const label& data
)
{
// Need to configure a new face.
face newFace(modFace);
OPstream::write
(
Pstream::blocking,
toID,
reinterpret_cast<const char*>(&data),
sizeof(label)
);
}
const T pointFive(0.5);
forAllConstIter(labelHashSet, pLabels, pIter)
// Parallel blocking receive
inline void pRead
(
const label fromID,
label& data
)
{
IPstream::read
(
Pstream::blocking,
fromID,
reinterpret_cast<char*>(&data),
sizeof(label)
);
}
// Parallel non-blocking send for fixed lists
template <class Type, unsigned Size>
inline void pWrite
(
const label toID,
const FixedList<Type, Size>& data
)
{
OPstream::write
(
Pstream::blocking,
toID,
reinterpret_cast<const char*>(&data[0]),
Size*sizeof(Type)
);
}
// Parallel non-blocking receive for fixed lists
template <class Type, unsigned Size>
inline void pRead
(
const label fromID,
FixedList<Type, Size>& data
)
{
IPstream::read
(
Pstream::blocking,
fromID,
reinterpret_cast<char*>(data.begin()),
Size*sizeof(Type)
);
}
// Parallel non-blocking send for lists
template <class Type>
inline void pWrite
(
const label toID,
const UList<Type>& data
)
{
OPstream::write
(
Pstream::nonBlocking,
toID,
reinterpret_cast<const char*>(&data[0]),
data.size()*sizeof(Type)
);
}
// Parallel non-blocking receive for lists
template <class Type>
inline void pRead
(
const label fromID,
UList<Type>& data
)
{
IPstream::read
(
Pstream::nonBlocking,
fromID,
reinterpret_cast<char*>(&data[0]),
data.size()*sizeof(Type)
);
}
// Wait for buffer transfer completion.
inline void waitForBuffers()
{
if (Pstream::parRun())
{
forAll(newFace, pI)
{
label nI = newFace.fcIndex(pI);
const Vector<T>& a = points[newFace[pI]];
const Vector<T>& b = points[pIter.key()];
const Vector<T>& c = points[newFace[nI]];
// Compute the normal.
Vector<T> newNorm = (pointFive * ((b - a)^(c - a)));
if ((refNorm & newNorm) > pTraits<T>::zero)
{
// Insert the point.
meshOps::insertLabel
(
pIter.key(),
newFace[pI],
newFace[nI],
newFace
);
break;
}
}
OPstream::waitRequests();
IPstream::waitRequests();
}
// Take over storage
modFace.transfer(newFace);
}