Merge remote-tracking branch 'origin/feature/parallelDynamicTopoFvMesh' into nextRelease

Conflicts:
	src/dynamicMesh/dynamicFvMesh/Make/files
	src/dynamicMesh/dynamicFvMesh/Make/options
	src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C
	src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/eMesh/eMeshDemandDrivenData.C
This commit is contained in:
Hrvoje Jasak 2013-10-25 14:30:27 +01:00
commit e6b91c08cc
51 changed files with 29586 additions and 4498 deletions

8
.gitignore vendored
View file

@ -106,4 +106,12 @@ src/lduSolvers/amg/amgPolicy/samgPolicy.H
src/lduSolvers/amg/amgPolicy/aamgPolicy.C src/lduSolvers/amg/amgPolicy/aamgPolicy.C
src/lduSolvers/amg/amgPolicy/aamgPolicy.H src/lduSolvers/amg/amgPolicy/aamgPolicy.H
# The following files are blacklisted because of a DMCA complaint by ANSYS.
src/lduSolvers/tools/PriorityArray.C
src/lduSolvers/tools/PriorityArray.H
src/lduSolvers/amg/amgPolicy/samgPolicy.C
src/lduSolvers/amg/amgPolicy/samgPolicy.H
src/lduSolvers/amg/amgPolicy/aamgPolicy.C
src/lduSolvers/amg/amgPolicy/aamgPolicy.H
# end-of-file # end-of-file

View file

@ -40,4 +40,18 @@ dynamicTopoFvMesh/edgeSwap.C
dynamicTopoFvMesh/edgeBisect.C dynamicTopoFvMesh/edgeBisect.C
dynamicTopoFvMesh/edgeCollapse.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
LIB = $(FOAM_LIBBIN)/libdynamicFvMesh LIB = $(FOAM_LIBBIN)/libdynamicFvMesh

View file

@ -41,109 +41,140 @@ SourceFiles
#ifndef changeMap_H #ifndef changeMap_H
#define changeMap_H #define changeMap_H
#include "objectMap.H"
#include "labelList.H" #include "labelList.H"
#include "FixedList.H" #include "dictionary.H"
#include "DynamicList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam namespace Foam
{ {
// Class forward declarations
class changeMap;
// * * * * * * * * Forward declaration of friend fuctions * * * * * * * * * //
Ostream& operator<<(Ostream&, const changeMap&);
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class changeMap Declaration Class changeMap Declaration
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
class changeMap class changeMap
:
public dictionary
{ {
// Sliver type index. // Entity index
// Type 1: Sliver label index_;
// Type 2: Cap
// Type 3: Spade // Coupled patch index
// Type 4: Wedge label pIndex_;
// Type
label 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. // Entities that were added during the operation.
// - Add a master mapping for slaves. DynamicList<objectMap> addedPoints_;
// - Masters get a map of -1. DynamicList<objectMap> addedEdges_;
List<FixedList<label,2> > addedPoints_; DynamicList<objectMap> addedFaces_;
List<FixedList<label,2> > addedEdges_; DynamicList<objectMap> addedCells_;
List<FixedList<label,2> > addedFaces_;
List<FixedList<label,2> > addedCells_; // Entities that were removed during the operation
DynamicList<label> removedPoints_;
DynamicList<label> removedEdges_;
DynamicList<label> removedFaces_;
DynamicList<label> removedCells_;
public: public:
// Constructor // Constructor
changeMap() changeMap()
: :
index_(-1),
pIndex_(-1),
type_(-1), type_(-1),
firstEdge_(-1), addedPoints_(5),
secondEdge_(-1), addedEdges_(5),
apexPoint_(-1), addedFaces_(5),
opposingFace_(-1), addedCells_(5),
addedPoints_(0), removedPoints_(5),
addedEdges_(0), removedEdges_(5),
addedFaces_(0), removedFaces_(5),
addedCells_(0) removedCells_(5)
{} {}
//- Access //- Access
// Entity index
inline label& index();
inline label index() const;
// Coupled patch index
inline label& patchIndex();
inline label patchIndex() const;
// Type // Type
inline label& type(); inline label& type();
inline label type() const; 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 // Added entities
inline void addPoint inline void addPoint
( (
const label pIndex, const label pIndex,
const label master = -1 const labelList& master = labelList()
); );
inline void addEdge inline void addEdge
( (
const label eIndex, const label eIndex,
const label master = -1 const labelList& master = labelList()
); );
inline void addFace inline void addFace
( (
const label fIndex, const label fIndex,
const label master = -1 const labelList& master = labelList()
); );
inline void addCell inline void addCell
( (
const label cIndex, const label cIndex,
const label master = -1 const labelList& master = labelList()
); );
// Return the list of added entities // Return the list of added entities
inline const List<FixedList<label,2> >& addedPointList() const; inline const List<objectMap>& addedPointList() const;
inline const List<FixedList<label,2> >& addedEdgeList() const; inline const List<objectMap>& addedEdgeList() const;
inline const List<FixedList<label,2> >& addedFaceList() const; inline const List<objectMap>& addedFaceList() const;
inline const List<FixedList<label,2> >& addedCellList() 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 //- Operators
inline void operator=(const changeMap& rhs); inline void operator=(const changeMap& rhs);
//- IOstream Operators
inline friend Ostream& operator<<(Ostream&, const changeMap&);
}; };

View file

@ -22,6 +22,17 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation, along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 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 namespace Foam
@ -29,6 +40,32 @@ namespace Foam
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * 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 // Type
inline label& changeMap::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 inline void changeMap::addPoint
( (
const label pIndex, const label pIndex,
const label master const labelList& master
) )
{ {
label curSize = addedPoints_.size(); addedPoints_.append(objectMap(pIndex, master));
addedPoints_.setSize(curSize + 1);
addedPoints_[curSize][0] = pIndex;
addedPoints_[curSize][1] = master;
} }
inline void changeMap::addEdge inline void changeMap::addEdge
( (
const label eIndex, const label eIndex,
const label master const labelList& master
) )
{ {
label curSize = addedEdges_.size(); addedEdges_.append(objectMap(eIndex, master));
addedEdges_.setSize(curSize + 1);
addedEdges_[curSize][0] = eIndex;
addedEdges_[curSize][1] = master;
} }
inline void changeMap::addFace inline void changeMap::addFace
( (
const label fIndex, const label fIndex,
const label master const labelList& master
) )
{ {
label curSize = addedFaces_.size(); addedFaces_.append(objectMap(fIndex, master));
addedFaces_.setSize(curSize + 1);
addedFaces_[curSize][0] = fIndex;
addedFaces_[curSize][1] = master;
} }
inline void changeMap::addCell inline void changeMap::addCell
( (
const label cIndex, const label cIndex,
const label master const labelList& master
) )
{ {
label curSize = addedCells_.size(); addedCells_.append(objectMap(cIndex, master));
addedCells_.setSize(curSize + 1);
addedCells_[curSize][0] = cIndex;
addedCells_[curSize][1] = master;
} }
// Return an added point inline const List<objectMap>&
inline const List<FixedList<label,2> >&
changeMap::addedPointList() const changeMap::addedPointList() const
{ {
return addedPoints_; return addedPoints_;
} }
// Return the list of added entities inline const List<objectMap>&
inline const List<FixedList<label,2> >&
changeMap::addedEdgeList() const changeMap::addedEdgeList() const
{ {
return addedEdges_; return addedEdges_;
} }
inline const List<FixedList<label,2> >& inline const List<objectMap>&
changeMap::addedFaceList() const changeMap::addedFaceList() const
{ {
return addedFaces_; return addedFaces_;
} }
inline const List<FixedList<label,2> >& inline const List<objectMap>&
changeMap::addedCellList() const changeMap::addedCellList() const
{ {
return addedCells_; 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) inline void changeMap::operator=(const changeMap& rhs)
{ {
// Copy base dictionary
dictionary::operator=(rhs);
index_ = rhs.index_;
pIndex_ = rhs.pIndex_;
type_ = rhs.type_; type_ = rhs.type_;
firstEdge_ = rhs.firstEdge_; // Copy maps
secondEdge_ = rhs.secondEdge_; addedPoints_.setSize(rhs.addedPoints_.size());
apexPoint_ = rhs.apexPoint_;
opposingFace_ = rhs.opposingFace_;
addedPoints_ = rhs.addedPoints_; forAll(addedPoints_, indexI)
addedEdges_ = rhs.addedEdges_; {
addedFaces_ = rhs.addedFaces_; addedPoints_[indexI].index() = rhs.addedPoints_[indexI].index();
addedCells_ = rhs.addedCells_;
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 } // End namespace Foam
// ************************************************************************* // // ************************************************************************* //

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
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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
// ************************************************************************* //

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
// ************************************************************************* //

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
// ************************************************************************* //

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
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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
// ************************************************************************* //

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
// ************************************************************************* //

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;
}
}
// ************************************************************************* //

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
// ************************************************************************* //

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;
}
}
// ************************************************************************* //

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
// ************************************************************************* //

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
// ************************************************************************* //

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
// ************************************************************************* //

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
// ************************************************************************* //

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
// ************************************************************************* //

View file

@ -37,6 +37,7 @@ SourceFiles
dynamicTopoFvMesh.C dynamicTopoFvMesh.C
dynamicTopoFvMeshI.H dynamicTopoFvMeshI.H
dynamicTopoFvMeshCheck.C dynamicTopoFvMeshCheck.C
dynamicTopoFvMeshCoupled.C
dynamicTopoFvMeshReOrder.C dynamicTopoFvMeshReOrder.C
dynamicTopoFvMeshMapping.C dynamicTopoFvMeshMapping.C
edgeBisect.C edgeBisect.C
@ -50,6 +51,7 @@ SourceFiles
#include "Switch.H" #include "Switch.H"
#include "tetMetric.H" #include "tetMetric.H"
#include "topoMapper.H"
#include "DynamicField.H" #include "DynamicField.H"
#include "threadHandler.H" #include "threadHandler.H"
#include "dynamicFvMesh.H" #include "dynamicFvMesh.H"
@ -64,8 +66,9 @@ class eMesh;
class Stack; class Stack;
class changeMap; class changeMap;
class objectMap; class objectMap;
class topoMapper; class coupledInfo;
class motionSolver; class motionSolver;
class convexSetAlgorithm;
class lengthScaleEstimator; class lengthScaleEstimator;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
@ -78,9 +81,18 @@ class dynamicTopoFvMesh
{ {
// Private data // Private data
//- Reference to base mesh
const dynamicTopoFvMesh& baseMesh_;
//- Number of boundary patches
label nPatches_;
//- Topology change flag //- Topology change flag
bool topoChangeFlag_; bool topoChangeFlag_;
//- Is this a subMesh?
bool isSubMesh_;
//- Dynamic mesh dictionary //- Dynamic mesh dictionary
IOdictionary dict_; IOdictionary dict_;
@ -94,9 +106,15 @@ class dynamicTopoFvMesh
//- Edge refinement switch //- Edge refinement switch
Switch edgeRefinement_; Switch edgeRefinement_;
//- Load motion solver
Switch loadMotionSolver_;
//- Switch for cell-bandwidth reduction //- Switch for cell-bandwidth reduction
Switch bandWidthReduction_; Switch bandWidthReduction_;
//- Coupled modification switch
mutable Switch coupledModification_;
//- Specify the re-meshing interval //- Specify the re-meshing interval
label interval_; label interval_;
@ -119,27 +137,21 @@ class dynamicTopoFvMesh
// support template typedefs, // support template typedefs,
// so this is a work-around // so this is a work-around
template <class T> template <class T>
class resizableList class resizable
{ {
public: 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> resizable<point>::FieldType oldPoints_, points_;
class resizableField resizable<face>::ListType faces_;
{ resizable<label>::ListType owner_, neighbour_;
public: resizable<cell>::ListType cells_;
typedef DynamicField<T> Type; resizable<edge>::ListType edges_;
}; resizable<labelList>::ListType pointEdges_;
resizable<labelList>::ListType edgeFaces_, faceEdges_;
resizableField<point>::Type oldPoints_, points_; resizable<scalar>::ListType lengthScale_;
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_;
//- Size information //- Size information
labelList oldPatchSizes_, patchSizes_; labelList oldPatchSizes_, patchSizes_;
@ -181,6 +193,7 @@ class dynamicTopoFvMesh
Map<labelList> faceParents_; Map<labelList> faceParents_;
List<scalarField> faceWeights_; List<scalarField> faceWeights_;
List<vectorField> faceCentres_; List<vectorField> faceCentres_;
Map<labelList> cellParents_; Map<labelList> cellParents_;
List<scalarField> cellWeights_; List<scalarField> cellWeights_;
List<vectorField> cellCentres_; List<vectorField> cellCentres_;
@ -201,9 +214,8 @@ class dynamicTopoFvMesh
labelHashSet deletedFaces_; labelHashSet deletedFaces_;
labelHashSet deletedCells_; labelHashSet deletedCells_;
//- List of flipped faces / modified old-points //- List of flipped faces
labelHashSet flipFaces_; labelHashSet flipFaces_;
Map<labelList> modPoints_;
//- Run-time statistics //- Run-time statistics
label maxModifications_; label maxModifications_;
@ -237,11 +249,42 @@ class dynamicTopoFvMesh
// in multi-threaded reOrdering // in multi-threaded reOrdering
FixedList<Mutex, 4> entityMutex_; 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 // Private Member Functions
//- Disallow default bitwise copy construct //- Disallow default bitwise copy construct
dynamicTopoFvMesh(const dynamicTopoFvMesh&); 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 //- Disallow default bitwise assignment
void operator=(const dynamicTopoFvMesh&); void operator=(const dynamicTopoFvMesh&);
@ -256,6 +299,7 @@ class dynamicTopoFvMesh
( (
const scalar matchTol, const scalar matchTol,
const bool skipMapping, const bool skipMapping,
const bool mappingOutput,
const label faceStart, const label faceStart,
const label faceSize, const label faceSize,
const label cellStart, const label cellStart,
@ -265,19 +309,13 @@ class dynamicTopoFvMesh
// Static equivalent for multiThreading // Static equivalent for multiThreading
static void computeMappingThread(void *argument); 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 // 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 // Initialize mesh edges and related connectivity lists
void initEdges(); void initEdges();
@ -316,11 +354,21 @@ class dynamicTopoFvMesh
// contains a triangular face. // contains a triangular face.
inline label getTriBoundaryEdge(const label fIndex) const; 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] // 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] // 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 // Check for bisection
inline bool checkBisection(const label index) const; inline bool checkBisection(const label index) const;
@ -349,6 +397,13 @@ class dynamicTopoFvMesh
// Initialize stacks // Initialize stacks
inline void initStacks(const labelHashSet& entities); 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 // Method to determine the boundary patch index for a given face
inline label whichPatch(const label index) const; inline label whichPatch(const label index) const;
@ -385,13 +440,6 @@ class dynamicTopoFvMesh
bool forceOp = false 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 // Method for the bisection of an edge in 3D
const changeMap const changeMap
bisectEdge bisectEdge
@ -420,6 +468,9 @@ class dynamicTopoFvMesh
bool forceOp = false bool forceOp = false
); );
// Slice the mesh at a particular location
void sliceMesh(const labelPair& pointPair);
// Utility method to check for invalid face-bisection. // Utility method to check for invalid face-bisection.
bool checkBisection bool checkBisection
( (
@ -429,8 +480,9 @@ class dynamicTopoFvMesh
) const; ) const;
// Utility method to check for invalid face-collapse. // Utility method to check for invalid face-collapse.
bool checkCollapse static bool checkCollapse
( (
const dynamicTopoFvMesh& mesh,
const labelList& triFaces, const labelList& triFaces,
const FixedList<label,2>& c0BdyIndex, const FixedList<label,2>& c0BdyIndex,
const FixedList<label,2>& c1BdyIndex, const FixedList<label,2>& c1BdyIndex,
@ -440,7 +492,7 @@ class dynamicTopoFvMesh
scalar& collapseQuality, scalar& collapseQuality,
const bool checkNeighbour, const bool checkNeighbour,
bool forceOp = false bool forceOp = false
) const; );
// Utility method to check for invalid edge-collapse. // Utility method to check for invalid edge-collapse.
bool checkCollapse bool checkCollapse
@ -449,7 +501,7 @@ class dynamicTopoFvMesh
const point& oldPoint, const point& oldPoint,
const label pointIndex, const label pointIndex,
const label cellIndex, const label cellIndex,
labelHashSet& cellsChecked, DynamicList<label>& cellsChecked,
scalar& collapseQuality, scalar& collapseQuality,
bool forceOp = false bool forceOp = false
) const; ) const;
@ -491,7 +543,9 @@ class dynamicTopoFvMesh
removeCells removeCells
( (
const labelList& cList, const labelList& cList,
const label patch const label patch,
const word& rcName,
bool checkOnly = false
); );
// Remove the specified cell from the mesh // Remove the specified cell from the mesh
@ -525,13 +579,19 @@ class dynamicTopoFvMesh
const Map<bool>& cellColors 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 // Insert the specified edge to the mesh
label insertEdge label insertEdge
( (
const label patch, const label patch,
const edge& newEdge, const edge& newEdge,
const labelList& edgeFaces, const labelList& edgeFaces
const labelList& edgePoints = labelList::null()
); );
// Remove the specified edge from the mesh // Remove the specified edge from the mesh
@ -555,12 +615,13 @@ class dynamicTopoFvMesh
const label pIndex const label pIndex
); );
// Utility method to build edgePoints for an edge [3D] // Utility method to build vertexHull for an edge [3D]
void buildEdgePoints void buildVertexHull
( (
const label eIndex, const label eIndex,
labelList& vertexHull,
const label checkIndex = 0 const label checkIndex = 0
); ) const;
// Utility to check whether points of an edge lie on a boundary // Utility to check whether points of an edge lie on a boundary
const FixedList<bool,2> checkEdgeBoundary(const label eIndex) const; const FixedList<bool,2> checkEdgeBoundary(const label eIndex) const;
@ -569,7 +630,20 @@ class dynamicTopoFvMesh
inline void setFlip(const label fIndex); inline void setFlip(const label fIndex);
// Utility method to compute the minimum quality of a vertex hull // 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 // Utility method to compute the quality of a vertex hull
// around an edge after bisection. // around an edge after bisection.
@ -580,7 +654,12 @@ class dynamicTopoFvMesh
scalar computeTrisectionQuality(const label fIndex) const; scalar computeTrisectionQuality(const label fIndex) const;
// Check whether the given edge is on a bounding curve // 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 // Allocate dynamic programming tables
void initTables void initTables
@ -606,14 +685,28 @@ class dynamicTopoFvMesh
bool fillTables bool fillTables
( (
const label eIndex, const label eIndex,
const scalar minQuality, scalar& minQuality,
labelList& m, labelList& m,
labelList& hullVertices,
PtrList<scalarListList>& Q, PtrList<scalarListList>& Q,
PtrList<labelListList>& K, PtrList<labelListList>& K,
PtrList<labelListList>& triangulations, PtrList<labelListList>& triangulations,
const label checkIndex = 0 const label checkIndex = 0
) const; ) 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 // Print out tables for debugging
void printTables void printTables
( (
@ -629,7 +722,9 @@ class dynamicTopoFvMesh
( (
const label eIndex, const label eIndex,
const scalar minQuality, const scalar minQuality,
const PtrList<labelListList>& K, const labelList& hullVertices,
PtrList<scalarListList>& Q,
PtrList<labelListList>& K,
PtrList<labelListList>& triangulations, PtrList<labelListList>& triangulations,
const label checkIndex = 0 const label checkIndex = 0
); );
@ -649,7 +744,8 @@ class dynamicTopoFvMesh
( (
const label eIndex, const label eIndex,
const labelList& hullVertices, const labelList& hullVertices,
const labelListList& triangulations const labelListList& triangulations,
bool output = false
) const; ) const;
// Check old-volumes for the configuration. // Check old-volumes for the configuration.
@ -781,6 +877,9 @@ class dynamicTopoFvMesh
// Reset the mesh and generate mapping information // Reset the mesh and generate mapping information
bool resetMesh(); bool resetMesh();
// Write out connectivity for an edge
void writeEdgeConnectivity(const label eIndex) const;
// Output an entity as a VTK file // Output an entity as a VTK file
void writeVTK void writeVTK
( (
@ -798,7 +897,9 @@ class dynamicTopoFvMesh
const labelList& cList, const labelList& cList,
const label primitiveType = 3, const label primitiveType = 3,
const bool useOldConnectivity = false, 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; ) const;
// Return the status report interval // Return the status report interval
@ -807,6 +908,18 @@ class dynamicTopoFvMesh
// Check the state of local connectivity lists // Check the state of local connectivity lists
void checkConnectivity(const label maxErrors = 0) const; 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. // Test an edge / face for proximity with other non-neighbouring faces.
// Return the scalar distance to the nearest face. // Return the scalar distance to the nearest face.
scalar testProximity scalar testProximity
@ -818,9 +931,165 @@ class dynamicTopoFvMesh
// Calculate the edge length-scale for the mesh // Calculate the edge length-scale for the mesh
void calculateLengthScale(bool dump = false); 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 // Read optional dictionary parameters
void readOptionalParameters(bool reRead = false); 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 // Dump cell-quality statistics
bool meshQuality(bool outputOption); bool meshQuality(bool outputOption);
@ -840,8 +1109,11 @@ public:
// Member Functions // Member Functions
//- Update mesh corresponding to the given map
virtual void updateMesh(const mapPolyMesh& mpm);
// Map all fields in time using a customized mapper // 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 // Update the mesh for motion / topology changes
// Return true if topology changes have occurred // Return true if topology changes have occurred

File diff suppressed because it is too large Load diff

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 // Check for edge bisection
inline bool dynamicTopoFvMesh::checkBisection inline bool dynamicTopoFvMesh::checkBisection
( (
@ -475,6 +209,12 @@ inline bool dynamicTopoFvMesh::checkBisection
if (twoDMesh_) if (twoDMesh_)
{ {
// If this entity was deleted, skip it.
if (faces_[index].empty())
{
return false;
}
// Fetch the edge // Fetch the edge
const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)]; const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
@ -493,6 +233,12 @@ inline bool dynamicTopoFvMesh::checkBisection
} }
else else
{ {
// If this entity was deleted, skip it.
if (edgeFaces_[index].empty())
{
return false;
}
// Fetch the edge // Fetch the edge
const edge& edgeToCheck = edges_[index]; const edge& edgeToCheck = edges_[index];
@ -531,6 +277,12 @@ inline bool dynamicTopoFvMesh::checkCollapse
if (twoDMesh_) if (twoDMesh_)
{ {
// If this entity was deleted, skip it.
if (faces_[index].empty())
{
return false;
}
// Fetch the edge // Fetch the edge
const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)]; const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
@ -549,6 +301,12 @@ inline bool dynamicTopoFvMesh::checkCollapse
} }
else else
{ {
// If this entity was deleted, skip it.
if (edgeFaces_[index].empty())
{
return false;
}
// Fetch the edge // Fetch the edge
const edge& edgeToCheck = edges_[index]; const edge& edgeToCheck = edges_[index];
@ -637,6 +395,15 @@ inline void dynamicTopoFvMesh::initStacks
{ {
forAll(faces_, faceI) forAll(faces_, faceI)
{ {
// For coupled meshes, avoid certain faces.
if (patchCoupling_.size() || procIndices_.size())
{
if (entities.found(faceI))
{
continue;
}
}
if (faces_[faceI].size() == 4) if (faces_[faceI].size() == 4)
{ {
stack(tID[tIndex]).insert(faceI); stack(tID[tIndex]).insert(faceI);
@ -649,6 +416,15 @@ inline void dynamicTopoFvMesh::initStacks
{ {
forAll(edges_, edgeI) forAll(edges_, edgeI)
{ {
// For coupled meshes, avoid certain edges.
if (patchCoupling_.size() || procIndices_.size())
{
if (entities.found(edgeI))
{
continue;
}
}
if (edgeFaces_[edgeI].size()) if (edgeFaces_[edgeI].size())
{ {
stack(tID[tIndex]).insert(edgeI); 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 // 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 // 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 else
{ {
@ -737,9 +541,11 @@ inline label dynamicTopoFvMesh::whichEdgePatch
// If not in any of the above, it's possible that the edge was added // 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 // 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 else
{ {
@ -799,13 +605,15 @@ inline void dynamicTopoFvMesh::setFlip(const label fIndex)
{ {
if (fIndex < nOldFaces_) if (fIndex < nOldFaces_)
{ {
if (flipFaces_.found(fIndex)) labelHashSet::iterator it = flipFaces_.find(fIndex);
if (it == flipFaces_.end())
{ {
flipFaces_.erase(fIndex); flipFaces_.insert(fIndex);
} }
else else
{ {
flipFaces_.insert(fIndex); flipFaces_.erase(it);
} }
} }
} }

View file

@ -37,14 +37,21 @@ Author
#include "dynamicTopoFvMesh.H" #include "dynamicTopoFvMesh.H"
#include "Time.H"
#include "meshOps.H" #include "meshOps.H"
#include "IOmanip.H" #include "IOmanip.H"
#include "triFace.H" #include "triFace.H"
#include "objectMap.H" #include "objectMap.H"
#include "faceSetAlgorithm.H"
#include "cellSetAlgorithm.H"
namespace Foam namespace Foam
{ {
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTemplateTypeNameAndDebugWithName(IOList<objectMap>, "objectMapIOList", 0);
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Compute mapping weights for modified entities // Compute mapping weights for modified entities
@ -52,77 +59,397 @@ void dynamicTopoFvMesh::computeMapping
( (
const scalar matchTol, const scalar matchTol,
const bool skipMapping, const bool skipMapping,
const bool mappingOutput,
const label faceStart, const label faceStart,
const label faceSize, const label faceSize,
const label cellStart, const label cellStart,
const label cellSize 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 // Compute cell mapping
for (label cellI = cellStart; cellI < (cellStart + cellSize); cellI++) for (label cellI = cellStart; cellI < (cellStart + cellSize); cellI++)
{ {
label cIndex = cellsFromCells_[cellI].index(); label cIndex = cellsFromCells_[cellI].index();
labelList& masterObjects = cellsFromCells_[cellI].masterObjects();
if (skipMapping) if (skipMapping)
{ {
// Set empty mapping parameters // Dummy map from cell[0]
const labelList& mo = cellParents_[cIndex]; masterObjects = labelList(1, 0);
cellWeights_[cellI].setSize(1, 1.0);
cellsFromCells_[cellI].masterObjects() = mo; cellCentres_[cellI].setSize(1, vector::zero);
cellWeights_[cellI].setSize(mo.size(), (1.0/(mo.size() + VSMALL)));
cellCentres_[cellI].setSize(mo.size(), vector::zero);
} }
else else
{ {
// Compute master objects for inverse-distance weighting // Obtain weighting factors for this cell.
computeParents cellAlgorithm.computeWeights
( (
cIndex, cIndex,
0,
cellParents_[cIndex], cellParents_[cIndex],
polyMesh::cellCells(), polyMesh::cellCells(),
polyMesh::cellCentres(), masterObjects,
3, cellWeights_[cellI],
cellsFromCells_[cellI].masterObjects() 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 // Compute face mapping
for (label faceI = faceStart; faceI < (faceStart + faceSize); faceI++) for (label faceI = faceStart; faceI < (faceStart + faceSize); faceI++)
{ {
label fIndex = facesFromFaces_[faceI].index(); label fIndex = facesFromFaces_[faceI].index();
label patchIndex = whichPatch(fIndex); labelList& masterObjects = facesFromFaces_[faceI].masterObjects();
// Skip mapping for internal faces. label patchIndex = whichPatch(fIndex);
if (patchIndex == -1) label neiProc = getNeighbourProcessor(patchIndex);
// Skip mapping for internal / processor faces.
if (patchIndex == -1 || neiProc > -1)
{ {
// Set dummy masters, so that the conventional // Set dummy masters, so that the conventional
// faceMapper doesn't incur a seg-fault. // faceMapper doesn't crash-and-burn
facesFromFaces_[faceI].masterObjects() = labelList(1, 0); masterObjects = labelList(1, 0);
continue; continue;
} }
if (skipMapping) if (skipMapping)
{ {
// Set empty mapping parameters // Dummy map from patch[0]
const labelList& mo = faceParents_[fIndex]; masterObjects = labelList(1, 0);
faceWeights_[faceI].setSize(1, 1.0);
facesFromFaces_[faceI].masterObjects() = mo; faceCentres_[faceI].setSize(1, vector::zero);
faceWeights_[faceI].setSize(mo.size(), (1.0/(mo.size() + VSMALL)));
faceCentres_[faceI].setSize(mo.size(), vector::zero);
} }
else else
{ {
// Compute master objects for inverse-distance weighting // Obtain weighting factors for this face.
computeParents faceAlgorithm.computeWeights
( (
fIndex, fIndex,
boundaryMesh()[patchIndex].start(),
faceParents_[fIndex], faceParents_[fIndex],
boundaryMesh()[patchIndex].faceFaces(), boundaryMesh()[patchIndex].faceFaces(),
boundaryMesh()[patchIndex].faceCentres(), masterObjects,
2, faceWeights_[faceI],
facesFromFaces_[faceI].masterObjects() 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 // Recast the pointers for the argument
scalar& matchTol = *(static_cast<scalar*>(thread->operator()(0))); scalar& matchTol = *(static_cast<scalar*>(thread->operator()(0)));
bool& skipMapping = *(static_cast<bool*>(thread->operator()(1))); bool& skipMapping = *(static_cast<bool*>(thread->operator()(1)));
label& faceStart = *(static_cast<label*>(thread->operator()(2))); bool& mappingOutput = *(static_cast<bool*>(thread->operator()(2)));
label& faceSize = *(static_cast<label*>(thread->operator()(3))); label& faceStart = *(static_cast<label*>(thread->operator()(3)));
label& cellStart = *(static_cast<label*>(thread->operator()(4))); label& faceSize = *(static_cast<label*>(thread->operator()(4)));
label& cellSize = *(static_cast<label*>(thread->operator()(5))); label& cellStart = *(static_cast<label*>(thread->operator()(5)));
label& cellSize = *(static_cast<label*>(thread->operator()(6)));
// Now calculate addressing // Now calculate addressing
mesh.computeMapping mesh.computeMapping
( (
matchTol, matchTol,
skipMapping, skipMapping,
mappingOutput,
faceStart, faceSize, faceStart, faceSize,
cellStart, cellSize cellStart, cellSize
); );
@ -169,7 +498,8 @@ void dynamicTopoFvMesh::computeMappingThread(void *argument)
void dynamicTopoFvMesh::threadedMapping void dynamicTopoFvMesh::threadedMapping
( (
scalar matchTol, scalar matchTol,
bool skipMapping bool skipMapping,
bool mappingOutput
) )
{ {
label nThreads = threader_->getNumThreads(); label nThreads = threader_->getNumThreads();
@ -177,7 +507,7 @@ void dynamicTopoFvMesh::threadedMapping
// If mapping is being skipped, issue a warning. // If mapping is being skipped, issue a warning.
if (skipMapping) if (skipMapping)
{ {
Info << " *** Mapping is being skipped *** " << endl; Info<< " *** Mapping is being skipped *** " << endl;
} }
// Check if single-threaded // Check if single-threaded
@ -187,6 +517,7 @@ void dynamicTopoFvMesh::threadedMapping
( (
matchTol, matchTol,
skipMapping, skipMapping,
mappingOutput,
0, facesFromFaces_.size(), 0, facesFromFaces_.size(),
0, cellsFromCells_.size() 0, cellsFromCells_.size()
); );
@ -212,8 +543,8 @@ void dynamicTopoFvMesh::threadedMapping
if (debug > 2) if (debug > 2)
{ {
Info << " Mapping Faces: " << index[0] << endl; Pout<< " Mapping Faces: " << index[0] << nl
Info << " Mapping Cells: " << index[1] << endl; << " Mapping Cells: " << index[1] << endl;
} }
forAll(index, indexI) forAll(index, indexI)
@ -234,8 +565,8 @@ void dynamicTopoFvMesh::threadedMapping
if (debug > 2) if (debug > 2)
{ {
Info << " Load starts: " << tStarts[indexI] << endl; Pout<< " Load starts: " << tStarts[indexI] << nl
Info << " Load sizes: " << tSizes[indexI] << endl; << " Load sizes: " << tSizes[indexI] << endl;
} }
} }
@ -243,7 +574,7 @@ void dynamicTopoFvMesh::threadedMapping
forAll(hdl, i) forAll(hdl, i)
{ {
// Size up the argument list // Size up the argument list
hdl[i].setSize(6); hdl[i].setSize(7);
// Set match tolerance // Set match tolerance
hdl[i].set(0, &matchTol); hdl[i].set(0, &matchTol);
@ -251,25 +582,26 @@ void dynamicTopoFvMesh::threadedMapping
// Set the skipMapping flag // Set the skipMapping flag
hdl[i].set(1, &skipMapping); hdl[i].set(1, &skipMapping);
// Set the mappingOutput flag
hdl[i].set(2, &mappingOutput);
// Set the start/size indices // Set the start/size indices
hdl[i].set(2, &(tStarts[0][i])); hdl[i].set(3, &(tStarts[0][i]));
hdl[i].set(3, &(tSizes[0][i])); hdl[i].set(4, &(tSizes[0][i]));
hdl[i].set(4, &(tStarts[1][i])); hdl[i].set(5, &(tStarts[1][i]));
hdl[i].set(5, &(tSizes[1][i])); hdl[i].set(6, &(tSizes[1][i]));
} }
// Prior to multi-threaded operation, // Prior to multi-threaded operation,
// force calculation of demand-driven data. // force calculation of demand-driven data.
polyMesh::cells(); polyMesh::cells();
primitiveMesh::cellCells(); primitiveMesh::cellCells();
primitiveMesh::cellCentres();
const polyBoundaryMesh& boundary = boundaryMesh(); const polyBoundaryMesh& boundary = boundaryMesh();
forAll(boundary, patchI) forAll(boundary, patchI)
{ {
boundary[patchI].faceFaces(); boundary[patchI].faceFaces();
boundary[patchI].faceCentres();
} }
// Execute threads in linear sequence // 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 // Set fill-in mapping information for a particular cell
void dynamicTopoFvMesh::setCellMapping void dynamicTopoFvMesh::setCellMapping
( (
@ -413,9 +621,9 @@ void dynamicTopoFvMesh::setCellMapping
{ {
if (debug > 3) if (debug > 3)
{ {
Info << "Inserting mapping cell: " << cIndex << nl Pout<< "Inserting mapping cell: " << cIndex
<< " mapCells: " << mapCells << " mapCells: " << mapCells
<< endl; << endl;
} }
// Insert index into the list, and overwrite if necessary // Insert index into the list, and overwrite if necessary
@ -445,7 +653,7 @@ void dynamicTopoFvMesh::setCellMapping
} }
// Update cell-parents information // Update cell-parents information
labelHashSet masterCells; DynamicList<label> masterCells(5);
forAll(mapCells, cellI) forAll(mapCells, cellI)
{ {
@ -456,7 +664,10 @@ void dynamicTopoFvMesh::setCellMapping
if (mapCells[cellI] < nOldCells_) if (mapCells[cellI] < nOldCells_)
{ {
masterCells.insert(mapCells[cellI]); if (findIndex(masterCells, mapCells[cellI]) == -1)
{
masterCells.append(mapCells[cellI]);
}
} }
else else
if (cellParents_.found(mapCells[cellI])) if (cellParents_.found(mapCells[cellI]))
@ -465,12 +676,15 @@ void dynamicTopoFvMesh::setCellMapping
forAll(nParents, cI) 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 patch = whichPatch(fIndex);
label neiProc = getNeighbourProcessor(patch);
if (debug > 3) if (debug > 3)
{ {
Info << "Inserting mapping face: " << fIndex << nl const polyBoundaryMesh& boundary = boundaryMesh();
<< " patch: " << patch << nl
<< " mapFaces: " << mapFaces word pName;
<< endl;
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; bool foundError = false;
@ -500,15 +734,16 @@ void dynamicTopoFvMesh::setFaceMapping
foundError = true; foundError = true;
} }
// Check to ensure that boundary faces map // Check to ensure that mapFaces is non-empty for physical patches
// only from other faces on the same patch if (patch > -1 && mapFaces.empty() && neiProc == -1)
if (patch > -1 && mapFaces.empty())
{ {
foundError = true; foundError = true;
} }
if (foundError) if (foundError)
{ {
writeVTK("mapFace_" + Foam::name(fIndex), fIndex, 2);
FatalErrorIn FatalErrorIn
( (
"\n" "\n"
@ -524,7 +759,8 @@ void dynamicTopoFvMesh::setFaceMapping
<< " 2. Mapping specified for an internal face, " << nl << " 2. Mapping specified for an internal face, " << nl
<< " when none was expected." << nl << nl << " when none was expected." << nl << nl
<< " Face: " << fIndex << nl << " Face: " << fIndex << nl
<< " Patch: " << patch << nl << " Patch: "
<< (patch > -1 ? boundaryMesh()[patch].name() : "Internal") << nl
<< " Owner: " << owner_[fIndex] << nl << " Owner: " << owner_[fIndex] << nl
<< " Neighbour: " << neighbour_[fIndex] << nl << " Neighbour: " << neighbour_[fIndex] << nl
<< " mapFaces: " << mapFaces << nl << " mapFaces: " << mapFaces << nl
@ -556,14 +792,14 @@ void dynamicTopoFvMesh::setFaceMapping
facesFromFaces_[index].masterObjects() = labelList(0); facesFromFaces_[index].masterObjects() = labelList(0);
} }
// For internal faces, set dummy maps / weights, and bail out // For internal / processor faces, bail out
if (patch == -1) if (patch == -1 || neiProc > -1)
{ {
return; return;
} }
// Update face-parents information // Update face-parents information
labelHashSet masterFaces; DynamicList<label> masterFaces(5);
forAll(mapFaces, faceI) forAll(mapFaces, faceI)
{ {
@ -574,7 +810,10 @@ void dynamicTopoFvMesh::setFaceMapping
if (mapFaces[faceI] < nOldFaces_) if (mapFaces[faceI] < nOldFaces_)
{ {
masterFaces.insert(mapFaces[faceI]); if (findIndex(masterFaces, mapFaces[faceI]) == -1)
{
masterFaces.append(mapFaces[faceI]);
}
} }
else else
if (faceParents_.found(mapFaces[faceI])) if (faceParents_.found(mapFaces[faceI]))
@ -583,12 +822,15 @@ void dynamicTopoFvMesh::setFaceMapping
forAll(nParents, fI) 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);
} }

View file

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

View file

@ -42,7 +42,6 @@ defineTypeNameAndDebug(eMesh, 0);
word eMesh::meshSubDir = "eMesh"; word eMesh::meshSubDir = "eMesh";
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void eMesh::clearGeom() const void eMesh::clearGeom() const

View file

@ -24,6 +24,8 @@ License
Description Description
Demand-driven edge mesh data
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "eMesh.H" #include "eMesh.H"
@ -179,12 +181,8 @@ void eMesh::calcEdgePoints() const
// If a boundary point is present, that must begin the list. // If a boundary point is present, that must begin the list.
// NOTE: Will work only on tetrahedral meshes! // NOTE: Will work only on tetrahedral meshes!
# ifdef FULLDEBUG
bool found; bool found;
# endif label faceIndex = -1, cellIndex = -1;
label faceIndex = -1;
label cellIndex = -1;
const labelList& owner = mesh_.faceOwner(); const labelList& owner = mesh_.faceOwner();
const labelList& neighbour = mesh_.faceNeighbour(); const labelList& neighbour = mesh_.faceNeighbour();
const cellList& cells = mesh_.cells(); const cellList& cells = mesh_.cells();
@ -249,9 +247,7 @@ void eMesh::calcEdgePoints() const
const cell& cellToCheck = cells[cellIndex]; const cell& cellToCheck = cells[cellIndex];
# ifdef FULLDEBUG
found = false; found = false;
# endif
// Assuming tet-cells, // Assuming tet-cells,
// Loop through edgeFaces and get the next face // Loop through edgeFaces and get the next face
@ -264,12 +260,7 @@ void eMesh::calcEdgePoints() const
) )
{ {
faceIndex = cellToCheck[0]; faceIndex = cellToCheck[0];
found = true; break;
# ifdef FULLDEBUG
found = true;
# endif
break;
} }
if if
@ -279,12 +270,7 @@ void eMesh::calcEdgePoints() const
) )
{ {
faceIndex = cellToCheck[1]; faceIndex = cellToCheck[1];
found = true; break;
# ifdef FULLDEBUG
found = true;
# endif
break;
} }
if if
@ -294,12 +280,7 @@ void eMesh::calcEdgePoints() const
) )
{ {
faceIndex = cellToCheck[2]; faceIndex = cellToCheck[2];
found = true; break;
# ifdef FULLDEBUG
found = true;
# endif
break;
} }
if if
@ -309,12 +290,7 @@ void eMesh::calcEdgePoints() const
) )
{ {
faceIndex = cellToCheck[3]; faceIndex = cellToCheck[3];
found = true; break;
# ifdef FULLDEBUG
found = true;
# endif
break;
} }
} }

View file

@ -80,16 +80,6 @@ ePatch::ePatch
size_(readLabel(dict.lookup("size"))) 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) ePatch::ePatch(const ePatch& p, const eBoundaryMesh& bm)
: :
patchIdentifier(p, p.index()), patchIdentifier(p, p.index()),

View file

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

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

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
// ************************************************************************* //

View file

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

View file

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

View file

@ -33,9 +33,8 @@ Author
University of Massachusetts Amherst University of Massachusetts Amherst
All rights reserved All rights reserved
\*----------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "IOmanip.H"
#include "topoMapper.H" #include "topoMapper.H"
#include "mapPolyMesh.H" #include "mapPolyMesh.H"
#include "topoCellMapper.H" #include "topoCellMapper.H"
@ -306,6 +305,7 @@ const List<vectorField>& topoCellMapper::intersectionCentres() const
return *centresPtr_; return *centresPtr_;
} }
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components // Construct from components
@ -319,6 +319,7 @@ topoCellMapper::topoCellMapper
mpm_(mpm), mpm_(mpm),
tMapper_(mapper), tMapper_(mapper),
direct_(false), direct_(false),
sizeBeforeMapping_(mpm.nOldCells()),
directAddrPtr_(NULL), directAddrPtr_(NULL),
interpolationAddrPtr_(NULL), interpolationAddrPtr_(NULL),
weightsPtr_(NULL), weightsPtr_(NULL),
@ -326,6 +327,18 @@ topoCellMapper::topoCellMapper
volumesPtr_(NULL), volumesPtr_(NULL),
centresPtr_(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 // Check for possibility of direct mapping
if if
( (
@ -363,7 +376,7 @@ label topoCellMapper::size() const
//- Return size before mapping //- Return size before mapping
label topoCellMapper::sizeBeforeMapping() const 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 * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void topoCellMapper::operator=(const topoCellMapper& rhs) void topoCellMapper::operator=(const topoCellMapper& rhs)

View file

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

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
// ************************************************************************* //

View file

@ -35,120 +35,18 @@ Author
\*----------------------------------------------------------------------------*/ \*----------------------------------------------------------------------------*/
#include "fvc.H"
#include "topoMapper.H" #include "topoMapper.H"
#include "fluxCorrector.H" #include "fluxCorrector.H"
#include "topoCellMapper.H" #include "topoCellMapper.H"
#include "leastSquaresGrad.H"
#include "topoSurfaceMapper.H" #include "topoSurfaceMapper.H"
#include "topoBoundaryMeshMapper.H" #include "topoBoundaryMeshMapper.H"
#include "fixedValueFvPatchFields.H"
namespace Foam namespace Foam
{ {
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * 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 //- Store gradients prior to mesh reset
void topoMapper::storeGradients() const void topoMapper::storeGradients() const
{ {
@ -157,11 +55,8 @@ void topoMapper::storeGradients() const
if (fvMesh::debug) if (fvMesh::debug)
{ {
Info << "Registered volScalarFields: " << endl; Info<< "Registered volScalarFields: " << scalarGrads() << endl;
Info << sGrads_.toc() << endl; Info<< "Registered volVectorFields: " << vectorGrads() << endl;
Info << "Registered volVectorFields: " << endl;
Info << vGrads_.toc() << endl;
} }
} }
@ -169,52 +64,78 @@ void topoMapper::storeGradients() const
//- Store geometric information //- Store geometric information
void topoMapper::storeGeometry() const 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_); volCentrePatches.set
deleteDemandDrivenData(cellCentresPtr_); (
patchI,
new fixedValueFvPatchField<vector>
(
mesh_.boundary()[patchI],
DimensionedField<vector, volMesh>::null()
)
);
patchAreasPtr_.clear(); // Slice field to patch (forced assignment)
patchCentresPtr_.clear(); volCentrePatches[patchI] ==
(
mesh_.boundaryMesh()[patchI].patchSlice(Cf)
);
} }
// Set the cell-volumes pointer.
cellVolumesPtr_ = new scalarField(mesh_.cellVolumes());
// Set the cell-centres pointer. // Set the cell-centres pointer.
cellCentresPtr_ = new vectorField(mesh_.cellCentres()); cellCentresPtr_ =
(
// Set patch-areas new volVectorField
patchAreasPtr_.setSize(mesh_.boundaryMesh().size());
forAll(mesh_.boundaryMesh(), patchI)
{
patchAreasPtr_.set
( (
patchI, IOobject
new scalarField
( (
mag(mesh_.boundaryMesh()[patchI].faceAreas()) "cellCentres",
) mesh_.time().timeName(),
); mesh_,
} IOobject::NO_READ,
IOobject::NO_WRITE,
// Set patch-centres. false
patchCentresPtr_.setSize(mesh_.boundaryMesh().size()); ),
mesh_,
forAll(mesh_.boundaryMesh(), patchI) dimLength,
{ SubField<vector>(Cv, mesh_.nCells()),
patchCentresPtr_.set volCentrePatches
( )
patchI, );
new vectorField
(
mesh_.boundaryMesh()[patchI].faceCentres()
)
);
}
} }
// * * * * * * * * * * * * * * * * 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 * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * 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 //- Fetch face weights
const List<scalarField>& topoMapper::faceWeights() const 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 //- Store mesh information for the mapping stage
void topoMapper::storeMeshInformation() const void topoMapper::storeMeshInformation() const
{ {
@ -325,20 +308,19 @@ void topoMapper::storeMeshInformation() const
storeGeometry(); storeGeometry();
} }
//- Return non-const access to cell centres
//- Return stored cell volume information volVectorField& topoMapper::volCentres() const
const scalarField& topoMapper::cellVolumes() const
{ {
if (!cellVolumesPtr_) if (!cellCentresPtr_)
{ {
FatalErrorIn FatalErrorIn
( (
"const scalarField& topoMapper::cellVolumes()" "const vectorField& topoMapper::volCentres() const"
) << nl << " Pointer has not been set. " ) << nl << " Pointer has not been set. "
<< abort(FatalError); << abort(FatalError);
} }
return *cellVolumesPtr_; return *cellCentresPtr_;
} }
@ -349,7 +331,7 @@ const vectorField& topoMapper::internalCentres() const
{ {
FatalErrorIn FatalErrorIn
( (
"const vectorField& topoMapper::internalCentres()" "const vectorField& topoMapper::internalCentres() const"
) << nl << " Pointer has not been set. " ) << nl << " Pointer has not been set. "
<< abort(FatalError); << 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 //- Return stored patch centre information
const vectorField& topoMapper::patchCentres(const label i) const const vectorField& topoMapper::patchCentres(const label i) const
{ {
if (!patchCentresPtr_.set(i)) if (!cellCentresPtr_)
{ {
FatalErrorIn FatalErrorIn
( (
"const vectorField& topoMapper::patchCentres" "const vectorField& topoMapper::patchCentres"
"(const label i) const" "(const label i) const"
) << nl << " Pointer has not been set at index: " << i ) << nl << " Pointer has not been set. index: " << i
<< abort(FatalError); << abort(FatalError);
} }
return patchCentresPtr_[i]; return (*cellCentresPtr_).boundaryField()[i];
} }
// Conservatively map all volFields in the registry //- Return names of stored scalar gradients
template <class Type> const wordList topoMapper::scalarGrads() const
void topoMapper::conservativeMapVolFields() const
{ {
// Define a few typedefs for convenience return sGrads_.toc();
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();
}
} }
// Conservatively map all surfaceFields in the registry //- Return names of stored vector gradients
template <class Type> const wordList topoMapper::vectorGrads() const
void topoMapper::conservativeMapSurfaceFields() const
{ {
// Define a few typedefs for convenience return vGrads_.toc();
typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfType; }
HashTable<const surfType*> fields(mesh_.lookupClass<surfType>());
// Store old-times before mapping //- Fetch the gradient field (template specialisation)
forAllIter(typename HashTable<const surfType*>, fields, fIter) template <>
volVectorField& topoMapper::gradient(const word& name) const
{
if (!sGrads_.found(name))
{ {
surfType& field = const_cast<surfType&>(*fIter()); FatalErrorIn
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
( (
field.name(), "volVectorField& topoMapper::gradient(const word& name) const"
field.internalField() ) << nl << " Gradient for: " << name
); << " has not been stored."
<< abort(FatalError);
// 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_[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 FatalErrorIn
( (
"const topoCellMapper& topoMapper::volMap()" "const topoCellMapper& topoMapper::volMap() const"
) << nl << " Volume mapper has not been set. " ) << nl << " Volume mapper has not been set. "
<< abort(FatalError); << abort(FatalError);
} }
@ -549,7 +447,7 @@ const topoSurfaceMapper& topoMapper::surfaceMap() const
{ {
FatalErrorIn FatalErrorIn
( (
"const topoSurfaceMapper& topoMapper::surfaceMap()" "const topoSurfaceMapper& topoMapper::surfaceMap() const"
) << nl << " Surface mapper has not been set. " ) << nl << " Surface mapper has not been set. "
<< abort(FatalError); << abort(FatalError);
} }
@ -565,7 +463,7 @@ const topoBoundaryMeshMapper& topoMapper::boundaryMap() const
{ {
FatalErrorIn FatalErrorIn
( (
"const topoBoundaryMeshMapper& topoMapper::boundaryMap()" "const topoBoundaryMeshMapper& topoMapper::boundaryMap() const"
) << nl << " Boundary mapper has not been set. " ) << nl << " Boundary mapper has not been set. "
<< abort(FatalError); << abort(FatalError);
} }
@ -581,7 +479,7 @@ const fluxCorrector& topoMapper::surfaceFluxCorrector() const
{ {
FatalErrorIn FatalErrorIn
( (
"const fluxCorrector& topoMapper::surfaceFluxCorrector()" "const fluxCorrector& topoMapper::surfaceFluxCorrector() const"
) << nl << " fluxCorrector has not been set. " ) << nl << " fluxCorrector has not been set. "
<< abort(FatalError); << abort(FatalError);
} }
@ -591,7 +489,7 @@ const fluxCorrector& topoMapper::surfaceFluxCorrector() const
//- Clear out member data //- Clear out member data
void topoMapper::clear() void topoMapper::clear() const
{ {
// Clear out mappers // Clear out mappers
cellMap_.clear(); cellMap_.clear();
@ -603,10 +501,7 @@ void topoMapper::clear()
vGrads_.clear(); vGrads_.clear();
// Wipe out geomtry information // Wipe out geomtry information
deleteDemandDrivenData(cellVolumesPtr_);
deleteDemandDrivenData(cellCentresPtr_); deleteDemandDrivenData(cellCentresPtr_);
patchAreasPtr_.clear();
patchCentresPtr_.clear();
// Clear maps // Clear maps
faceWeights_.clear(); faceWeights_.clear();
@ -614,6 +509,16 @@ void topoMapper::clear()
faceCentres_.clear(); faceCentres_.clear();
cellCentres_.clear(); cellCentres_.clear();
// Clear sizes / offsets
cellSizes_.clear();
cellStarts_.clear();
faceSizes_.clear();
faceStarts_.clear();
patchSizes_.clear();
patchStarts_.clear();
} }

View file

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

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
// ************************************************************************* //

View file

@ -33,9 +33,10 @@ Author
University of Massachusetts Amherst University of Massachusetts Amherst
All rights reserved All rights reserved
\*----------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "IOmanip.H" #include "IOmanip.H"
#include "meshOps.H"
#include "topoMapper.H" #include "topoMapper.H"
#include "mapPolyMesh.H" #include "mapPolyMesh.H"
#include "topoPatchMapper.H" #include "topoPatchMapper.H"
@ -73,9 +74,17 @@ void topoPatchMapper::calcInsertedFaceAddressing() const
} }
// Information from the old patch // 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 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 // Allocate for inserted face labels and addressing
label nInsertedFaces = 0; label nInsertedFaces = 0;
@ -108,6 +117,19 @@ void topoPatchMapper::calcInsertedFaceAddressing() const
{ {
if (fffI.masterObjects().empty()) 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 FatalErrorIn
( (
"void topoPatchMapper::" "void topoPatchMapper::"
@ -130,23 +152,26 @@ void topoPatchMapper::calcInsertedFaceAddressing() const
// Make an entry for addressing // Make an entry for addressing
labelList& addr = insertedAddressing[nInsertedFaces]; labelList& addr = insertedAddressing[nInsertedFaces];
// Renumber addressing to patch. // Check for illegal addressing
// Also, check mapping for hits into
// other patches / internal faces.
addr = fffI.masterObjects(); addr = fffI.masterObjects();
forAll(addr, faceI) forAll(addr, faceI)
{ {
if if (addr[faceI] < 0 || addr[faceI] >= oldPatchSize)
(
addr[faceI] >= oldPatchStart
&& addr[faceI] < oldPatchEnd
)
{
addr[faceI] -= oldPatchStart;
}
else
{ {
// 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 FatalErrorIn
( (
"void topoPatchMapper::" "void topoPatchMapper::"
@ -154,10 +179,13 @@ void topoPatchMapper::calcInsertedFaceAddressing() const
) )
<< "Addressing into another patch is not allowed." << "Addressing into another patch is not allowed."
<< nl << " Patch face index: " << faceI << 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 << " addr[faceI]: " << addr[faceI]
<< nl << " oldPatchStart: " << oldPatchStart << nl << " oldPatchStart: " << oldPatchStart
<< nl << " oldPatchSize: " << oldPatchSize << nl << " oldPatchSize: " << oldPatchSize
<< nl << " oldPatchEnd: " << oldPatchEnd
<< abort(FatalError); << abort(FatalError);
} }
} }
@ -185,21 +213,31 @@ void topoPatchMapper::calcAddressing() const
} }
// Information from the old patch // 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 oldPatchStart = mpm_.oldPatchStarts()[patch_.index()];
const label oldPatchEnd = oldPatchStart + oldPatchSize; const label oldPatchEnd = oldPatchStart + oldPatchSize;
// Assemble the maps: slice to patch // Assemble the maps: slice to patch
if (direct()) if (direct())
{ {
if (oldPatchSize == 0)
{
// Nothing to map from. Return empty.
directAddrPtr_ = new labelList(0);
return;
}
// Direct mapping - slice to size // Direct mapping - slice to size
directAddrPtr_ = new labelList(patch_.patch().patchSlice(mpm_.faceMap())); directAddrPtr_ =
(
new labelList(patch_.patch().patchSlice(mpm_.faceMap()))
);
labelList& addr = *directAddrPtr_; labelList& addr = *directAddrPtr_;
// Shift to local patch indices. // Shift to local patch indices.
// Also, check mapping for hits into other patches / internal faces. // Also, check mapping for hits into other patches / internal faces.
forAll (addr, faceI) forAll(addr, faceI)
{ {
if if
( (
@ -211,22 +249,55 @@ void topoPatchMapper::calcAddressing() const
} }
else 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 FatalErrorIn
( (
"void topoPatchMapper::calcAddressing() const" "void topoPatchMapper::calcAddressing() const"
) )
<< "Addressing into another patch is not allowed." << "Addressing into another patch is not allowed."
<< nl << " Patch: " << patch_.name()
<< nl << " Patch index: " << patch_.index()
<< nl << " Patch face index: " << faceI << nl << " Patch face index: " << faceI
<< nl << " addr[faceI]: " << addr[faceI] << nl << " addr[faceI]: " << addr[faceI]
<< nl << " oldPatchStart: " << oldPatchStart << nl << " oldPatchStart: " << oldPatchStart
<< nl << " oldPatchSize: " << oldPatchSize << nl << " oldPatchSize: " << oldPatchSize
<< nl << " oldPatchEnd: " << oldPatchEnd << nl << " oldPatchEnd: " << oldPatchEnd
<< nl << " nInserted: " << insertedObjectLabels().size()
<< abort(FatalError); << abort(FatalError);
} }
} }
} }
else else
{ {
if (oldPatchSize == 0)
{
// Nothing to map from. Return empty.
interpolationAddrPtr_ = new labelListList(0);
return;
}
// Interpolative addressing // Interpolative addressing
interpolationAddrPtr_ = new labelListList(patchSize(), labelList(0)); interpolationAddrPtr_ = new labelListList(patchSize(), labelList(0));
labelListList& addr = *interpolationAddrPtr_; labelListList& addr = *interpolationAddrPtr_;
@ -267,6 +338,8 @@ void topoPatchMapper::calcAddressing() const
"void topoPatchMapper::calcAddressing() const" "void topoPatchMapper::calcAddressing() const"
) )
<< "Addressing into another patch is not allowed." << "Addressing into another patch is not allowed."
<< nl << " Patch: " << patch_.name()
<< nl << " Patch index: " << patch_.index()
<< nl << " Patch face index: " << faceI << nl << " Patch face index: " << faceI
<< nl << " faceMap[faceI]: " << oldFace << nl << " faceMap[faceI]: " << oldFace
<< nl << " oldPatchStart: " << oldPatchStart << nl << " oldPatchStart: " << oldPatchStart
@ -284,6 +357,17 @@ void topoPatchMapper::calcAddressing() const
{ {
if (addr[faceI].empty()) 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 FatalErrorIn
( (
"void topoPatchMapper::calcAddressing() const" "void topoPatchMapper::calcAddressing() const"
@ -316,6 +400,13 @@ void topoPatchMapper::calcInverseDistanceWeights() const
// Fetch interpolative addressing // Fetch interpolative addressing
const labelListList& addr = addressing(); const labelListList& addr = addressing();
if (addr.empty())
{
// Nothing to map from. Return empty.
weightsPtr_ = new scalarListList(0);
return;
}
// Allocate memory // Allocate memory
weightsPtr_ = new scalarListList(patchSize()); weightsPtr_ = new scalarListList(patchSize());
scalarListList& w = *weightsPtr_; scalarListList& w = *weightsPtr_;
@ -386,12 +477,20 @@ void topoPatchMapper::calcIntersectionWeightsAndCentres() const
// Fetch interpolative addressing // Fetch interpolative addressing
const labelListList& addr = addressing(); const labelListList& addr = addressing();
// Allocate memory if (addr.empty())
areasPtr_ = new List<scalarField>(patchSize(), scalarField(0)); {
List<scalarField>& a = *areasPtr_; // Nothing to map from. Return empty.
areasPtr_ = new List<scalarList>(0);
centresPtr_ = new List<vectorList>(0);
return;
}
centresPtr_ = new List<vectorField>(patchSize(), vectorField(0)); // Allocate memory
List<vectorField>& x = *centresPtr_; 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 // Obtain stored face-centres
const vectorField& faceCentres = tMapper_.patchCentres(patch_.index()); const vectorField& faceCentres = tMapper_.patchCentres(patch_.index());
@ -418,8 +517,8 @@ void topoPatchMapper::calcIntersectionWeightsAndCentres() const
// Check if this is indeed a mapped face // Check if this is indeed a mapped face
if (mo.size() == 1 && x[faceI].empty() && a[faceI].empty()) if (mo.size() == 1 && x[faceI].empty() && a[faceI].empty())
{ {
x[faceI] = vectorField(1, faceCentres[mo[0]]); x[faceI] = vectorList(1, faceCentres[mo[0]]);
a[faceI] = scalarField(1, 1.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()) if (direct())
{ {
FatalErrorIn FatalErrorIn
( (
"const List<scalarField>& " "const List<scalarList>& "
"topoPatchMapper::intersectionWeights() const" "topoPatchMapper::intersectionWeights() const"
) << "Requested interpolative weights for a direct mapper." ) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError); << 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()) if (direct())
{ {
FatalErrorIn FatalErrorIn
( (
"const List<vectorField>& " "const List<vectorList>& "
"topoPatchMapper::intersectionCentres() const" "topoPatchMapper::intersectionCentres() const"
) << "Requested interpolative weights for a direct mapper." ) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError); << abort(FatalError);
@ -497,6 +596,8 @@ topoPatchMapper::topoPatchMapper
mpm_(mpm), mpm_(mpm),
tMapper_(mapper), tMapper_(mapper),
direct_(false), direct_(false),
sizeBeforeMapping_(0),
conservative_(false),
directAddrPtr_(NULL), directAddrPtr_(NULL),
interpolationAddrPtr_(NULL), interpolationAddrPtr_(NULL),
weightsPtr_(NULL), weightsPtr_(NULL),
@ -506,6 +607,39 @@ topoPatchMapper::topoPatchMapper
areasPtr_(NULL), areasPtr_(NULL),
centresPtr_(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 // Check for the possibility of direct mapping
if (insertedObjects()) if (insertedObjects())
{ {
@ -543,12 +677,7 @@ label topoPatchMapper::size() const
//- Return size before mapping //- Return size before mapping
label topoPatchMapper::sizeBeforeMapping() const label topoPatchMapper::sizeBeforeMapping() const
{ {
if (patch_.type() == "empty") return sizeBeforeMapping_;
{
return 0;
}
return mpm_.oldPatchSizes()[patch_.index()];
} }
@ -616,6 +745,16 @@ const scalarListList& topoPatchMapper::weights() const
<< abort(FatalError); << abort(FatalError);
} }
if (conservative_)
{
if (!areasPtr_ && !centresPtr_)
{
calcIntersectionWeightsAndCentres();
}
return *areasPtr_;
}
if (!weightsPtr_) if (!weightsPtr_)
{ {
calcInverseDistanceWeights(); 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 * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void topoPatchMapper::operator=(const topoPatchMapper& rhs) void topoPatchMapper::operator=(const topoPatchMapper& rhs)

View file

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

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
// ************************************************************************* //

View file

@ -33,7 +33,7 @@ Author
University of Massachusetts Amherst University of Massachusetts Amherst
All rights reserved All rights reserved
\*----------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "topoMapper.H" #include "topoMapper.H"
#include "mapPolyMesh.H" #include "mapPolyMesh.H"
@ -157,11 +157,24 @@ topoSurfaceMapper::topoSurfaceMapper
mpm_(mpm), mpm_(mpm),
tMapper_(mapper), tMapper_(mapper),
direct_(false), direct_(false),
sizeBeforeMapping_(mpm.nOldInternalFaces()),
directAddrPtr_(NULL), directAddrPtr_(NULL),
interpolationAddrPtr_(NULL), interpolationAddrPtr_(NULL),
weightsPtr_(NULL), weightsPtr_(NULL),
insertedFaceLabelsPtr_(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 // Calculate the insertedFaces list
calcInsertedFaceLabels(); calcInsertedFaceLabels();
@ -188,7 +201,7 @@ label topoSurfaceMapper::size() const
//- Return size before mapping //- Return size before mapping
label topoSurfaceMapper::sizeBeforeMapping() const 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 * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void topoSurfaceMapper::operator=(const topoSurfaceMapper& rhs) void topoSurfaceMapper::operator=(const topoSurfaceMapper& rhs)

View file

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

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
// ************************************************************************* //

View file

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

View file

@ -34,8 +34,8 @@ Author
All rights reserved All rights reserved
SourceFiles SourceFiles
meshOpsI.H
meshOps.C meshOps.C
meshOpsTemplates.C
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
@ -43,7 +43,6 @@ SourceFiles
#define meshOps_H #define meshOps_H
#include "Map.H" #include "Map.H"
#include "line.H"
#include "point.H" #include "point.H"
#include "label.H" #include "label.H"
#include "scalar.H" #include "scalar.H"
@ -51,8 +50,10 @@ SourceFiles
#include "cellList.H" #include "cellList.H"
#include "edgeList.H" #include "edgeList.H"
#include "faceList.H" #include "faceList.H"
#include "triangle.H" #include "triPointRef.H"
#include "tetPointRef.H"
#include "vectorField.H" #include "vectorField.H"
#include "linePointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -93,7 +94,7 @@ namespace meshOps
// Utility method to build a hull of cells // Utility method to build a hull of cells
// connected to the edge (for 2D simplical meshes) // connected to the edge (for 2D simplical meshes)
void constructPrismHull inline void constructPrismHull
( (
const label eIndex, const label eIndex,
const UList<face>& faces, const UList<face>& faces,
@ -107,7 +108,7 @@ namespace meshOps
// Utility method to build a hull of cells (and faces) // Utility method to build a hull of cells (and faces)
// around an edge (for 3D simplical meshes) // around an edge (for 3D simplical meshes)
void constructHull inline void constructHull
( (
const label eIndex, const label eIndex,
const UList<face>& faces, const UList<face>& faces,
@ -117,7 +118,7 @@ namespace meshOps
const UList<label>& neighbour, const UList<label>& neighbour,
const UList<labelList>& faceEdges, const UList<labelList>& faceEdges,
const UList<labelList>& edgeFaces, const UList<labelList>& edgeFaces,
const UList<labelList>& edgePoints, const labelList& hullVertices,
labelList& hullEdges, labelList& hullEdges,
labelList& hullFaces, labelList& hullFaces,
labelList& hullCells, labelList& hullCells,
@ -147,12 +148,12 @@ namespace meshOps
( (
const label cIndex, const label cIndex,
const label fIndex, const label fIndex,
const UList<face> faces, const UList<face>& faces,
const UList<cell> cells const UList<cell>& cells
); );
// Given a cell index, find the centroid / volume of a cell. // Given a cell index, find the centroid / volume of a cell.
inline bool cellCentreAndVolume inline void cellCentreAndVolume
( (
const label cIndex, const label cIndex,
const vectorField& points, const vectorField& points,
@ -160,33 +161,28 @@ namespace meshOps
const UList<cell>& cells, const UList<cell>& cells,
const UList<label>& owner, const UList<label>& owner,
vector& centre, vector& centre,
scalar& volume, scalar& volume
bool checkClosed = false
); );
// Determine whether a segment intersects a triangular face // Determine whether a segment intersects a triangular face
template <class T>
inline bool segmentTriFaceIntersection inline bool segmentTriFaceIntersection
( (
const triangle<Vector<T>, const Vector<T>&>& faceToCheck, const triPointRef& faceToCheck,
const line<Vector<T>, const Vector<T>&>& edgeToCheck, const linePointRef& edgeToCheck,
const T& matchTol, vector& intPoint
Vector<T>& intPoint
); );
// Determine whether the particular point lies // Determine whether the particular point lies
// inside the given triangular face // inside the given triangular face
template <class T>
inline bool pointInTriFace inline bool pointInTriFace
( (
const triangle<Vector<T>, const Vector<T>&>& faceToCheck, const triPointRef& faceToCheck,
const Vector<T>& checkPoint, const vector& checkPoint,
const T& matchTol,
bool testCoplanarity bool testCoplanarity
); );
// Dijkstra's algorithm for the shortest path problem // Dijkstra's algorithm for the shortest path problem
bool Dijkstra inline bool Dijkstra
( (
const Map<point>& points, const Map<point>& points,
const Map<edge>& edges, const Map<edge>& edges,
@ -195,17 +191,6 @@ namespace meshOps
Map<label>& pi 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 // Method to insert a label between two labels in a list
inline void insertLabel inline void insertLabel
( (
@ -247,23 +232,75 @@ namespace meshOps
List<Type>& list 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, // Select a list of elements from connectivity,
// and output to a VTK format // and output to a VTK format
void writeVTK inline void writeVTK
( (
const polyMesh& mesh, const polyMesh& mesh,
const word& name, const word& name,
const labelList& cList, const labelList& cList,
const label primitiveType, const label primitiveType,
const UList<point>& meshPoints, const UList<point>& meshPoints,
const UList<edge>& edges, const UList<edge>& edges = List<edge>(0),
const UList<face>& faces, const UList<face>& faces = List<face>(0),
const UList<cell>& cells, const UList<cell>& cells = List<cell>(0),
const UList<label>& owner 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 // Actual routine to write out the VTK file
void writeVTK inline void writeVTK
( (
const polyMesh& mesh, const polyMesh& mesh,
const word& name, const word& name,
@ -274,7 +311,10 @@ namespace meshOps
const labelListList& cpList = labelListList(0), const labelListList& cpList = labelListList(0),
const label primitiveType = 3, const label primitiveType = 3,
const Map<label>& reversePointMap = Map<label>(), 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 } // End namespace meshOps
@ -285,7 +325,9 @@ namespace meshOps
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "meshOpsI.H" #ifdef NoRepository
# include "meshOpsTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -216,8 +216,8 @@ inline label tetApexPoint
( (
const label cIndex, const label cIndex,
const label fIndex, const label fIndex,
const UList<face> faces, const UList<face>& faces,
const UList<cell> cells const UList<cell>& cells
) )
{ {
label apexPoint = -1; label apexPoint = -1;
@ -241,13 +241,13 @@ inline label tetApexPoint
if (!foundApex) if (!foundApex)
{ {
Info << "fIndex: " << fIndex << ":: " << faces[fIndex] << endl; Pout<< "fIndex: " << fIndex << ":: " << faces[fIndex] << endl;
Info << "cIndex: " << cIndex << ":: " << cellToCheck << endl; Pout<< "cIndex: " << cIndex << ":: " << cellToCheck << endl;
forAll(cellToCheck, faceI) forAll(cellToCheck, faceI)
{ {
Info << '\t' << cellToCheck[faceI] << ":: " Pout<< '\t' << cellToCheck[faceI] << ":: "
<< faces[cellToCheck[faceI]] << endl; << faces[cellToCheck[faceI]] << endl;
} }
FatalErrorIn FatalErrorIn
@ -257,8 +257,8 @@ inline label tetApexPoint
"(\n" "(\n"
" const label cIndex,\n" " const label cIndex,\n"
" const label fIndex,\n" " const label fIndex,\n"
" const UList<face> faces,\n" " const UList<face>& faces,\n"
" const UList<cell> cells\n" " const UList<cell>& cells\n"
")\n" ")\n"
) )
<< "Could not find an apex point in cell: " << cIndex << "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. // Given a cell index, find the centroid / volume of a cell.
// - If checking is enabled, return whether the cell is closed // - If checking is enabled, return whether the cell is closed
inline bool cellCentreAndVolume inline void cellCentreAndVolume
( (
const label cIndex, const label cIndex,
const vectorField& points, const vectorField& points,
@ -280,8 +280,7 @@ inline bool cellCentreAndVolume
const UList<cell>& cells, const UList<cell>& cells,
const UList<label>& owner, const UList<label>& owner,
vector& centre, vector& centre,
scalar& volume, scalar& volume
bool checkClosed
) )
{ {
// Reset inputs // Reset inputs
@ -291,18 +290,12 @@ inline bool cellCentreAndVolume
const cell& cellToCheck = cells[cIndex]; const cell& cellToCheck = cells[cIndex];
// Average face-centres to get an estimate centroid // Average face-centres to get an estimate centroid
vector cEst(vector::zero), fCentre(vector::zero), fArea(vector::zero); vector cEst(vector::zero), fCentre(vector::zero);
vector sumClosed(vector::zero), sumMagClosed(vector::zero);
forAll(cellToCheck, faceI) forAll(cellToCheck, faceI)
{ {
const face& checkFace = faces[cellToCheck[faceI]]; const face& checkFace = faces[cellToCheck[faceI]];
if (checkFace.empty())
{
continue;
}
cEst += checkFace.centre(points); cEst += checkFace.centre(points);
} }
@ -312,93 +305,82 @@ inline bool cellCentreAndVolume
{ {
const face& checkFace = faces[cellToCheck[faceI]]; const face& checkFace = faces[cellToCheck[faceI]];
if (checkFace.empty()) if (checkFace.size() == 3)
{ {
continue; tetPointRef tpr
}
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 =
( (
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 tetVol *= -1.0;
return false; }
// 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 centre /= volume + VSMALL;
return true;
} }
// Determine whether a segment intersects a triangular face // Determine whether a segment intersects a triangular face
template <class T>
inline bool segmentTriFaceIntersection inline bool segmentTriFaceIntersection
( (
const triangle<Vector<T>, const Vector<T>&>& faceToCheck, const triPointRef& faceToCheck,
const line<Vector<T>, const Vector<T>&>& edgeToCheck, const linePointRef& edgeToCheck,
const T& matchTol, vector& intPoint
Vector<T>& intPoint
) )
{ {
// Fetch references // Fetch references
const Vector<T>& p1 = edgeToCheck.start(); const vector& p1 = edgeToCheck.start();
const Vector<T>& p2 = edgeToCheck.end(); const vector& p2 = edgeToCheck.end();
const Vector<T>& a = faceToCheck.a(); // Define face-normal
const Vector<T>& b = faceToCheck.b(); vector n = faceToCheck.normal();
const Vector<T>& c = faceToCheck.c();
// Define custom face-normal
Vector<T> n = (T(0.5) * ((b - a)^(c - a)));
n /= mag(n) + VSMALL; n /= mag(n) + VSMALL;
// Compute uValue // Compute uValue
T numerator = n & (a - p1); scalar numerator = n & (faceToCheck.a() - p1);
T denominator = n & (p2 - p1); scalar denominator = n & (p2 - p1);
// Check if the edge is parallel to the face // Check if the edge is parallel to the face
if (mag(denominator) < VSMALL) if (mag(denominator) < VSMALL)
@ -406,17 +388,16 @@ inline bool segmentTriFaceIntersection
return false; return false;
} }
T u = (numerator / denominator); scalar u = (numerator / denominator);
T tolerance = (matchTol * mag(p2 - p1));
// Check for intersection along line. // Check for intersection along line.
if ((u > tolerance) && (u < (pTraits<T>::one - tolerance))) if ((u >= 0.0) && (u <= 1.0))
{ {
// Compute point of intersection // Compute point of intersection
intPoint = p1 + u*(p2 - p1); intPoint = p1 + u*(p2 - p1);
// Also make sure that intPoint lies within the face // Also make sure that intPoint lies within the face
if (pointInTriFace(faceToCheck, intPoint, matchTol, false)) if (pointInTriFace(faceToCheck, intPoint, false))
{ {
return true; return true;
} }
@ -429,49 +410,45 @@ inline bool segmentTriFaceIntersection
// Determine whether the particular point lies // Determine whether the particular point lies
// inside the given triangular face // inside the given triangular face
template <class T>
inline bool pointInTriFace inline bool pointInTriFace
( (
const triangle<Vector<T>, const Vector<T>&>& faceToCheck, const triPointRef& faceToCheck,
const Vector<T>& cP, const vector& cP,
const T& matchTol,
bool testCoplanarity bool testCoplanarity
) )
{ {
// Copy inputs // Copy inputs
const Vector<T>& a = faceToCheck.a(); const vector& a = faceToCheck.a();
const Vector<T>& b = faceToCheck.b(); const vector& b = faceToCheck.b();
const Vector<T>& c = faceToCheck.c(); const vector& c = faceToCheck.c();
const T pointFive(0.5);
// Compute the normal // 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; return false;
} }
if ( ((pointFive * ((c - b)^(cP - b))) & nf) < pTraits<T>::zero) if ( ((0.5 * ((c - b)^(cP - b))) & nf) < 0.0)
{ {
return false; return false;
} }
if ( ((pointFive * ((a - c)^(cP - c))) & nf) < pTraits<T>::zero) if ( ((0.5 * ((a - c)^(cP - c))) & nf) < 0.0)
{ {
return false; return false;
} }
if (testCoplanarity) if (testCoplanarity)
{ {
Vector<T> ftp(a - cP); vector ftp(a - cP);
// Normalize vectors // Normalize vectors
nf /= mag(nf) + VSMALL; nf /= mag(nf) + VSMALL;
ftp /= mag(ftp) + VSMALL; ftp /= mag(ftp) + VSMALL;
if (mag(ftp & nf) > matchTol) if (mag(ftp & nf) > VSMALL)
{ {
return false; return false;
} }
@ -482,53 +459,120 @@ inline bool pointInTriFace
} }
// Method to insert labels in a face, so that // Parallel blocking send
// right-handedness is preserved. inline void pWrite
template <class T>
inline void insertPointLabels
( (
const Vector<T>& refNorm, const label toID,
const Field<Vector<T> >& points, const label& data
const labelHashSet& pLabels,
face& modFace
) )
{ {
// Need to configure a new face. OPstream::write
face newFace(modFace); (
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) OPstream::waitRequests();
{ IPstream::waitRequests();
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;
}
}
} }
// Take over storage
modFace.transfer(newFace);
} }