Adding initial commit of dynamicTopoFvMesh, sans mapping capabilities.

This commit is contained in:
Sandeep Menon 2010-10-17 20:26:26 -04:00
parent 34693bd858
commit eeb78f3f70
42 changed files with 30957 additions and 3 deletions

View file

@ -0,0 +1,115 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
Stack
Description
Thread-safe stack implementation using an internal DynamicList
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
StackI.H
\*---------------------------------------------------------------------------*/
#ifndef Stack_H
#define Stack_H
#include "label.H"
#include "DynamicList.H"
#include "multiThreader.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class Stack Declaration
\*---------------------------------------------------------------------------*/
class Stack
{
// Private data
//- Internal list
DynamicList<label> stack_;
//- Mutex for multi-threading
Mutex stackMutex_;
public:
// Constructor
Stack(){}
// Member functions for access to the stack
//- Push items on to the stack
inline void push(const label index);
//- Insert item onto stack
inline void insert(const label index);
//- Pop an item off the stack
inline label pop();
//- Remove a specific index off the stack
inline void remove(const label index);
//- Return if a stack is empty or not
inline bool empty();
//- Return the size of the stack
inline label size();
//- Clear out the stack
inline void clear();
//- Print out the stack
inline void print();
//- Access the stack as a list
inline label operator[](const label index);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "StackI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,158 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
Stack
Description
Member functions of the Stack class
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#include "ListOps.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Push items on to the stack
inline void Stack::push(const label index)
{
stackMutex_.lock();
if (findIndex(stack_, index) == -1)
{
stack_.append(index);
}
stackMutex_.unlock();
}
//- Insert item onto stack (no checking)
inline void Stack::insert(const label index)
{
stack_.append(index);
}
// Pop an item off the stack
inline label Stack::pop()
{
stackMutex_.lock();
const label index = stack_.remove();
stackMutex_.unlock();
return index;
}
// Remove a specific item off the stack
inline void Stack::remove(const label index)
{
stackMutex_.lock();
label loc = findIndex(stack_, index);
if (loc != -1)
{
// Create a new list
labelList newList(stack_.size() - 1);
label n = 0;
// Copy items upto location
for(label i = 0; i < loc; i++)
{
newList[n++] = stack_[i];
}
// Copy items from location
for(label i = (loc + 1); i < stack_.size(); i++)
{
newList[n++] = stack_[i];
}
// Overwrite
stack_ = newList;
}
stackMutex_.unlock();
}
// Return if the stack is empty or not
inline bool Stack::empty()
{
return (stack_.size() == 0);
}
//- Return the size of the stack
inline label Stack::size()
{
return stack_.size();
}
//- Clear out the stack
inline void Stack::clear()
{
stackMutex_.lock();
stack_.clear();
stackMutex_.unlock();
}
//- Print out the stack
inline void Stack::print()
{
Info << stack_ << endl;
}
//- Access the stack as a list
inline label Stack::operator[]
(
const label index
)
{
return stack_[index];
}
} // End namespace Foam
// ************************************************************************* //

View file

@ -23,4 +23,22 @@ turboFvMesh/turboFvMesh.C
tetMetrics/tetMetric.C
tetMetrics/tetMetrics.C
lengthScaleEstimator/lengthScaleEstimator.C
eMesh = dynamicTopoFvMesh/eMesh
$(eMesh)/eMesh.C
$(eMesh)/eMeshDemandDrivenData.C
$(eMesh)/eBoundaryMesh/eBoundaryMesh.C
ePatches = $(eMesh)/ePatches
$(ePatches)/ePatch/ePatch.C
$(ePatches)/ePatch/newEPatch.C
dynamicTopoFvMesh/meshOps.C
dynamicTopoFvMesh/dynamicTopoFvMesh.C
dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C
dynamicTopoFvMesh/dynamicTopoFvMeshReOrder.C
dynamicTopoFvMesh/dynamicTopoFvMeshMapping.C
dynamicTopoFvMesh/edgeSwap.C
dynamicTopoFvMesh/edgeBisect.C
dynamicTopoFvMesh/edgeCollapse.C
LIB = $(FOAM_LIBBIN)/libdynamicFvMesh

View file

@ -6,7 +6,9 @@ EXE_INC = \
-I$(LIB_SRC)/dynamicMesh/meshMotion/tetDecompositionMotionSolver/lnInclude \
-I$(LIB_SRC)/tetDecompositionFiniteElement/lnInclude \
$(WM_DECOMP_INC) \
-I$(LIB_SRC)/dynamicMesh/meshMotion/fvMotionSolver/lnInclude
-I$(LIB_SRC)/dynamicMesh/meshMotion/fvMotionSolver/lnInclude \
-I$(LIB_SRC)/dynamicMesh/meshMotion/RBFMotionSolver/lnInclude \
-I$(LIB_SRC)/dynamicMesh/meshMotion/mesquiteMotionSolver/lnInclude
LIB_LIBS = \
-ltriSurface \
@ -14,4 +16,6 @@ LIB_LIBS = \
-ldynamicMesh \
-lfiniteVolume \
$(WM_DECOMP_LIBS) \
-lfvMotionSolver
-lfvMotionSolver \
-lRBFMotionSolver \
-lmesquiteMotionSolver

View file

@ -0,0 +1,162 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
changeMap
Description
Accumulate topology change statistics
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
changeMapI.H
\*---------------------------------------------------------------------------*/
#ifndef changeMap_H
#define changeMap_H
#include "labelList.H"
#include "FixedList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class changeMap Declaration
\*---------------------------------------------------------------------------*/
class changeMap
{
// Sliver type index.
// Type 1: Sliver
// Type 2: Cap
// Type 3: Spade
// Type 4: Wedge
label type_;
// Data specific to sliver-type cells
label firstEdge_;
label secondEdge_;
// Data specific to cap-type cells
label apexPoint_;
label opposingFace_;
// Entities that were added during the operation.
// - Add a master mapping for slaves.
// - Masters get a map of -1.
List<FixedList<label,2> > addedPoints_;
List<FixedList<label,2> > addedEdges_;
List<FixedList<label,2> > addedFaces_;
List<FixedList<label,2> > addedCells_;
public:
// Constructor
changeMap()
:
type_(-1),
firstEdge_(-1),
secondEdge_(-1),
apexPoint_(-1),
opposingFace_(-1),
addedPoints_(0),
addedEdges_(0),
addedFaces_(0),
addedCells_(0)
{}
//- Access
// Type
inline label& type();
inline label type() const;
// For sliver-type cells, opposite edges
// are identified for removal.
inline label& firstEdge();
inline label& secondEdge();
// For cap-type cells, the face requiring splitting
// is identified for removal.
inline label& apexPoint();
inline label& opposingFace();
// Added entities
inline void addPoint
(
const label pIndex,
const label master = -1
);
inline void addEdge
(
const label eIndex,
const label master = -1
);
inline void addFace
(
const label fIndex,
const label master = -1
);
inline void addCell
(
const label cIndex,
const label master = -1
);
// Return the list of added entities
inline const List<FixedList<label,2> >& addedPointList() const;
inline const List<FixedList<label,2> >& addedEdgeList() const;
inline const List<FixedList<label,2> >& addedFaceList() const;
inline const List<FixedList<label,2> >& addedCellList() const;
//- Operators
inline void operator=(const changeMap& rhs);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "changeMapI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,181 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 * * * * * * * * * * * * * //
// Type
inline label& changeMap::type()
{
return type_;
}
inline label changeMap::type() const
{
return type_;
}
// For sliver-type cells, opposite edges
// are identified for removal.
inline label& changeMap::firstEdge()
{
return firstEdge_;
}
inline label& changeMap::secondEdge()
{
return secondEdge_;
}
// For cap-type cells, the face requiring splitting
// is identified for removal.
inline label& changeMap::apexPoint()
{
return apexPoint_;
}
inline label& changeMap::opposingFace()
{
return opposingFace_;
}
// Added entities
inline void changeMap::addPoint
(
const label pIndex,
const label master
)
{
label curSize = addedPoints_.size();
addedPoints_.setSize(curSize + 1);
addedPoints_[curSize][0] = pIndex;
addedPoints_[curSize][1] = master;
}
inline void changeMap::addEdge
(
const label eIndex,
const label master
)
{
label curSize = addedEdges_.size();
addedEdges_.setSize(curSize + 1);
addedEdges_[curSize][0] = eIndex;
addedEdges_[curSize][1] = master;
}
inline void changeMap::addFace
(
const label fIndex,
const label master
)
{
label curSize = addedFaces_.size();
addedFaces_.setSize(curSize + 1);
addedFaces_[curSize][0] = fIndex;
addedFaces_[curSize][1] = master;
}
inline void changeMap::addCell
(
const label cIndex,
const label master
)
{
label curSize = addedCells_.size();
addedCells_.setSize(curSize + 1);
addedCells_[curSize][0] = cIndex;
addedCells_[curSize][1] = master;
}
// Return an added point
inline const List<FixedList<label,2> >&
changeMap::addedPointList() const
{
return addedPoints_;
}
// Return the list of added entities
inline const List<FixedList<label,2> >&
changeMap::addedEdgeList() const
{
return addedEdges_;
}
inline const List<FixedList<label,2> >&
changeMap::addedFaceList() const
{
return addedFaces_;
}
inline const List<FixedList<label,2> >&
changeMap::addedCellList() const
{
return addedCells_;
}
inline void changeMap::operator=(const changeMap& rhs)
{
type_ = rhs.type_;
firstEdge_ = rhs.firstEdge_;
secondEdge_ = rhs.secondEdge_;
apexPoint_ = rhs.apexPoint_;
opposingFace_ = rhs.opposingFace_;
addedPoints_ = rhs.addedPoints_;
addedEdges_ = rhs.addedEdges_;
addedFaces_ = rhs.addedFaces_;
addedCells_ = rhs.addedCells_;
}
} // End namespace Foam
// ************************************************************************* //

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,845 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
dynamicTopoFvMesh
Description
An implementation of dynamic changes to mesh-topology
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
dynamicTopoFvMesh.C
dynamicTopoFvMeshI.H
dynamicTopoFvMeshCheck.C
dynamicTopoFvMeshReOrder.C
dynamicTopoFvMeshMapping.C
edgeBisect.C
edgeCollapse.C
edgeSwap.C
\*---------------------------------------------------------------------------*/
#ifndef dynamicTopoFvMesh_H
#define dynamicTopoFvMesh_H
#include "Switch.H"
#include "tetMetric.H"
#include "threadHandler.H"
#include "dynamicFvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
class eMesh;
class Stack;
class changeMap;
class objectMap;
class topoMapper;
class motionSolver;
class lengthScaleEstimator;
/*---------------------------------------------------------------------------*\
Class dynamicTopoFvMesh Declaration
\*---------------------------------------------------------------------------*/
class dynamicTopoFvMesh
:
public dynamicFvMesh
{
// Private data
//- Topology change flag
bool topoChangeFlag_;
//- Dynamic mesh dictionary
IOdictionary dict_;
//- Should all options be mandatorily specified?
// Handy for first-time use.
Switch mandatory_;
//- Mesh characteristics [2D/3D]
Switch twoDMesh_;
//- Edge refinement switch
Switch edgeRefinement_;
//- Switch for cell-bandwidth reduction
Switch bandWidthReduction_;
//- Specify the re-meshing interval
label interval_;
//- Edge-mesh
autoPtr<eMesh> eMeshPtr_;
//- Field mapper
autoPtr<topoMapper> mapper_;
//- Mesh motion solver
autoPtr<motionSolver> motionSolver_;
//- Length scale estimator
autoPtr<lengthScaleEstimator> lengthEstimator_;
//- Lists that dynamically resize during topo-changes
// - Since resizes happen infrequently,
// scale up by 10% to save memory.
// - Current C++ standard doesn't
// support template typedefs,
// so this is a work-around
template <class T>
class resizableList
{
public:
typedef DynamicList<T, 0, 11, 10> Type;
};
resizableList<point>::Type oldPoints_, points_;
resizableList<face>::Type faces_;
resizableList<label>::Type owner_, neighbour_;
resizableList<cell>::Type cells_;
resizableList<edge>::Type edges_;
resizableList<labelList>::Type pointEdges_, edgePoints_;
resizableList<labelList>::Type edgeFaces_, faceEdges_;
resizableList<scalar>::Type lengthScale_;
//- Size information
labelList oldPatchSizes_, patchSizes_;
labelList oldPatchStarts_, patchStarts_;
labelList oldEdgePatchSizes_, edgePatchSizes_;
labelList oldEdgePatchStarts_, edgePatchStarts_;
labelList oldPatchNMeshPoints_, patchNMeshPoints_;
label nOldPoints_, nPoints_;
label nOldEdges_, nEdges_;
label nOldFaces_, nFaces_;
label nOldCells_, nCells_;
label nOldInternalFaces_, nInternalFaces_;
label nOldInternalEdges_, nInternalEdges_;
//- Maps [Old-to-New]
labelList reversePointMap_;
labelList reverseEdgeMap_;
labelList reverseFaceMap_;
labelList reverseCellMap_;
//- Maps [New-to-Old]
labelList pointMap_;
labelList edgeMap_;
labelList faceMap_;
labelList cellMap_;
//- Maps for the renumbering of added entities
Map<label> addedPointRenumbering_;
Map<label> addedEdgeRenumbering_;
Map<label> addedFaceRenumbering_;
Map<label> addedCellRenumbering_;
Map<label> addedFacePatches_;
Map<label> addedEdgePatches_;
Map<label> addedPointZones_;
Map<label> addedFaceZones_;
Map<label> addedCellZones_;
// Information for field-mapping
Map<labelList> faceParents_;
List<scalarField> faceWeights_;
List<vectorField> faceCentres_;
Map<labelList> cellParents_;
List<scalarField> cellWeights_;
List<vectorField> cellCentres_;
// Information for mapPolyMesh
List<objectMap> pointsFromPoints_;
List<objectMap> facesFromPoints_;
List<objectMap> facesFromEdges_;
List<objectMap> facesFromFaces_;
List<objectMap> cellsFromPoints_;
List<objectMap> cellsFromEdges_;
List<objectMap> cellsFromFaces_;
List<objectMap> cellsFromCells_;
//- Maps to keep track of entities deleted after addition
labelHashSet deletedPoints_;
labelHashSet deletedEdges_;
labelHashSet deletedFaces_;
labelHashSet deletedCells_;
//- List of flipped faces / modified old-points
labelHashSet flipFaces_;
Map<labelList> modPoints_;
//- Run-time statistics
label maxModifications_;
FixedList<label, 8> statistics_;
//- Sliver exudation
scalar sliverThreshold_;
Map<scalar> thresholdSlivers_;
//- Specific to proximity-based refinement
List<labelPair> slicePairs_;
//- 3D Mesh Flipping data
label maxTetsPerEdge_;
scalar swapDeviation_;
Switch allowTableResize_;
labelList noSwapPatchIDs_;
//- Stack-list of entities to be checked for topo-changes.
List<Stack> entityStack_;
//- Support for multithreading
autoPtr<IOmultiThreader> threader_;
// Pointer-list of thread-handlers
typedef threadHandler<dynamicTopoFvMesh> meshHandler;
PtrList<meshHandler> handlerPtr_;
// Entity mutexes used for synchronization
// in multi-threaded reOrdering
FixedList<Mutex, 4> entityMutex_;
// Private Member Functions
//- Disallow default bitwise copy construct
dynamicTopoFvMesh(const dynamicTopoFvMesh&);
//- Disallow default bitwise assignment
void operator=(const dynamicTopoFvMesh&);
// Perform a Delaunay test on an internal face
bool testDelaunay(const label fIndex) const;
//- Quality metric for tetrahedra in 3D
tetMetric::tetMetricReturnType tetMetric_;
// Compute mapping weights for modified entities
void computeMapping
(
const scalar matchTol,
const bool skipMapping,
const label faceStart,
const label faceSize,
const label cellStart,
const label cellSize
);
// Static equivalent for multiThreading
static void computeMappingThread(void *argument);
// Routine to invoke threaded mapping
void threadedMapping(scalar matchTol, bool skipMapping);
// Initialize mesh edges and related connectivity lists
void initEdges();
// Load the mesh-quality metric from the library
void loadMetric();
// Load the mesh-motion solver
void loadMotionSolver();
// Load the field mapper
void loadFieldMapper();
// Load the length scale estimator
void loadLengthScaleEstimator();
// Initialize the threading environment
void initializeThreadingEnvironment(const label specThreads = -1);
// Return a non-const reference to the lengthScaleEstimator
inline lengthScaleEstimator& lengthEstimator();
// Return a const reference to the lengthScaleEstimator
inline const lengthScaleEstimator& lengthEstimator() const;
// Return a const reference to the multiThreader
inline const multiThreader& threader() const;
// Return a reference to the entity mutexes
inline const Mutex& entityMutex(const label entity) const;
// Return the edge index for a provided edge
inline label getEdgeIndex(const edge& edgeToCheck) const;
// Given a boundary face, pick out a boundary edge that
// contains a triangular face.
inline label getTriBoundaryEdge(const label fIndex) const;
// Return length-scale at an face-location in the mesh [2D]
inline scalar faceLengthScale(const label fIndex) const;
// Return length-scale at an edge-location in the mesh [3D]
inline scalar edgeLengthScale(const label eIndex) const;
// Check for bisection
inline bool checkBisection(const label index) const;
// Check for collapse
inline bool checkCollapse(const label index) const;
// MultiThreaded topology modifier
void threadedTopoModifier();
// 2D Edge-swapping engine
static void swap2DEdges(void *argument);
// 3D Edge-swapping engine
static void swap3DEdges(void *argument);
// Edge refinement engine
static void edgeRefinementEngine(void *argument);
// Return the entity stack for a particular thread
inline Stack& stack(const label threadID);
// Return the integer ID for a given thread
inline label self() const;
// Initialize stacks
inline void initStacks(const labelHashSet& entities);
// Method to determine the boundary patch index for a given face
inline label whichPatch(const label index) const;
// Method to determine the boundary patch index for a given edge
inline label whichEdgePatch(const label index) const;
// Report topo-modification status
inline label status(const label type) const;
// Method for the swapping of a quad-face in 2D
const changeMap
swapQuadFace
(
const label fIndex
);
// Method for the bisection of a quad-face in 2D
const changeMap
bisectQuadFace
(
const label fIndex,
const changeMap& masterMap,
bool checkOnly = false,
bool forceOp = false
);
// Method for the collapse of a quad-face in 2D
const changeMap
collapseQuadFace
(
const label fIndex,
label overRideCase = -1,
bool checkOnly = 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
const changeMap
bisectEdge
(
const label eIndex,
bool checkOnly = false,
bool forceOp = false
);
// Method for the trisection of a face in 3D
const changeMap
trisectFace
(
const label fIndex,
bool checkOnly = false,
bool forceOp = false
);
// Method for the collapse of an edge in 3D
const changeMap
collapseEdge
(
const label eIndex,
label overRideCase = -1,
bool checkOnly = false,
bool forceOp = false
);
// Utility method to check for invalid face-bisection.
bool checkBisection
(
const label fIndex,
const label bFaceIndex,
bool forceOp = false
) const;
// Utility method to check for invalid face-collapse.
bool checkCollapse
(
const labelList& triFaces,
const FixedList<label,2>& c0BdyIndex,
const FixedList<label,2>& c1BdyIndex,
const FixedList<label,2>& pointIndex,
const FixedList<point,2>& newPoint,
const FixedList<point,2>& oldPoint,
scalar& collapseQuality,
const bool checkNeighbour,
bool forceOp = false
) const;
// Utility method to check for invalid edge-collapse.
bool checkCollapse
(
const point& newPoint,
const point& oldPoint,
const label pointIndex,
const label cellIndex,
labelHashSet& cellsChecked,
scalar& collapseQuality,
bool forceOp = false
) const;
// Remove 2D sliver cells from the mesh
void remove2DSlivers();
// Identify the sliver type in 3D
const changeMap identifySliverType(const label cIndex) const;
// Remove sliver cells
void removeSlivers();
// Insert the specified cell to the mesh
label insertCell
(
const cell& newCell,
const scalar lengthScale,
const label zoneID = -1
);
// Set fill-in mapping information for a particular cell
void setCellMapping
(
const label cIndex,
const labelList& mapCells,
bool addEntry = true
);
// Set fill-in mapping information for a particular face
void setFaceMapping
(
const label fIndex,
const labelList& mapFaces = labelList(0)
);
// Remove the specified cells from the mesh
const changeMap
removeCells
(
const labelList& cList,
const label patch
);
// Remove the specified cell from the mesh
void removeCell
(
const label cIndex
);
// Insert the specified face to the mesh
label insertFace
(
const label patch,
const face& newFace,
const label newOwner,
const label newNeighbour,
const label zoneID = -1
);
// Remove the specified face from the mesh
void removeFace
(
const label fIndex
);
// Split a set of internal faces into boundary faces,
// and add them to the specified patch.
void splitInternalFaces
(
const label patchIndex,
const labelList& internalFaces,
const Map<bool>& cellColors
);
// Insert the specified edge to the mesh
label insertEdge
(
const label patch,
const edge& newEdge,
const labelList& edgeFaces,
const labelList& edgePoints = labelList::null()
);
// Remove the specified edge from the mesh
void removeEdge
(
const label eIndex
);
// Insert the specified point to the mesh
label insertPoint
(
const point& newPoint,
const point& oldPoint,
const labelList& mapPoints,
const label zoneID = -1
);
// Remove the specified point from the mesh
void removePoint
(
const label pIndex
);
// Utility method to build edgePoints for an edge [3D]
void buildEdgePoints
(
const label eIndex,
const label checkIndex = 0
);
// Utility to check whether points of an edge lie on a boundary
const FixedList<bool,2> checkEdgeBoundary(const label eIndex) const;
// Set a particular face index as flipped.
inline void setFlip(const label fIndex);
// Utility method to compute the minimum quality of a vertex hull
scalar computeMinQuality(const label eIndex) const;
// Utility method to compute the quality of a vertex hull
// around an edge after bisection.
scalar computeBisectionQuality(const label eIndex) const;
// Utility method to compute the quality of cells
// around a face after trisection.
scalar computeTrisectionQuality(const label fIndex) const;
// Check whether the given edge is on a bounding curve
bool checkBoundingCurve(const label eIndex) const;
// Allocate dynamic programming tables
void initTables
(
labelList& m,
PtrList<scalarListList>& Q,
PtrList<labelListList>& K,
PtrList<labelListList>& triangulations,
const label checkIndex = -1
) const;
// Check triangulation quality for an edge index
bool checkQuality
(
const label eIndex,
const labelList& m,
const PtrList<scalarListList>& Q,
const scalar minQuality,
const label checkIndex = 0
) const;
// Utility method to fill the dynamic programming tables
bool fillTables
(
const label eIndex,
const scalar minQuality,
labelList& m,
PtrList<scalarListList>& Q,
PtrList<labelListList>& K,
PtrList<labelListList>& triangulations,
const label checkIndex = 0
) const;
// Print out tables for debugging
void printTables
(
const labelList& m,
const PtrList<scalarListList>& Q,
const PtrList<labelListList>& K,
const label checkIndex = 0
) const;
// Remove the edge according to the swap sequence
const changeMap
removeEdgeFlips
(
const label eIndex,
const scalar minQuality,
const PtrList<labelListList>& K,
PtrList<labelListList>& triangulations,
const label checkIndex = 0
);
// Extract triangulations from the programming table
void extractTriangulation
(
const label i,
const label j,
const labelListList& K,
label& numTriangulations,
labelListList& triangulations
) const;
// Identify the 3-2 swap from the triangulation sequence
label identify32Swap
(
const label eIndex,
const labelList& hullVertices,
const labelListList& triangulations
) const;
// Check old-volumes for the configuration.
bool checkTriangulationVolumes
(
const label eIndex,
const labelList& hullVertices,
const labelListList& triangulations
) const;
// Routine to check whether the triangulation at the
// index lies on the boundary of the vertex ring.
bool boundaryTriangulation
(
const label index,
label& isolatedVertex,
labelListList& triangulations
) const;
// Routine to perform 2-3 swaps
const changeMap
swap23
(
const label isolatedVertex,
const label eIndex,
const label triangulationIndex,
const label numTriangulations,
const labelListList& triangulations,
const labelList& hullVertices,
const labelList& hullFaces,
const labelList& hullCells
);
// Routine to perform 3-2 or 2-2 swaps
const changeMap
swap32
(
const label eIndex,
const label triangulationIndex,
const label numTriangulations,
const labelListList& triangulations,
const labelList& hullVertices,
const labelList& hullFaces,
const labelList& hullCells
);
// Reorder points after a topology change
void reOrderPoints
(
pointField& points,
pointField& preMotionPoints,
labelListList& pointZoneMap,
bool threaded = false
);
// Static equivalent for multi-threading
static void reOrderPointsThread(void *argument);
// Reorder edges after a topology change
void reOrderEdges
(
edgeList& edges,
labelListList& edgeFaces,
labelListList& faceEdges,
bool threaded = false
);
// Static equivalent for multi-threading
static void reOrderEdgesThread(void *argument);
// Reorder faces in upper-triangular order after a topology change
void reOrderFaces
(
faceList& faces,
labelList& owner,
labelList& neighbour,
labelListList& faceEdges,
labelListList& faceZoneFaceMap,
bool threaded = false
);
// Static equivalent for multi-threading
static void reOrderFacesThread(void *argument);
// Reorder & renumber cells with bandwidth
// reduction after a topology change
void reOrderCells
(
labelListList& cellZoneMap,
bool threaded = false
);
// Static equivalent for multi-threading
static void reOrderCellsThread(void *argument);
// Reorder the mesh in upper-triangular order,
// and generate mapping information
void reOrderMesh
(
pointField& points,
pointField& preMotionPoints,
edgeList& edges,
faceList& faces,
labelList& owner,
labelList& neighbour,
labelListList& faceEdges,
labelListList& edgeFaces,
labelListList& pointZoneMap,
labelListList& faceZoneFaceMap,
labelListList& cellZoneMap
);
// Invoke reOrdering with multiple threads
void threadedMeshReOrdering
(
pointField& points,
pointField& preMotionPoints,
edgeList& edges,
faceList& faces,
labelList& owner,
labelList& neighbour,
labelListList& faceEdges,
labelListList& edgeFaces,
labelListList& pointZoneMap,
labelListList& faceZoneFaceMap,
labelListList& cellZoneMap
);
// Reset the mesh and generate mapping information
bool resetMesh();
// Output an entity as a VTK file
void writeVTK
(
const word& name,
const label entity,
const label primitiveType = 3,
const bool useOldConnectivity = false,
const bool useOldPoints = 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 bool useOldPoints = false
) const;
// Return the status report interval
scalar reportInterval() const;
// Check the state of local connectivity lists
void checkConnectivity(const label maxErrors = 0) const;
// Test an edge / face for proximity with other non-neighbouring faces.
// Return the scalar distance to the nearest face.
scalar testProximity
(
const label index,
label& proximityFace
) const;
// Calculate the edge length-scale for the mesh
void calculateLengthScale(bool dump = false);
// Read optional dictionary parameters
void readOptionalParameters(bool reRead = false);
// Dump cell-quality statistics
bool meshQuality(bool outputOption);
public:
// Declare the name of the class and its debug switch
TypeName("dynamicTopoFvMesh");
// Constructors
//- Construct from existing IOobject
explicit dynamicTopoFvMesh(const IOobject& io);
// Destructor
virtual ~dynamicTopoFvMesh();
// Member Functions
// Map all fields in time using a customized mapper
virtual void mapFields(const mapPolyMesh& meshMap) const;
// Update the mesh for motion / topology changes
// Return true if topology changes have occurred
virtual bool update();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "dynamicTopoFvMeshI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,823 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
dynamicTopoFvMesh
Description
An implementation of dynamic changes to mesh-topology
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#include "Stack.H"
#include "meshOps.H"
#include "tetPointRef.H"
#include "linePointRef.H"
#include "lengthScaleEstimator.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Return a non-const reference to the lengthScaleEstimator
inline lengthScaleEstimator&
dynamicTopoFvMesh::lengthEstimator()
{
if (!lengthEstimator_.valid())
{
FatalErrorIn
(
"inline lengthScaleEstimator& "
"dynamicTopoFvMesh::lengthEstimator()"
) << nl
<< " Invalid request. Length scale estimator was not loaded. "
<< abort(FatalError);
}
return lengthEstimator_();
}
// Return a const reference to the lengthScaleEstimator
inline const lengthScaleEstimator&
dynamicTopoFvMesh::lengthEstimator() const
{
if (!lengthEstimator_.valid())
{
FatalErrorIn
(
"inline const lengthScaleEstimator& "
"dynamicTopoFvMesh::lengthEstimator() const"
) << nl
<< " Invalid request. Length scale estimator was not loaded. "
<< abort(FatalError);
}
return lengthEstimator_();
}
// Return a const reference to the multiThreader
inline const multiThreader&
dynamicTopoFvMesh::threader() const
{
if (!threader_.valid())
{
FatalErrorIn
(
"inline const multiThreader& "
"dynamicTopoFvMesh::threader() const"
) << nl
<< " Invalid request. multiThreader was not loaded. "
<< abort(FatalError);
}
return threader_();
}
// Return a reference to the entity mutexes.
// The index 'entity' ranges from 0 to 3 for point/edge/face/cell.
inline const Mutex& dynamicTopoFvMesh::entityMutex
(
const label entity
) const
{
return entityMutex_[entity];
}
// Return the edge index for a provided edge
inline label dynamicTopoFvMesh::getEdgeIndex
(
const edge& edgeToCheck
) const
{
if (twoDMesh_)
{
// No efficient search method. Loop through all edges.
forAll(edges_, edgeI)
{
if (edges_[edgeI] == edgeToCheck)
{
return edgeI;
}
}
}
else
{
// Look througg pointEdges list
const labelList& pEdges = pointEdges_[edgeToCheck.start()];
forAll(pEdges, edgeI)
{
if (edges_[pEdges[edgeI]] == edgeToCheck)
{
return pEdges[edgeI];
}
}
}
// Probably not a valid edge.
FatalErrorIn
(
"inline label dynamicTopoFvMesh::getEdgeIndex"
"(const edge& edgeToCheck) const"
)
<< "Could not find an appropriate edge index for edge:"
<< edgeToCheck
<< abort(FatalError);
return -1;
}
// Given any quad face, pick out a boundary edge that
// contains a triangular face. For 2D simplical meshes only.
inline label dynamicTopoFvMesh::getTriBoundaryEdge
(
const label fIndex
) const
{
const labelList& fEdges = faceEdges_[fIndex];
forAll(fEdges, edgeI)
{
// Obtain edgeFaces for this edge
const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
forAll(eFaces, faceI)
{
if (faces_[eFaces[faceI]].size() == 3)
{
// Found a triangular face. Return this edge.
return fEdges[edgeI];
}
}
}
// This bit should never happen.
FatalErrorIn
(
"inline label dynamicTopoFvMesh::getTriBoundaryEdge"
"(const label fIndex) const"
)
<< "Cannot find a triangular face bordering face: "
<< fIndex << " :: " << faces_[fIndex]
<< abort(FatalError);
return -1;
}
// 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.
meshOps::faceNormal
(
faces_[eFaces[faceI]],
points_,
fNorm[count]
);
// Normalize it.
fNorm[count] /= mag(fNorm[count]);
count++;
}
}
scalar deviation = (fNorm[0] & fNorm[1]);
scalar refDeviation = lengthEstimator().curvatureDeviation();
if (mag(deviation) < refDeviation)
{
// Fetch the edge
const edge& edgeToCheck = edges_[eIndex];
// Get the edge-length.
scalar length =
(
linePointRef
(
points_[edgeToCheck.start()],
points_[edgeToCheck.end()]
).mag()
);
if (debug > 3 && self() == 0)
{
Info << "Deviation: " << deviation << nl
<< "curvatureDeviation: " << refDeviation
<< ", Edge: " << eIndex << ", Length: " << length
<< ", Scale: " << scale << nl
<< " Half-length: " << (0.5*length) << nl
<< " MinRatio: "
<< (lengthEstimator().ratioMin()*scale)
<< endl;
}
scale =
(
Foam::min
(
scale,
((length - SMALL)/lengthEstimator().ratioMax())
)
);
}
}
// Limit scales if necessary
lengthEstimator().limitScale(scale);
}
return scale;
}
// Check for edge bisection
inline bool dynamicTopoFvMesh::checkBisection
(
const label index
) const
{
scalar length = 0.0, scale = 0.0;
if (twoDMesh_)
{
// Fetch the edge
const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
// Measure the boundary edge-length of the face in question
length =
(
linePointRef
(
points_[edgeToCheck.start()],
points_[edgeToCheck.end()]
).mag()
);
// Determine the length-scale at this face
scale = faceLengthScale(index);
}
else
{
// Fetch the edge
const edge& edgeToCheck = edges_[index];
// Measure the edge-length
length =
(
linePointRef
(
points_[edgeToCheck.start()],
points_[edgeToCheck.end()]
).mag()
);
// Determine the length-scale at this point in the mesh
scale = edgeLengthScale(index);
}
if (length > lengthEstimator().ratioMax() * scale)
{
return true;
}
else
{
return false;
}
}
// Check for edge collapse
inline bool dynamicTopoFvMesh::checkCollapse
(
const label index
) const
{
scalar length = 0.0, scale = 0.0;
if (twoDMesh_)
{
// Fetch the edge
const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
// Measure the boundary edge-length of the face in question
length =
(
linePointRef
(
points_[edgeToCheck.start()],
points_[edgeToCheck.end()]
).mag()
);
// Determine the length-scale at this face
scale = faceLengthScale(index);
}
else
{
// Fetch the edge
const edge& edgeToCheck = edges_[index];
// Measure the edge-length
length =
(
linePointRef
(
points_[edgeToCheck.start()],
points_[edgeToCheck.end()]
).mag()
);
// Determine the length-scale at this point in the mesh
scale = edgeLengthScale(index);
}
if (length < lengthEstimator().ratioMin() * scale)
{
return true;
}
else
{
return false;
}
}
// Return the entity stack
inline Stack& dynamicTopoFvMesh::stack
(
const label threadID
)
{
return entityStack_[threadID];
}
// Return the integer ID for a given thread
// Return zero for single-threaded operation
inline label dynamicTopoFvMesh::self() const
{
if (threader_->multiThreaded())
{
for (label i = 1; i <= threader_->getNumThreads(); i++)
{
if (handlerPtr_[i].self())
{
return i;
}
}
}
return 0;
}
// Initialize edge-stacks
inline void dynamicTopoFvMesh::initStacks
(
const labelHashSet& entities
)
{
forAll(entityStack_, stackI)
{
entityStack_[stackI].clear();
}
// Prepare a filling sequence based on threading operation
label tIndex = 0;
labelList tID(threader().getNumThreads());
if (threader_->multiThreaded())
{
forAll(tID, tI)
{
tID[tI] = (tI + 1);
}
}
else
{
tID = 0;
}
if (twoDMesh_)
{
forAll(faces_, faceI)
{
if (faces_[faceI].size() == 4)
{
stack(tID[tIndex]).insert(faceI);
tIndex = tID.fcIndex(tIndex);
}
}
}
else
{
forAll(edges_, edgeI)
{
if (edgeFaces_[edgeI].size())
{
stack(tID[tIndex]).insert(edgeI);
tIndex = tID.fcIndex(tIndex);
}
}
}
}
// Method to determine the old boundary patch index for a given face
// Similar to the polyBoundaryMesh routine, but works on local information
inline label dynamicTopoFvMesh::whichPatch
(
const label index
) const
{
if (index < nOldInternalFaces_)
{
return -1;
}
for (label i = 0; i < boundaryMesh().size(); i++)
{
if
(
index >= oldPatchStarts_[i]
&& index < oldPatchStarts_[i] + oldPatchSizes_[i]
)
{
return i;
}
}
// If not in any of the above, it's possible that the face was added
// at the end of the list. Check addedFacePatches_ for the patch info
if (addedFacePatches_.found(index))
{
return addedFacePatches_[index];
}
else
{
FatalErrorIn
(
"label dynamicTopoFvMesh::whichPatch"
"(const label index) const"
)
<< "Cannot find patch information for face index: "
<< index << nl
<< " It appears that face ordering is"
<< " inconsistent with patch information." << nl
<< " Mesh info: " << nl
<< " nOldInternalFaces: " << nOldInternalFaces_ << nl
<< " oldPatchStarts: " << oldPatchStarts_ << nl
<< " oldPatchSizes: " << oldPatchSizes_ << nl
<< abort(FatalError);
}
return -2;
}
// Method to determine the old boundary patch index for a given edge
inline label dynamicTopoFvMesh::whichEdgePatch
(
const label index
) const
{
if (index < nOldInternalEdges_)
{
return -1;
}
for (label i = 0; i < boundaryMesh().size(); i++)
{
if
(
index >= oldEdgePatchStarts_[i]
&& index < oldEdgePatchStarts_[i] + oldEdgePatchSizes_[i]
)
{
return i;
}
}
// If not in any of the above, it's possible that the edge was added
// at the end of the list. Check addedEdgePatches_ for the patch info
if (addedEdgePatches_.found(index))
{
return addedEdgePatches_[index];
}
else
{
FatalErrorIn
(
"label dynamicTopoFvMesh::whichEdgePatch"
"(const label index) const"
)
<< "Cannot find patch information for edge index: "
<< index << nl
<< " It appears that edge ordering is"
<< " inconsistent with patch information." << nl
<< " Mesh info: " << nl
<< " nOldInternalEdges: " << nOldInternalEdges_ << nl
<< " oldEdgePatchStarts: " << oldEdgePatchStarts_ << nl
<< " oldEdgePatchSizes: " << oldEdgePatchSizes_ << nl
<< abort(FatalError);
}
return -2;
}
// Report topo-modification status
// - Enumerants are as follows
// [0] nModifications
// [1] Total swaps (internal and surface)
// [2] Surface swaps
// [3] Total bisections
// [4] Total collapses
// [5] Surface bisections
// [6] Surface collapses
// [7] Slivers removed
inline label dynamicTopoFvMesh::status(const label type) const
{
if (type < statistics_.size())
{
return statistics_[type];
}
else
{
FatalErrorIn
(
"inline label dynamicTopoFvMesh::status"
"(const label type) const"
)
<< " Unknown type: " << type << nl
<< abort(FatalError);
}
return 0;
}
// Set a particular face index as flipped.
inline void dynamicTopoFvMesh::setFlip(const label fIndex)
{
if (fIndex < nOldFaces_)
{
if (flipFaces_.found(fIndex))
{
flipFaces_.erase(fIndex);
}
else
{
flipFaces_.insert(fIndex);
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,405 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
dynamicTopoFvMesh
Description
Functions specific to conservative mapping
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
#include "dynamicTopoFvMesh.H"
#include "meshOps.H"
#include "IOmanip.H"
#include "triFace.H"
#include "objectMap.H"
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Compute mapping weights for modified entities
void dynamicTopoFvMesh::computeMapping
(
const scalar matchTol,
const bool skipMapping,
const label faceStart,
const label faceSize,
const label cellStart,
const label cellSize
)
{
}
// Static equivalent for multiThreading
void dynamicTopoFvMesh::computeMappingThread(void *argument)
{
// Recast the argument
meshHandler *thread = static_cast<meshHandler*>(argument);
if (thread->slave())
{
thread->sendSignal(meshHandler::START);
}
dynamicTopoFvMesh& mesh = thread->reference();
// Recast the pointers for the argument
scalar& matchTol = *(static_cast<scalar*>(thread->operator()(0)));
bool& skipMapping = *(static_cast<bool*>(thread->operator()(1)));
label& faceStart = *(static_cast<label*>(thread->operator()(2)));
label& faceSize = *(static_cast<label*>(thread->operator()(3)));
label& cellStart = *(static_cast<label*>(thread->operator()(4)));
label& cellSize = *(static_cast<label*>(thread->operator()(5)));
// Now calculate addressing
mesh.computeMapping
(
matchTol,
skipMapping,
faceStart, faceSize,
cellStart, cellSize
);
if (thread->slave())
{
thread->sendSignal(meshHandler::STOP);
}
}
// Routine to invoke threaded mapping
void dynamicTopoFvMesh::threadedMapping
(
scalar matchTol,
bool skipMapping
)
{
label nThreads = threader_->getNumThreads();
// If mapping is being skipped, issue a warning.
if (skipMapping)
{
Info << " *** Mapping is being skipped *** " << endl;
}
// Check if single-threaded
if (nThreads == 1)
{
computeMapping
(
matchTol,
skipMapping,
0, facesFromFaces_.size(),
0, cellsFromCells_.size()
);
return;
}
// Set one handler per thread
PtrList<meshHandler> hdl(nThreads);
forAll(hdl, i)
{
hdl.set(i, new meshHandler(*this, threader()));
}
// Simple load-balancing scheme
FixedList<label, 2> index(-1);
FixedList<labelList, 2> tStarts(labelList(nThreads, 0));
FixedList<labelList, 2> tSizes(labelList(nThreads, 0));
index[0] = facesFromFaces_.size();
index[1] = cellsFromCells_.size();
if (debug > 2)
{
Info << " Mapping Faces: " << index[0] << endl;
Info << " Mapping Cells: " << index[1] << endl;
}
forAll(index, indexI)
{
label j = 0, total = 0;
while (index[indexI]--)
{
tSizes[indexI][(j = tSizes[indexI].fcIndex(j))]++;
}
for (label i = 1; i < tStarts[indexI].size(); i++)
{
tStarts[indexI][i] = tSizes[indexI][i-1] + total;
total += tSizes[indexI][i-1];
}
if (debug > 2)
{
Info << " Load starts: " << tStarts[indexI] << endl;
Info << " Load sizes: " << tSizes[indexI] << endl;
}
}
// Set the argument list for each thread
forAll(hdl, i)
{
// Size up the argument list
hdl[i].setSize(6);
// Set match tolerance
hdl[i].set(0, &matchTol);
// Set the skipMapping flag
hdl[i].set(1, &skipMapping);
// Set the start/size indices
hdl[i].set(2, &(tStarts[0][i]));
hdl[i].set(3, &(tSizes[0][i]));
hdl[i].set(4, &(tStarts[1][i]));
hdl[i].set(5, &(tSizes[1][i]));
}
// Prior to multi-threaded operation,
// force calculation of demand-driven data.
polyMesh::cells();
primitiveMesh::cellCells();
const polyBoundaryMesh& boundary = boundaryMesh();
forAll(boundary, patchI)
{
boundary[patchI].faceFaces();
}
// Execute threads in linear sequence
executeThreads(identity(nThreads), hdl, &computeMappingThread);
}
// Set fill-in mapping information for a particular cell
void dynamicTopoFvMesh::setCellMapping
(
const label cIndex,
const labelList& mapCells,
bool addEntry
)
{
if (addEntry)
{
if (debug > 3)
{
Info << "Inserting mapping cell: " << cIndex << nl
<< " mapCells: " << mapCells
<< endl;
}
// Insert index into the list, and overwrite if necessary
label index = -1;
forAll(cellsFromCells_, indexI)
{
if (cellsFromCells_[indexI].index() == cIndex)
{
index = indexI;
break;
}
}
if (index == -1)
{
meshOps::sizeUpList
(
objectMap(cIndex, labelList(0)),
cellsFromCells_
);
}
else
{
cellsFromCells_[index].masterObjects() = labelList(0);
}
}
// Update cell-parents information
labelHashSet masterCells;
forAll(mapCells, cellI)
{
if (mapCells[cellI] < 0)
{
continue;
}
if (mapCells[cellI] < nOldCells_)
{
masterCells.insert(mapCells[cellI]);
}
else
if (cellParents_.found(mapCells[cellI]))
{
const labelList& nParents = cellParents_[mapCells[cellI]];
forAll(nParents, cI)
{
masterCells.insert(nParents[cI]);
}
}
}
cellParents_.set(cIndex, masterCells.toc());
}
// Set fill-in mapping information for a particular face
void dynamicTopoFvMesh::setFaceMapping
(
const label fIndex,
const labelList& mapFaces
)
{
label patch = whichPatch(fIndex);
if (debug > 3)
{
Info << "Inserting mapping face: " << fIndex << nl
<< " patch: " << patch << nl
<< " mapFaces: " << mapFaces
<< endl;
}
bool foundError = false;
// Check to ensure that internal faces are not mapped
// from any faces in the mesh
if (patch == -1 && mapFaces.size())
{
foundError = true;
}
// Check to ensure that boundary faces map
// only from other faces on the same patch
if (patch > -1 && mapFaces.empty())
{
foundError = true;
}
if (foundError)
{
FatalErrorIn
(
"\n"
"void dynamicTopoFvMesh::setFaceMapping\n"
"(\n"
" const label fIndex,\n"
" const labelList& mapFaces\n"
")"
)
<< nl << " Incompatible mapping. " << nl
<< " Possible reasons: " << nl
<< " 1. No mapping specified for a boundary face; " << nl
<< " 2. Mapping specified for an internal face, " << nl
<< " when none was expected." << nl << nl
<< " Face: " << fIndex << nl
<< " Patch: " << patch << nl
<< " Owner: " << owner_[fIndex] << nl
<< " Neighbour: " << neighbour_[fIndex] << nl
<< " mapFaces: " << mapFaces << nl
<< abort(FatalError);
}
// Insert addressing into the list, and overwrite if necessary
label index = -1;
forAll(facesFromFaces_, indexI)
{
if (facesFromFaces_[indexI].index() == fIndex)
{
index = indexI;
break;
}
}
if (index == -1)
{
meshOps::sizeUpList
(
objectMap(fIndex, labelList(0)),
facesFromFaces_
);
}
else
{
facesFromFaces_[index].masterObjects() = labelList(0);
}
// For internal faces, set dummy maps / weights, and bail out
if (patch == -1)
{
return;
}
// Update face-parents information
labelHashSet masterFaces;
forAll(mapFaces, faceI)
{
if (mapFaces[faceI] < 0)
{
continue;
}
if (mapFaces[faceI] < nOldFaces_)
{
masterFaces.insert(mapFaces[faceI]);
}
else
if (faceParents_.found(mapFaces[faceI]))
{
const labelList& nParents = faceParents_[mapFaces[faceI]];
forAll(nParents, fI)
{
masterFaces.insert(nParents[fI]);
}
}
}
faceParents_.set(fIndex, masterFaces.toc());
}
} // End namespace Foam
// ************************************************************************* //

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,295 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
Description
\*---------------------------------------------------------------------------*/
#include "eBoundaryMesh.H"
#include "eMesh.H"
#include "primitiveMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(eBoundaryMesh, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from dictionary
eBoundaryMesh::eBoundaryMesh
(
const IOobject& io,
const eMesh& mesh
)
:
ePatchList(),
regIOobject(io),
mesh_(mesh)
{
if
(
(readOpt() == IOobject::MUST_READ) ||
(readOpt() == IOobject::READ_IF_PRESENT && headerOk())
)
{
readFromInputStream();
}
}
// Construct given size. Patches will be set later
eBoundaryMesh::eBoundaryMesh
(
const IOobject& io,
const eMesh& em,
const label size
)
:
ePatchList(size),
regIOobject(io),
mesh_(em)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Read from input stream
void eBoundaryMesh::readFromInputStream()
{
ePatchList& patches = *this;
// Read polyPatchList
Istream& is = readStream(typeName);
PtrList<entry> patchEntries(is);
patches.setSize(patchEntries.size());
forAll(patches, patchI)
{
patches.set
(
patchI,
ePatch::New
(
patchEntries[patchI].keyword(),
patchEntries[patchI].dict(),
patchI,
*this
)
);
}
// Check state of IOstream
is.check
(
"eBoundaryMesh::eBoundaryMesh"
"(const IOobject&, const faMesh&)"
);
close();
}
// Return the mesh reference
const eMesh& eBoundaryMesh::mesh() const
{
return mesh_;
}
// Return a list of patch types
wordList eBoundaryMesh::types() const
{
const ePatchList& patches = *this;
wordList t(patches.size());
forAll (patches, patchI)
{
t[patchI] = patches[patchI].type();
}
return t;
}
// Return a list of patch names
wordList eBoundaryMesh::names() const
{
const ePatchList& patches = *this;
wordList t(patches.size());
forAll (patches, patchI)
{
t[patchI] = patches[patchI].name();
}
return t;
}
// Return a list of patch sizes
labelList eBoundaryMesh::patchSizes() const
{
const ePatchList& patches = *this;
labelList t(patches.size());
forAll (patches, patchI)
{
t[patchI] = patches[patchI].size();
}
return t;
}
// Return a list of patch starts
labelList eBoundaryMesh::patchStarts() const
{
const ePatchList& patches = *this;
labelList t(patches.size());
forAll (patches, patchI)
{
t[patchI] = patches[patchI].start();
}
return t;
}
//- Return patch index for a given edge label
label eBoundaryMesh::whichPatch(const label edgeIndex) const
{
// Find out which patch the current edge belongs to by comparing label
// with patch start labels.
// If the face is internal, return -1;
// if it is off the end of the list, abort
if (edgeIndex >= mesh().nEdges())
{
FatalErrorIn
(
"eBoundaryMesh::whichPatch(const label edgeIndex) const"
) << "given label greater than the number of geometric edges"
<< abort(FatalError);
}
if (edgeIndex < mesh().nInternalEdges())
{
return -1;
}
forAll (*this, patchI)
{
const ePatch& bp = operator[](patchI);
if
(
edgeIndex >= bp.start()
&& edgeIndex < bp.start() + bp.size()
)
{
return patchI;
}
}
// If not in any of above, it is trouble!
FatalErrorIn
(
"label eBoundaryMesh::whichPatch(const label edgeIndex) const"
) << "Cannot find edge " << edgeIndex << " in any of the patches "
<< names() << nl
<< "It seems your patches are not consistent with the mesh :"
<< " internalEdges:" << mesh().nInternalEdges()
<< " total number of edges:" << mesh().nEdges()
<< abort(FatalError);
return -1;
}
label eBoundaryMesh::findPatchID(const word& patchName) const
{
const ePatchList& patches = *this;
forAll (patches, patchI)
{
if (patches[patchI].name() == patchName)
{
return patchI;
}
}
// Patch not found
return -1;
}
// writeData member function required by regIOobject
bool eBoundaryMesh::writeData(Ostream& os) const
{
const ePatchList& patches = *this;
os << patches.size() << nl << token::BEGIN_LIST << incrIndent << nl;
forAll(patches, patchi)
{
os << indent << patches[patchi].name() << nl
<< indent << token::BEGIN_BLOCK << nl
<< incrIndent << patches[patchi] << decrIndent
<< indent << token::END_BLOCK << endl;
}
os << decrIndent << token::END_LIST;
// Check state of IOstream
os.check("eBoundaryMesh::writeData(Ostream& os) const");
return os.good();
}
Ostream& operator<<(Ostream& os, const eBoundaryMesh& bm)
{
bm.writeData(os);
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,144 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
eBoundaryMesh
SourceFiles
eBoundaryMesh.C
\*---------------------------------------------------------------------------*/
#ifndef eBoundaryMesh_H
#define eBoundaryMesh_H
#include "ePatchList.H"
#include "wordList.H"
#include "pointField.H"
#include "regIOobject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
class eMesh;
/*---------------------------------------------------------------------------*\
Class eBoundaryMesh Declaration
\*---------------------------------------------------------------------------*/
class eBoundaryMesh
:
public ePatchList,
public regIOobject
{
// private data
//- Reference to mesh
const eMesh& mesh_;
//- Disallow construct as copy
eBoundaryMesh(const eBoundaryMesh&);
//- Disallow assignment
void operator=(const eBoundaryMesh&);
public:
//- Runtime type information
TypeName("eBoundaryMesh");
// Constructors
//- Construct from dictionary
eBoundaryMesh
(
const IOobject& io,
const eMesh& mesh
);
//- Construct given size
eBoundaryMesh
(
const IOobject& io,
const eMesh& mesh,
const label size
);
// Destructor - default
// Member functions
// Access
//- Return the mesh reference
const eMesh& mesh() const;
//- Return a list of patch types
wordList types() const;
//- Return a list of patch names
wordList names() const;
//- Return a list of patch sizes
labelList patchSizes() const;
//- Return a list of patch starts
labelList patchStarts() const;
//- Return patch index for a given edge label
label whichPatch(const label edgeIndex) const;
//- Find patch index given a name
label findPatchID(const word& patchName) const;
//- writeData member function required by regIOobject
bool writeData(Ostream&) const;
// Input / Output
//- Read from input stream
void readFromInputStream();
// Ostream operator
friend Ostream& operator<<(Ostream&, const eBoundaryMesh&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,482 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
Description
Mesh needed to do edge-based addressing.
\*---------------------------------------------------------------------------*/
#include "eMesh.H"
#include "demandDrivenData.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(eMesh, 0);
word eMesh::meshSubDir = "eMesh";
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void eMesh::clearGeom() const
{
if (debug)
{
Info<< "void eMesh::clearGeom() const : "
<< "Clearing geometry" << endl;
}
}
void eMesh::clearAddressing() const
{
if (debug)
{
Info<< "void eMesh::clearAddressing() const : "
<< "Clearing addressing" << endl;
}
edges_.clear();
deleteDemandDrivenData(pePtr_);
deleteDemandDrivenData(epPtr_);
deleteDemandDrivenData(fePtr_);
deleteDemandDrivenData(efPtr_);
}
// Helper function to isolate points on triangular faces
label eMesh::findIsolatedPoint(const face& f, const edge& e) const
{
// Check the first point
if ( f[0] != e.start() && f[0] != e.end() )
{
return f[0];
}
// Check the second point
if ( f[1] != e.start() && f[1] != e.end() )
{
return f[1];
}
// Check the third point
if ( f[2] != e.start() && f[2] != e.end() )
{
return f[2];
}
return -1;
}
//- Helper function to determine the orientation of a triangular face
label eMesh::edgeDirection(const face& f, const edge& e) const
{
if
(
(f[0] == e.start() && f[1] == e.end())
|| (f[1] == e.start() && f[2] == e.end())
|| (f[2] == e.start() && f[0] == e.end())
)
{
// Return counter-clockwise
return 1;
}
else
{
// Return clockwise
return -1;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
eMesh::eMesh(const polyMesh& pMesh, const word& subDir)
:
objectRegistry(pMesh.time()),
mesh_(pMesh),
edges_
(
IOobject
(
"edges",
mesh_.facesInstance(),
subDir,
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
)
),
boundary_
(
IOobject
(
"edgeBoundary",
mesh_.facesInstance(),
subDir,
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
*this
),
pePtr_(NULL),
epPtr_(NULL),
fePtr_(NULL),
efPtr_(NULL)
{
if (debug)
{
Info << "eMesh::eMesh(const polyMesh&, const word&) : "
<< "Creating eMesh from polyMesh"
<< endl;
}
// Re-initialize / override meshSubDir
meshSubDir = subDir;
// Try to read from disk.
if (edges_.headerOk() && boundary_.headerOk())
{
// Set sizes
nEdges_ = edges_.size();
nInternalEdges_ = boundary_[0].start();
}
else
{
// Could not read ordered edges, so calculate it instead.
calcOrderedEdgeList();
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
eMesh::~eMesh()
{
clearOut();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
fileName eMesh::meshDir() const
{
return mesh_.dbDir();
}
fileName eMesh::meshSubDirectory() const
{
return mesh_.dbDir()/meshSubDir;
}
const Time& eMesh::time() const
{
return mesh_.time();
}
label eMesh::nEdges() const
{
return nEdges_;
}
label eMesh::nInternalEdges() const
{
return nInternalEdges_;
}
const edgeList& eMesh::edges() const
{
return edges_;
}
void eMesh::addEdgePatches(const List<ePatch*>& p)
{
if (debug)
{
Info << "void eMesh::addEdgePatches(const List<ePatch*>& p) : "
<< "Adding patches to eMesh" << endl;
}
if (boundary().size() > 0)
{
FatalErrorIn
(
"void eMesh::addEdgePatches(const List<ePatch*>& p)"
)
<< "Boundary already exists"
<< abort(FatalError);
}
boundary_.setSize(p.size());
forAll(p, patchI)
{
boundary_.set(patchI, p[patchI]);
}
}
const objectRegistry& eMesh::db() const
{
return mesh_.db();
}
const eBoundaryMesh& eMesh::boundary() const
{
return boundary_;
}
void eMesh::resetPrimitives
(
edgeList& edges,
labelListList& faceEdges,
labelListList& edgeFaces,
const labelList& patchSizes,
const labelList& patchStarts,
const bool reUse,
const bool storePrimitives
)
{
// Clear out geometry and addressing
clearOut();
// Initialize pointers for storage
if (storePrimitives)
{
fePtr_ =
(
new labelListIOList
(
IOobject
(
"faceEdges",
mesh_.facesInstance(),
meshSubDir,
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
)
)
);
efPtr_ =
(
new labelListIOList
(
IOobject
(
"edgeFaces",
mesh_.facesInstance(),
meshSubDir,
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
)
)
);
if (reUse)
{
edges_.transfer(edges);
fePtr_->transfer(faceEdges);
efPtr_->transfer(edgeFaces);
}
else
{
edges_ = edges;
fePtr_->operator=(faceEdges);
efPtr_->operator=(edgeFaces);
}
}
// Reset patch sizes and starts
forAll(boundary_, patchI)
{
boundary_[patchI] = ePatch
(
boundary_[patchI].name(),
patchSizes[patchI],
patchStarts[patchI],
patchI,
boundary_
);
}
// Reset size information
nEdges_ = edges.size();
nInternalEdges_ = boundary_[0].start();
// Set mesh files as changed
setInstance(time().timeName());
}
//- Clear demand-driven data
void eMesh::clearOut() const
{
clearGeom();
clearAddressing();
}
//- Set the instance for mesh files
void eMesh::setInstance(const fileName& inst)
{
if (debug)
{
Info<< "void eMesh::setInstance(const fileName& inst) : "
<< "Resetting file instance to " << inst << endl;
}
if (edges_.size())
{
edges_.writeOpt() = IOobject::AUTO_WRITE;
edges_.instance() = inst;
}
if (efPtr_)
{
efPtr_->writeOpt() = IOobject::AUTO_WRITE;
efPtr_->instance() = inst;
}
if (fePtr_)
{
fePtr_->writeOpt() = IOobject::AUTO_WRITE;
fePtr_->instance() = inst;
}
// Write out boundary information only
// if all others are being written out
if (edges_.size() && efPtr_ && fePtr_)
{
boundary_.writeOpt() = IOobject::AUTO_WRITE;
boundary_.instance() = inst;
}
}
const labelListList& eMesh::pointEdges() const
{
if (!pePtr_)
{
calcPointEdges();
}
return *pePtr_;
}
const labelListList& eMesh::edgePoints() const
{
if (!epPtr_)
{
calcEdgePoints();
}
return *epPtr_;
}
const labelListList& eMesh::faceEdges() const
{
if (!fePtr_)
{
calcFaceEdges();
}
return *fePtr_;
}
const labelListList& eMesh::edgeFaces() const
{
if (!efPtr_)
{
calcEdgeFaces();
}
return *efPtr_;
}
bool eMesh::write() const
{
if (edges_.size())
{
edges_.write();
}
if (efPtr_)
{
efPtr_->write();
}
if (fePtr_)
{
fePtr_->write();
}
if (edges_.size() && efPtr_ && fePtr_)
{
boundary_.write();
}
return false;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
bool eMesh::operator!=(const eMesh& m) const
{
return &m != this;
}
bool eMesh::operator==(const eMesh& m) const
{
return &m == this;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,251 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
eMesh
Description
Mesh class to store edge-based connectivity structures.
SourceFiles
eMesh.C
\*---------------------------------------------------------------------------*/
#ifndef eMesh_H
#define eMesh_H
#include "Time.H"
#include "polyMesh.H"
#include "edgeIOList.H"
#include "labelIOList.H"
#include "eBoundaryMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
/*---------------------------------------------------------------------------*\
Class eMesh Declaration
\*---------------------------------------------------------------------------*/
class eMesh
:
public objectRegistry
{
// Private data
const polyMesh& mesh_;
//- Ordered edge-list
mutable edgeIOList edges_;
//- Boundary mesh
mutable eBoundaryMesh boundary_;
// Primitive size data
//- Number of edges
mutable label nEdges_;
//- Number of internal edges
mutable label nInternalEdges_;
// Demand-driven data
//- Mapping between ordered edge-data and regular edges.
labelList reverseEdgeMap_;
//- Point-edges
mutable labelListList* pePtr_;
//- Edge-points
mutable labelListList* epPtr_;
//- Face-edges
mutable labelListIOList* fePtr_;
//- Edge-faces
mutable labelListIOList* efPtr_;
// Private Member Functions
//- Helper function to isolate points on triangular faces
label findIsolatedPoint(const face& f, const edge& e) const;
//- Helper function to determine the orientation of a triangular face
label edgeDirection(const face& f, const edge& e) const;
//- Disallow default bitwise copy construct
eMesh(const eMesh&);
//- Disallow default bitwise assignment
void operator=(const eMesh&);
// Private member functions to calculate demand driven data
//- Calculate ordered edges (and edgeFaces)
void calcOrderedEdgeList();
//- Calculate ordered edgePoints
void calcEdgePoints() const;
//- Calculate point-edges
void calcPointEdges() const;
//- Calculate face-edges
void calcFaceEdges() const;
//- Calculate edge-faces
void calcEdgeFaces() const;
//- Clear geometry
void clearGeom() const;
//- Clear addressing
void clearAddressing() const;
public:
// Public typedefs
typedef eMesh Mesh;
typedef eBoundaryMesh BoundaryMesh;
//- Runtime type information
TypeName("eMesh");
//- Return the mesh sub-directory name (usually "eMesh")
static word meshSubDir;
// Constructors
//- Construct from IOobject and polyMesh reference
eMesh(const polyMesh& m, const word& subDir = eMesh::meshSubDir);
// Destructor
virtual ~eMesh();
// Member Functions
// Helpers
//- Add boundary patches. Constructor helper
void addEdgePatches(const List<ePatch*> &);
// Database
//- Return reference to the mesh database
virtual const objectRegistry& db() const;
//- Return the base mesh directory (dbDir())
fileName meshDir() const;
//- Return the local mesh directory (dbDir()/meshSubDir)
fileName meshSubDirectory() const;
//- Return reference to time
const Time& time() const;
//- Mesh size parameters
label nEdges() const;
label nInternalEdges() const;
// Primitive mesh data
//- Return constant reference to the ordered edge-list
const edgeList& edges() const;
// Access
//- Return constant reference to boundary mesh
const eBoundaryMesh& boundary() const;
// Demand-driven data
//- Return constant reference to the pointEdges list
const labelListList& pointEdges() const;
//- Return constant reference to the edgePoints list
//- Valid for 3D tetrahedral meshes only.
const labelListList& edgePoints() const;
//- Return constant reference to the faceEdges list
const labelListList& faceEdges() const;
//- Return constant reference to the edgeFaces list
const labelListList& edgeFaces() const;
// Reset primitive data
void resetPrimitives
(
edgeList& edges,
labelListList& faceEdges,
labelListList& edgeFaces,
const labelList& patchSizes,
const labelList& patchStarts,
const bool reUse,
const bool storePrimitives
);
//- Clear demand-driven data
void clearOut() const;
//- Set the instance for mesh files
void setInstance(const fileName& inst);
//- Write mesh
virtual bool write() const;
// Member Operators
bool operator!=(const eMesh& m) const;
bool operator==(const eMesh& m) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,467 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
Description
\*---------------------------------------------------------------------------*/
#include "eMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Calculate ordered edges
void eMesh::calcOrderedEdgeList()
{
if (debug)
{
Info<< "void eMesh::calcOrderedEdges() const : "
<< "Calculating ordered edges" << endl;
}
if (edges_.size() || boundary_.size())
{
FatalErrorIn
(
"void eMesh::calcOrderedEdges() const"
) << "Ordered edges already allocated."
<< abort(FatalError);
}
// Set size first.
nEdges_ = mesh_.nEdges();
nInternalEdges_ = 0;
boundary_.setSize(mesh_.boundaryMesh().size());
// Allocate lists for re-ordering
labelList edgePatch(nEdges_, -1);
labelList edgePatchStarts(mesh_.boundaryMesh().size(), -1);
labelList edgePatchSizes(mesh_.boundaryMesh().size(), 0);
reverseEdgeMap_.setSize(nEdges_);
// Obtain connectivity from primitive mesh
const edgeList& edges = mesh_.edges();
const labelListList& fEdges = mesh_.faceEdges();
// Edge-patches are the same as faces
for (label i = mesh_.nInternalFaces(); i < mesh_.nFaces(); i++)
{
const labelList& fEdge = fEdges[i];
forAll(fEdge, edgeI)
{
edgePatch[fEdge[edgeI]] = mesh_.boundaryMesh().whichPatch(i);
}
}
// Loop through edgePatch and renumber internal edges
forAll(edgePatch, edgeI)
{
if (edgePatch[edgeI] == -1)
{
reverseEdgeMap_[edgeI] = nInternalEdges_++;
}
else
{
edgePatchSizes[edgePatch[edgeI]]++;
}
}
// Calculate patch-starts
label startCount = nInternalEdges_;
forAll(edgePatchStarts, patchI)
{
edgePatchStarts[patchI] = startCount;
startCount += edgePatchSizes[patchI];
}
// Now renumber boundary edges
labelList patchCount(edgePatchStarts);
forAll(edgePatch, edgeI)
{
if (edgePatch[edgeI] > -1)
{
reverseEdgeMap_[edgeI] = patchCount[edgePatch[edgeI]]++;
}
}
// Renumber and fill in edges
edges_.setSize(nEdges_, edge(-1,-1));
forAll(edges, edgeI)
{
edges_[reverseEdgeMap_[edgeI]] = edges[edgeI];
}
// Now set the boundary, copy name. (type is default)
forAll(boundary_, patchI)
{
boundary_.set
(
patchI,
ePatch::New
(
ePatch::typeName_(),
mesh_.boundaryMesh()[patchI].name(),
edgePatchSizes[patchI],
edgePatchStarts[patchI],
patchI,
boundary_
)
);
}
// Force calculation of other data...
calcEdgeFaces();
calcFaceEdges();
}
void eMesh::calcEdgePoints() const
{
if (debug)
{
Info<< "void eMesh::calcEdgePoints() const : "
<< "Calculating ordered edgePoints" << endl;
}
if (epPtr_)
{
FatalErrorIn
(
"void eMesh::calcEdgePoints() const"
) << "edgePoints already allocated."
<< abort(FatalError);
}
if (!efPtr_->size())
{
FatalErrorIn
(
"void eMesh::calcEdgePoints() const"
) << "edgeFaces doesn't exist."
<< abort(FatalError);
}
// Size the list
epPtr_ = new labelListList(nEdges_);
labelListList& edgePoints = *epPtr_;
// EdgePoints are ordered such that points appear in
// counter-clockwise direction around point[0] of the edge.
// If a boundary point is present, that must begin the list.
// NOTE: Will work only on tetrahedral meshes!
bool found;
label faceIndex = -1, cellIndex = -1;
const labelList& owner = mesh_.faceOwner();
const labelList& neighbour = mesh_.faceNeighbour();
const cellList& cells = mesh_.cells();
const faceList& faces = mesh_.faces();
const labelListList& eFaces = this->edgeFaces();
for (label eIndex=0; eIndex < nEdges_; eIndex++)
{
const edge& e = edges_[eIndex];
const labelList& eFace = eFaces[eIndex];
// Size the list
edgePoints[eIndex].setSize(eFace.size(), -1);
if (eIndex < nInternalEdges_)
{
// Pick the first face and start with that
faceIndex = eFace[0];
}
else
{
// Need to find a properly oriented start-face
forAll(eFace, faceI)
{
if
(
(!mesh_.isInternalFace(eFace[faceI]))
&& (edgeDirection(faces[eFace[faceI]], e) == 1)
)
{
faceIndex = eFace[faceI];
break;
}
}
}
// Shuffle vertices to appear in CCW order
for (label j=0; j < eFace.size(); j++)
{
const face& f = faces[faceIndex];
// Add the isolated point
edgePoints[eIndex][j] = findIsolatedPoint(f, e);
// Figure out how this edge is oriented.
if (edgeDirection(f, e) == 1)
{
// Counter-clockwise. Pick the owner.
cellIndex = owner[faceIndex];
}
else
if (mesh_.isInternalFace(faceIndex))
{
// Clockwise. Pick the neighbour.
cellIndex = neighbour[faceIndex];
}
else
{
// Looks like we've hit a boundary face. Break out.
break;
}
const cell& cellToCheck = cells[cellIndex];
found = false;
// Assuming tet-cells,
// Loop through edgeFaces and get the next face
forAll(eFace, faceI)
{
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[0]
)
{
faceIndex = cellToCheck[0];
found = true; break;
}
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[1]
)
{
faceIndex = cellToCheck[1];
found = true; break;
}
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[2]
)
{
faceIndex = cellToCheck[2];
found = true; break;
}
if
(
eFace[faceI] != faceIndex
&& eFace[faceI] == cellToCheck[3]
)
{
faceIndex = cellToCheck[3];
found = true; break;
}
}
# ifdef FULLDEBUG
if (!found)
{
// Something's terribly wrong
FatalErrorIn
(
"void eMesh::calcEdgePoints() const"
)
<< " Failed to determine a vertex ring. " << nl
<< " edgeFaces connectivity is inconsistent. " << nl
<< "Edge: " << eIndex << ":: " << e << nl
<< "edgeFaces: " << eFace
<< abort(FatalError);
}
# endif
}
}
}
void eMesh::calcPointEdges() const
{
if (debug)
{
Info<< "void eMesh::calcPointEdges() const : "
<< "Calculating PointEdges" << endl;
}
if (pePtr_)
{
FatalErrorIn
(
"void eMesh::calcPointEdges() const"
) << "pePtr_ already allocated"
<< abort(FatalError);
}
pePtr_ = new labelListList(mesh_.nPoints());
invertManyToMany(mesh_.nPoints(), edges(), *pePtr_);
}
void eMesh::calcFaceEdges() const
{
if (debug)
{
Info<< "void eMesh::calcFaceEdges() const : "
<< "Calculating FaceEdges" << endl;
}
if (fePtr_)
{
FatalErrorIn
(
"void eMesh::calcFaceEdges() const"
) << "fePtr_ already allocated"
<< abort(FatalError);
}
fePtr_ =
(
new labelListIOList
(
IOobject
(
"faceEdges",
mesh_.facesInstance(),
meshSubDir,
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
mesh_.nFaces()
)
);
labelListIOList& faceEdges = *fePtr_;
// If the read was successful, return.
if (faceEdges.headerOk())
{
return;
}
// If edgeFaces exists, use that.
if (efPtr_)
{
invertManyToMany(mesh_.nFaces(), edgeFaces(), faceEdges);
}
else
{
FatalErrorIn
(
"void eMesh::calcFaceEdges() const"
) << "Cannot calculate faceEdges."
<< abort(FatalError);
}
}
void eMesh::calcEdgeFaces() const
{
if (debug)
{
Info<< "void eMesh::calcEdgeFaces() const : "
<< "Calculating EdgeFaces" << endl;
}
if (efPtr_)
{
FatalErrorIn
(
"void eMesh::calcEdgeFaces() const"
) << "efPtr_ already allocated."
<< abort(FatalError);
}
efPtr_ =
(
new labelListIOList
(
IOobject
(
"edgeFaces",
mesh_.facesInstance(),
meshSubDir,
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
nEdges_
)
);
labelListIOList& edgeFaces = *efPtr_;
// If the read was successful, return.
if (edgeFaces.headerOk())
{
return;
}
if (fePtr_)
{
// If faceEdges exists, use that.
invertManyToMany(nEdges_, faceEdges(), edgeFaces);
}
else
{
// Obtain connectivity from primitive mesh.
if (!reverseEdgeMap_.size())
{
FatalErrorIn
(
"void eMesh::calcEdgeFaces() const"
) << "reverseEdgeMap has not been allocated."
<< abort(FatalError);
}
const labelListList& eFaces = mesh_.edgeFaces();
forAll(eFaces, edgeI)
{
edgeFaces[reverseEdgeMap_[edgeI]] = eFaces[edgeI];
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,138 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
Description
\*---------------------------------------------------------------------------*/
#include "ePatch.H"
#include "addToRunTimeSelectionTable.H"
#include "eBoundaryMesh.H"
#include "eMesh.H"
#include "polyMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(ePatch, 0);
defineRunTimeSelectionTable(ePatch, word);
defineRunTimeSelectionTable(ePatch, dictionary);
addToRunTimeSelectionTable(ePatch, ePatch, word);
addToRunTimeSelectionTable(ePatch, ePatch, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::ePatch::ePatch
(
const word& name,
const label size,
const label start,
const label index,
const eBoundaryMesh& bm
)
:
patchIdentifier(name, index),
boundaryMesh_(bm),
start_(start),
size_(size)
{}
// Construct from dictionary
ePatch::ePatch
(
const word& name,
const dictionary& dict,
const label index,
const eBoundaryMesh& bm
)
:
patchIdentifier(name, dict, index),
boundaryMesh_(bm),
start_(readLabel(dict.lookup("start"))),
size_(readLabel(dict.lookup("size")))
{}
ePatch::ePatch(const ePatch& p, const eBoundaryMesh& bm)
:
patchIdentifier(p, p.index()),
boundaryMesh_(bm),
start_(p.start()),
size_(p.size())
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
ePatch::~ePatch()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const eBoundaryMesh& ePatch::boundaryMesh() const
{
return boundaryMesh_;
}
void ePatch::write(Ostream& os) const
{
os.writeKeyword("type") << type() << token::END_STATEMENT << nl;
os.writeKeyword("start") << start() << token::END_STATEMENT << nl;
os.writeKeyword("size") << size() << token::END_STATEMENT << nl;
patchIdentifier::write(os);
}
//- Assignment operator
void ePatch::operator=(const ePatch& p)
{
patchIdentifier::operator=(p);
start_ = p.start_;
size_ = p.size_;
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const ePatch& p)
{
p.write(os);
os.check("Ostream& operator<<(Ostream& f, const faPatch& p)");
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,265 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
ePatch
SourceFiles
ePatch.C
newEPatch.C
\*---------------------------------------------------------------------------*/
#ifndef ePatch_H
#define ePatch_H
#include "patchIdentifier.H"
#include "labelList.H"
#include "pointField.H"
#include "typeInfo.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class eBoundaryMesh;
/*---------------------------------------------------------------------------*\
Class ePatch Declaration
\*---------------------------------------------------------------------------*/
class ePatch
:
public patchIdentifier
{
private:
// Private data
//- Reference to boundary mesh
const eBoundaryMesh& boundaryMesh_;
//- Start index of the patch
label start_;
//- Size of the patch
label size_;
// Demand-driven private data
// Private Member Functions
//- Disallow construct as copy
ePatch(const ePatch&);
protected:
// The ePatch geometry initialisation is called by eBoundaryMesh
friend class eBoundaryMesh;
//- Initialise the calculation of the patch geometry
virtual void initGeometry()
{}
//- Calculate the patch geometry
virtual void calcGeometry()
{}
//- Initialise the patches for moving points
virtual void initMovePoints(const pointField&)
{}
//- Correct patch after moving points
virtual void movePoints(const pointField&)
{}
//- Initialise the update of the patch topology
virtual void initUpdateMesh()
{}
//- Update of the patch topology
virtual void updateMesh()
{}
public:
typedef eBoundaryMesh BoundaryMesh;
//- Runtime type information
TypeName("patch");
// Declare run-time constructor selection tables
declareRunTimeSelectionTable
(
autoPtr,
ePatch,
word,
(
const word& name,
const label size,
const label start,
const label index,
const eBoundaryMesh& bm
),
(name, size, start, index, bm)
);
declareRunTimeSelectionTable
(
autoPtr,
ePatch,
dictionary,
(
const word& name,
const dictionary& dict,
const label index,
const eBoundaryMesh& bm
),
(name, dict, index, bm)
);
// Constructors
//- Construct from components
ePatch
(
const word& name,
const label size,
const label start,
const label index,
const eBoundaryMesh& bm
);
//- Construct from dictionary
ePatch
(
const word& name,
const dictionary& dict,
const label index,
const eBoundaryMesh& bm
);
//- Construct as copy, resetting the boundary mesh
ePatch(const ePatch&, const eBoundaryMesh&);
// Selectors
//- Return a pointer to a new patch created on freestore from
// components
static autoPtr<ePatch> New
(
const word& patchType,
const word& name,
const label size,
const label start,
const label index,
const eBoundaryMesh& bm
);
//- Return a pointer to a new patch created
// on freestore from dictionary
static autoPtr<ePatch> New
(
const word& name,
const dictionary& dict,
const label index,
const eBoundaryMesh& bm
);
// Destructor
virtual ~ePatch();
// Member Functions
//- Return the index of this patch in the boundaryMesh
label index() const
{
return patchIdentifier::index();
}
//- Return boundaryMesh reference
const eBoundaryMesh& boundaryMesh() const;
//- Return true if this patch is coupled
virtual bool coupled() const
{
return false;
}
//- Patch start in edge list
label start() const
{
return start_;
}
//- Patch size
virtual label size() const
{
return size_;
}
//- Slice list to patch
template<class T>
typename List<T>::subList patchSlice(const List<T>& l) const
{
return typename List<T>::subList(l, size(), start());
}
//- Write
virtual void write(Ostream&) const;
//- Assignment operator
void operator=(const ePatch&);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const ePatch&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,56 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
Type
ePatchList
Description
Container classes for ePatch
\*---------------------------------------------------------------------------*/
#ifndef ePatchList_H
#define ePatchList_H
#include "ePatch.H"
#include "PtrList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
typedef PtrList<ePatch> ePatchList;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,118 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
Description
\*---------------------------------------------------------------------------*/
#include "ePatch.H"
#include "dictionary.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
autoPtr<ePatch> ePatch::New
(
const word& patchType,
const word& name,
const label size,
const label start,
const label index,
const eBoundaryMesh& bm
)
{
if (debug)
{
Info<< "ePatch::New(const word&, const word&, const label, "
"const label, const label, const eBoundaryMesh&) : "
"constructing ePatch"
<< endl;
}
wordConstructorTable::iterator cstrIter =
wordConstructorTablePtr_->find(patchType);
if (cstrIter == wordConstructorTablePtr_->end())
{
FatalErrorIn
(
"ePatch::New(const word&, const word&, const label, "
"const label, const label, const eBoundaryMesh&) "
) << "Unknown ePatch type " << patchType << " for patch " << name
<< endl << endl
<< "Valid ePatch types are :" << endl
<< wordConstructorTablePtr_->toc()
<< exit(FatalError);
}
return autoPtr<ePatch>(cstrIter()(name, size, start, index, bm));
}
autoPtr<ePatch> ePatch::New
(
const word& name,
const dictionary& dict,
const label index,
const eBoundaryMesh& bm
)
{
if (debug)
{
Info<< "ePatch::New(const word&, const dictionary&, const label, "
"const eBoundaryMesh&) : constructing ePatch"
<< endl;
}
word patchType(dict.lookup("type"));
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(patchType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalIOErrorIn
(
"ePatch::New(const word&, const dictionary&, "
"const label, const eBoundaryMesh&)",
dict
) << "Unknown ePatch type " << patchType << endl << endl
<< "Valid ePatch types are :" << endl
<< dictionaryConstructorTablePtr_->toc()
<< exit(FatalIOError);
}
return autoPtr<ePatch>(cstrIter()(name, dict, index, bm));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

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,173 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
fluxCorrector
Description
Implementation of the fluxCorrector base class
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
#include "fluxCorrector.H"
#include "dlLibraryTable.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(fluxCorrector, 0);
defineRunTimeSelectionTable(fluxCorrector, mesh);
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
autoPtr<fluxCorrector> fluxCorrector::New
(
const fvMesh& mesh,
const dictionary& dict
)
{
// Check if an optional entry was specified
if (dict.found("fluxCorrector"))
{
word correctorTypeName(dict.lookup("fluxCorrector"));
// Open any supplied libraries in dictionary
dlLibraryTable::open
(
dict,
"fluxCorrectorLibs",
meshConstructorTablePtr_
);
if (!meshConstructorTablePtr_)
{
FatalErrorIn
(
"autoPtr<fluxCorrector> fluxCorrector::New"
"(const fvMesh& mesh, const dictionary& dict)"
) << "fluxCorrector table is empty"
<< exit(FatalError);
}
correctorTypeName = word(dict.lookup("fluxCorrector"));
meshConstructorTable::iterator cstrIter =
meshConstructorTablePtr_->find(correctorTypeName);
if (cstrIter == meshConstructorTablePtr_->end())
{
FatalErrorIn
(
"autoPtr<fluxCorrector> fluxCorrector::New"
"(const fvMesh& mesh, const dictionary& dict)"
) << "Unknown fluxCorrector type " << correctorTypeName
<< endl << endl
<< "Valid fluxCorrector types are: " << endl
<< meshConstructorTablePtr_->toc()
<< exit(FatalError);
}
return autoPtr<fluxCorrector>(cstrIter()(mesh, dict));
}
// Return the default fluxCorrector
return autoPtr<fluxCorrector>(new fluxCorrector(mesh, dict));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Return reference to mesh
const fvMesh& fluxCorrector::mesh() const
{
return mesh_;
}
//- Return reference to dictionary
const dictionary& fluxCorrector::dict() const
{
return dict_;
}
//- Is flux-correction required?
bool fluxCorrector::required() const
{
// Support for cases which do not involve a flow-solver.
// Notify the user, just in case.
Info << " ~~~ No flux correction ~~~ " << endl;
return false;
}
//- Interpolate fluxes to a specified list of faces
void fluxCorrector::interpolateFluxes(const labelList& faces) const
{
if (required())
{
// Throw an error stating that the scheme hasn't been implemented
notImplemented
(
"void fluxCorrector::interpolateFluxes"
"(const labelList& faces) const"
);
}
}
//- Update fluxes in the registry, if required
void fluxCorrector::updateFluxes() const
{
if (required())
{
// Throw an error stating that the scheme hasn't been implemented
notImplemented("void fluxCorrector::updateFluxes() const");
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void fluxCorrector::operator=(const fluxCorrector& rhs)
{
// Check for assignment to self
if (this == &rhs)
{
FatalErrorIn("fluxCorrector::operator=(const fluxCorrector&)")
<< "Attempted assignment to self"
<< abort(FatalError);
}
}
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,153 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
fluxCorrector
Description
Virtual base class that deals with flux-correction after topo-changes.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
fluxCorrector.C
\*---------------------------------------------------------------------------*/
#ifndef fluxCorrector_H
#define fluxCorrector_H
#include "fvMesh.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward class declarations
/*---------------------------------------------------------------------------*\
Class fluxCorrector Declaration
\*---------------------------------------------------------------------------*/
class fluxCorrector
{
// Private data
//- Reference to fvMesh
const fvMesh& mesh_;
//- Reference to dictionary
const dictionary& dict_;
// Private Member Functions
//- Disallow default bitwise copy construct
fluxCorrector(const fluxCorrector&);
//- Disallow default bitwise assignment
void operator=(const fluxCorrector&);
public:
//- Runtime type information
TypeName("fluxCorrector");
// Declare run-time constructor selection tables
declareRunTimeSelectionTable
(
autoPtr,
fluxCorrector,
mesh,
(
const fvMesh& mesh,
const dictionary& dict
),
(mesh, dict)
);
// Selectors
//- Select constructed from fvMesh
static autoPtr<fluxCorrector> New
(
const fvMesh& mesh,
const dictionary& dict
);
// Constructors
//- Construct from fvMesh and dictionary
fluxCorrector(const fvMesh& mesh, const dictionary& dict)
:
mesh_(mesh),
dict_(dict)
{}
// Destructor
virtual ~fluxCorrector()
{}
// Member Functions
//- Return reference to mesh
const fvMesh& mesh() const;
//- Return reference to dictionary
const dictionary& dict() const;
//- Is flux-correction required?
virtual bool required() const;
//- Interpolate fluxes to a specified list of faces
virtual void interpolateFluxes(const labelList& faces) const;
//- Update fluxes in the registry
virtual void updateFluxes() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#ifdef NoRepository
#include "fluxCorrector.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,109 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
topoBoundaryMeshMapper
Description
This object provides mapping and fill-in information for boundary data
between the two meshes after the topological change. It is
constructed from mapPolyMesh.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#ifndef topoBoundaryMeshMapper_H
#define topoBoundaryMeshMapper_H
#include "PtrList.H"
#include "topoPatchMapper.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class topoBoundaryMeshMapper Declaration
\*---------------------------------------------------------------------------*/
class topoBoundaryMeshMapper
:
public PtrList<topoPatchMapper>
{
// Private Member Functions
//- Disallow default bitwise copy construct
topoBoundaryMeshMapper(const topoBoundaryMeshMapper&);
//- Disallow default bitwise assignment
void operator=(const topoBoundaryMeshMapper&);
public:
// Constructors
//- Construct from components
topoBoundaryMeshMapper
(
const fvMesh& mesh,
const mapPolyMesh& mpm,
const topoMapper& mapper
)
:
PtrList<topoPatchMapper>(mesh.boundary().size())
{
const fvBoundaryMesh& patches = mesh.boundary();
forAll (patches, patchI)
{
set
(
patchI,
new topoPatchMapper
(
patches[patchI],
mpm,
mapper
)
);
}
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,570 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
topoCellMapper
Description
Implementation of the topoCellMapper class
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
#include "IOmanip.H"
#include "topoMapper.H"
#include "mapPolyMesh.H"
#include "topoCellMapper.H"
#include "demandDrivenData.H"
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
//- Clear out local storage
void topoCellMapper::clearOut()
{
deleteDemandDrivenData(directAddrPtr_);
deleteDemandDrivenData(interpolationAddrPtr_);
deleteDemandDrivenData(weightsPtr_);
deleteDemandDrivenData(insertedCellLabelsPtr_);
deleteDemandDrivenData(volumesPtr_);
deleteDemandDrivenData(centresPtr_);
}
//- Calculate addressing for interpolative mapping
void topoCellMapper::calcAddressing() const
{
if (directAddrPtr_ || interpolationAddrPtr_)
{
FatalErrorIn("void topoCellMapper::calcAddressing() const")
<< "Addressing already calculated."
<< abort(FatalError);
}
// Allocate for inserted cell labels
label nInsertedCells = 0;
insertedCellLabelsPtr_ = new labelList(mesh_.nCells(), -1);
labelList& insertedCells = *insertedCellLabelsPtr_;
if (direct())
{
// Direct addressing, no weights
directAddrPtr_ = new labelList(mpm_.cellMap());
}
else
{
// Interpolative addressing
interpolationAddrPtr_ = new labelListList(mesh_.nCells());
labelListList& addr = *interpolationAddrPtr_;
const List<objectMap>& cfc = mpm_.cellsFromCellsMap();
forAll (cfc, cfcI)
{
// Get addressing
const labelList& mo = cfc[cfcI].masterObjects();
label cellI = cfc[cfcI].index();
if (addr[cellI].size() > 0)
{
FatalErrorIn("void topoCellMapper::calcAddressing() const")
<< "Master cell " << cellI
<< " mapped from cells " << mo
<< " is already destination for mapping."
<< abort(FatalError);
}
// Set master objects
addr[cellI] = mo;
}
// Do mapped cells. Note that this can already be set by cellsFromCells
// so check if addressing size still zero.
const labelList& cm = mpm_.cellMap();
forAll (cm, cellI)
{
// Mapped from a single cell
if (cm[cellI] > -1 && addr[cellI].empty())
{
addr[cellI] = labelList(1, cm[cellI]);
}
// Check for inserted cells without any addressing
if (cm[cellI] < 0 && addr[cellI].empty())
{
insertedCells[nInsertedCells++] = cellI;
}
}
}
// Shorten inserted cells to actual size
insertedCells.setSize(nInsertedCells);
if (nInsertedCells)
{
FatalErrorIn("void topoCellMapper::calcAddressing() const")
<< " Found " << nInsertedCells << " which are"
<< " not mapped from any parent cells." << nl
<< " List: " << nl
<< insertedCells
<< abort(FatalError);
}
}
//- Calculate inverse-distance weights for interpolative mapping
void topoCellMapper::calcInverseDistanceWeights() const
{
if (weightsPtr_)
{
FatalErrorIn
(
"void topoCellMapper::calcInverseDistanceWeights() const"
)
<< "Weights already calculated."
<< abort(FatalError);
}
// Fetch interpolative addressing
const labelListList& addr = addressing();
// Allocate memory
weightsPtr_ = new scalarListList(size());
scalarListList& w = *weightsPtr_;
// Obtain cell-centre information from old/new meshes
const vectorField& oldCentres = tMapper_.internalCentres();
const vectorField& newCentres = mesh_.cellCentres();
forAll(addr, cellI)
{
const labelList& mo = addr[cellI];
// Do mapped cells
if (mo.size() == 1)
{
w[cellI] = scalarList(1, 1.0);
}
else
{
// Map from masters, inverse-distance weights
scalar totalWeight = 0.0;
w[cellI] = scalarList(mo.size(), 0.0);
forAll (mo, oldCellI)
{
w[cellI][oldCellI] =
(
1.0/stabilise
(
magSqr
(
newCentres[cellI]
- oldCentres[mo[oldCellI]]
),
VSMALL
)
);
totalWeight += w[cellI][oldCellI];
}
// Normalize weights
scalar normFactor = (1.0/totalWeight);
forAll (mo, oldCellI)
{
w[cellI][oldCellI] *= normFactor;
}
}
}
}
//- Calculate intersection weights for conservative mapping
void topoCellMapper::calcIntersectionWeightsAndCentres() const
{
if (volumesPtr_ || centresPtr_)
{
FatalErrorIn
(
"void topoCellMapper::"
"calcIntersectionWeightsAndCentres() const"
)
<< "Weights already calculated."
<< abort(FatalError);
}
// Fetch interpolative addressing
const labelListList& addr = addressing();
// Allocate memory
volumesPtr_ = new List<scalarField>(size(), scalarField(0));
List<scalarField>& v = *volumesPtr_;
centresPtr_ = new List<vectorField>(size(), vectorField(0));
List<vectorField>& x = *centresPtr_;
// Obtain stored cell-centres
const vectorField& cellCentres = tMapper_.internalCentres();
// Fetch maps
const List<objectMap>& cfc = mpm_.cellsFromCellsMap();
const List<vectorField>& mapCellCentres = tMapper_.cellCentres();
const List<scalarField>& mapCellWeights = tMapper_.cellWeights();
// Fill in maps first
forAll(cfc, indexI)
{
x[cfc[indexI].index()] = mapCellCentres[indexI];
v[cfc[indexI].index()] = mapCellWeights[indexI];
}
// Now do mapped cells
forAll(addr, cellI)
{
const labelList& mo = addr[cellI];
// Check if this is indeed a mapped cell
if (mo.size() == 1 && x[cellI].empty() && v[cellI].empty())
{
x[cellI] = vectorField(1, cellCentres[mo[0]]);
v[cellI] = scalarField(1, 1.0);
}
}
}
const List<scalarField>& topoCellMapper::intersectionWeights() const
{
if (direct())
{
FatalErrorIn
(
"const List<scalarField>& "
"topoCellMapper::intersectionWeights() const"
) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError);
}
if (!volumesPtr_)
{
calcIntersectionWeightsAndCentres();
}
return *volumesPtr_;
}
const List<vectorField>& topoCellMapper::intersectionCentres() const
{
if (direct())
{
FatalErrorIn
(
"const List<vectorField>& "
"topoCellMapper::intersectionCentres() const"
) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError);
}
if (!centresPtr_)
{
calcIntersectionWeightsAndCentres();
}
return *centresPtr_;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
topoCellMapper::topoCellMapper
(
const mapPolyMesh& mpm,
const topoMapper& mapper
)
:
mesh_(mpm.mesh()),
mpm_(mpm),
tMapper_(mapper),
direct_(false),
directAddrPtr_(NULL),
interpolationAddrPtr_(NULL),
weightsPtr_(NULL),
insertedCellLabelsPtr_(NULL),
volumesPtr_(NULL),
centresPtr_(NULL)
{
// Check for possibility of direct mapping
if
(
(min(mpm_.cellMap()) > -1)
&& mpm_.cellsFromPointsMap().empty()
&& mpm_.cellsFromEdgesMap().empty()
&& mpm_.cellsFromFacesMap().empty()
&& mpm_.cellsFromCellsMap().empty()
)
{
direct_ = true;
}
else
{
direct_ = false;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
topoCellMapper::~topoCellMapper()
{
clearOut();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Return size
label topoCellMapper::size() const
{
return mesh_.nCells();
}
//- Return size before mapping
label topoCellMapper::sizeBeforeMapping() const
{
return mpm_.nOldCells();
}
//- Is the mapping direct
bool topoCellMapper::direct() const
{
return direct_;
}
//- Return direct addressing
const unallocLabelList& topoCellMapper::directAddressing() const
{
if (!direct())
{
FatalErrorIn
(
"const unallocLabelList& "
"topoCellMapper::directAddressing() const"
) << "Requested direct addressing for an interpolative mapper."
<< abort(FatalError);
}
if (!directAddrPtr_)
{
calcAddressing();
}
return *directAddrPtr_;
}
//- Return interpolation addressing
const labelListList& topoCellMapper::addressing() const
{
if (direct())
{
FatalErrorIn
(
"const labelListList& "
"topoCellMapper::addressing() const"
) << "Requested interpolative addressing for a direct mapper."
<< abort(FatalError);
}
if (!interpolationAddrPtr_)
{
calcAddressing();
}
return *interpolationAddrPtr_;
}
//- Return weights
const scalarListList& topoCellMapper::weights() const
{
if (direct())
{
FatalErrorIn
(
"const scalarListList& "
"topoCellMapper::weights() const"
) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError);
}
if (!weightsPtr_)
{
calcInverseDistanceWeights();
}
return *weightsPtr_;
}
//- Are there any inserted cells
bool topoCellMapper::insertedObjects() const
{
return insertedObjectLabels().size();
}
//- Return list of inserted cells
const labelList& topoCellMapper::insertedObjectLabels() const
{
if (!insertedCellLabelsPtr_)
{
calcAddressing();
}
return *insertedCellLabelsPtr_;
}
//- Conservatively map the internal field
template <class Type, class gradType>
void topoCellMapper::mapInternalField
(
const word& fieldName,
const Field<gradType>& gF,
Field<Type>& iF
) const
{
if (iF.size() != sizeBeforeMapping() || gF.size() != sizeBeforeMapping())
{
FatalErrorIn
(
"\n\n"
"void topoCellMapper::mapInternalField<Type>\n"
"(\n"
" const word& fieldName,\n"
" const Field<gradType>& gF,\n"
" Field<Type>& iF\n"
") const\n"
) << "Incompatible size before mapping." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << iF.size() << nl
<< " Gradient Field size: " << gF.size() << nl
<< " map size: " << sizeBeforeMapping() << nl
<< abort(FatalError);
}
// Fetch addressing
const labelListList& cAddressing = addressing();
const List<scalarField>& wC = intersectionWeights();
const List<vectorField>& xC = intersectionCentres();
// Fetch geometry
const vectorField& centres = tMapper_.internalCentres();
// Compute the integral of the source field
Type intSource = sum(iF * tMapper_.cellVolumes());
// Copy the original field
Field<Type> fieldCpy(iF);
// Resize to current dimensions
iF.setSize(size());
// Map the internal field
forAll(iF, cellI)
{
const labelList& addr = cAddressing[cellI];
iF[cellI] = pTraits<Type>::zero;
forAll(addr, cellJ)
{
// Accumulate volume-weighted Taylor-series interpolate
iF[cellI] +=
(
wC[cellI][cellJ] *
(
fieldCpy[addr[cellJ]]
+ (
gF[addr[cellJ]] &
(
xC[cellI][cellJ]
- centres[addr[cellJ]]
)
)
)
);
}
}
// Compute the integral of the target field
Type intTarget = sum(iF * mesh_.cellVolumes());
if (polyMesh::debug)
{
int oldP = Info().precision();
// Compare the global integral
Info << " Field : " << fieldName
<< " integral errors : "
<< setprecision(15)
<< " source : " << mag(intSource)
<< " target : " << mag(intTarget)
<< " norm : "
<< (mag(intTarget - intSource) / (mag(intSource) + VSMALL))
<< setprecision(oldP)
<< endl;
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void topoCellMapper::operator=(const topoCellMapper& rhs)
{
// Check for assignment to self
if (this == &rhs)
{
FatalErrorIn("void topoCellMapper::operator=")
<< "Attempted assignment to self"
<< abort(FatalError);
}
}
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,188 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
topoCellMapper
Description
This object provides mapping and fill-in information for internal
cell data between the two meshes after topological changes.
It is constructed from mapPolyMesh.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
topoCellMapper.C
\*---------------------------------------------------------------------------*/
#ifndef topoCellMapper_H
#define topoCellMapper_H
#include "morphFieldMapper.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
class topoMapper;
class mapPolyMesh;
/*---------------------------------------------------------------------------*\
Class topoCellMapper Declaration
\*---------------------------------------------------------------------------*/
class topoCellMapper
:
public morphFieldMapper
{
// Private data
//- Reference to polyMesh
const polyMesh& mesh_;
//- Reference to mapPolyMesh
const mapPolyMesh& mpm_;
//- Reference to the topoMapper
const topoMapper& tMapper_;
//- Is the mapping direct
bool direct_;
// Demand-driven private data
//- Direct addressing
mutable labelList* directAddrPtr_;
//- Interpolated addressing
mutable labelListList* interpolationAddrPtr_;
//- Inverse-distance weights
mutable scalarListList* weightsPtr_;
//- Inserted cells
mutable labelList* insertedCellLabelsPtr_;
//- Interpolation volumes
mutable List<scalarField>* volumesPtr_;
//- Interpolation centres
mutable List<vectorField>* centresPtr_;
// Private Member Functions
//- Disallow default bitwise copy construct
topoCellMapper(const topoCellMapper&);
//- Disallow default bitwise assignment
void operator=(const topoCellMapper&);
//- Calculate addressing for interpolative mapping
void calcAddressing() const;
//- Calculate inverse-distance weights for interpolative mapping
void calcInverseDistanceWeights() const;
//- Calculate intersection weights for conservative mapping
void calcIntersectionWeightsAndCentres() const;
//- Return intersection volume weights
const List<scalarField>& intersectionWeights() const;
//- Return intersection volume centres
const List<vectorField>& intersectionCentres() const;
//- Clear out local storage
void clearOut();
public:
// Constructors
//- Construct from mapPolyMesh
topoCellMapper
(
const mapPolyMesh& mpm,
const topoMapper& mapper
);
// Destructor
virtual ~topoCellMapper();
// Member Functions
//- Return size
virtual label size() const;
//- Return size before mapping
virtual label sizeBeforeMapping() const;
//- Is the mapping direct
virtual bool direct() const;
//- Return direct addressing
virtual const unallocLabelList& directAddressing() const;
//- Return interpolation addressing
virtual const labelListList& addressing() const;
//- Return weights
virtual const scalarListList& weights() const;
//- Are there any inserted cells
virtual bool insertedObjects() const;
//- Return list of inserted cells
virtual const labelList& insertedObjectLabels() const;
//- Conservatively map the internal field
template <class Type, class gradType>
void mapInternalField
(
const word& fieldName,
const Field<gradType>& gF,
Field<Type>& iF
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#ifdef NoRepository
# include "topoCellMapper.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,637 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
topoMapper
Description
Implementation of topoMapper
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
#include "fvc.H"
#include "topoMapper.H"
#include "fluxCorrector.H"
#include "topoCellMapper.H"
#include "leastSquaresGrad.H"
#include "topoSurfaceMapper.H"
#include "topoBoundaryMeshMapper.H"
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Store gradients of fields on the mesh prior to topology changes
template <class Type, class gradType>
void topoMapper::storeGradients
(
HashTable<autoPtr<gradType> >& gradTable
) const
{
// Define a few typedefs for convenience
typedef typename outerProduct<vector, Type>::type gCmptType;
typedef GeometricField<Type, fvPatchField, volMesh> volType;
typedef GeometricField<gCmptType, fvPatchField, volMesh> gVolType;
// Fetch all fields from registry
HashTable<const volType*> fields
(
mesh_.objectRegistry::lookupClass<volType>()
);
forAllConstIter(typename HashTable<const volType*>, fields, fIter)
{
const volType& field = *fIter();
// Compute the gradient.
tmp<gVolType> tGrad;
// If the fvSolution dictionary contains an entry,
// use that, otherwise, default to leastSquares
word gradName("grad(" + field.name() + ')');
if (mesh_.schemesDict().subDict("gradSchemes").found(gradName))
{
tGrad = fvc::grad(field);
}
else
{
tGrad = fv::leastSquaresGrad<Type>(mesh_).grad(field);
}
// Make a new entry, but don't register the field.
gradTable.insert
(
field.name(),
autoPtr<gradType>
(
new gradType
(
IOobject
(
tGrad().name(),
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
tGrad()
)
)
);
}
}
//- Fetch the gradient field (template specialisation)
template <>
const volVectorField& topoMapper::gradient(const word& name) const
{
if (!sGrads_.found(name))
{
FatalErrorIn
(
"const volVectorField& "
"topoMapper::gradient(const word& name) const"
) << nl << " Gradient for: " << name
<< " has not been stored."
<< abort(FatalError);
}
return sGrads_[name]();
}
//- Fetch the gradient field (template specialisation)
template <>
const volTensorField& topoMapper::gradient(const word& name) const
{
if (!vGrads_.found(name))
{
FatalErrorIn
(
"const volTensorField& "
"topoMapper::gradient(const word& name) const"
) << nl << " Gradient for: " << name
<< " has not been stored."
<< abort(FatalError);
}
return vGrads_[name]();
}
//- Store gradients prior to mesh reset
void topoMapper::storeGradients() const
{
storeGradients<scalar>(sGrads_);
storeGradients<vector>(vGrads_);
if (fvMesh::debug)
{
Info << "Registered volScalarFields: " << endl;
Info << sGrads_.toc() << endl;
Info << "Registered volVectorFields: " << endl;
Info << vGrads_.toc() << endl;
}
}
//- Store geometric information
void topoMapper::storeGeometry() const
{
if (cellCentresPtr_)
{
deleteDemandDrivenData(cellVolumesPtr_);
deleteDemandDrivenData(cellCentresPtr_);
patchAreasPtr_.clear();
patchCentresPtr_.clear();
}
// Set the cell-volumes pointer.
cellVolumesPtr_ = new scalarField(mesh_.cellVolumes());
// Set the cell-centres pointer.
cellCentresPtr_ = new vectorField(mesh_.cellCentres());
// Set patch-areas
patchAreasPtr_.setSize(mesh_.boundaryMesh().size());
forAll(mesh_.boundaryMesh(), patchI)
{
patchAreasPtr_.set
(
patchI,
new scalarField
(
mag(mesh_.boundaryMesh()[patchI].faceAreas())
)
);
}
// Set patch-centres.
patchCentresPtr_.setSize(mesh_.boundaryMesh().size());
forAll(mesh_.boundaryMesh(), patchI)
{
patchCentresPtr_.set
(
patchI,
new vectorField
(
mesh_.boundaryMesh()[patchI].faceCentres()
)
);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
topoMapper::~topoMapper()
{
clear();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Return reference to the mesh
const fvMesh& topoMapper::mesh() const
{
return mesh_;
}
//- Return reference to objectRegistry storing fields.
const objectRegistry& topoMapper::thisDb() const
{
return mesh_;
}
//- Set mapping information
void topoMapper::setMapper(const mapPolyMesh& mpm) const
{
if
(
cellMap_.valid() ||
surfaceMap_.valid() ||
boundaryMap_.valid()
)
{
FatalErrorIn
(
"void topoMapper::setMapper() const"
) << nl << " Mapper has already been set. "
<< abort(FatalError);
}
// Set pointers
cellMap_.set(new topoCellMapper(mpm, *this));
surfaceMap_.set(new topoSurfaceMapper(mpm, *this));
boundaryMap_.set(new topoBoundaryMeshMapper(mesh(), mpm, *this));
}
//- Set face weighting information
void topoMapper::setFaceWeights
(
const Xfer<List<scalarField> >& weights,
const Xfer<List<vectorField> >& centres
) const
{
faceWeights_.transfer(weights());
faceCentres_.transfer(centres());
}
//- Set cell weighting information
void topoMapper::setCellWeights
(
const Xfer<List<scalarField> >& weights,
const Xfer<List<vectorField> >& centres
) const
{
cellWeights_.transfer(weights());
cellCentres_.transfer(centres());
}
//- Fetch face weights
const List<scalarField>& topoMapper::faceWeights() const
{
return faceWeights_;
}
//- Fetch cell weights
const List<scalarField>& topoMapper::cellWeights() const
{
return cellWeights_;
}
//- Fetch face centres
const List<vectorField>& topoMapper::faceCentres() const
{
return faceCentres_;
}
//- Fetch cell centres
const List<vectorField>& topoMapper::cellCentres() const
{
return cellCentres_;
}
//- Store mesh information for the mapping stage
void topoMapper::storeMeshInformation() const
{
// Store field-gradients
storeGradients();
// Store geometry
storeGeometry();
}
//- Return stored cell volume information
const scalarField& topoMapper::cellVolumes() const
{
if (!cellVolumesPtr_)
{
FatalErrorIn
(
"const scalarField& topoMapper::cellVolumes()"
) << nl << " Pointer has not been set. "
<< abort(FatalError);
}
return *cellVolumesPtr_;
}
//- Return stored cell centre information
const vectorField& topoMapper::internalCentres() const
{
if (!cellCentresPtr_)
{
FatalErrorIn
(
"const vectorField& topoMapper::internalCentres()"
) << nl << " Pointer has not been set. "
<< abort(FatalError);
}
return *cellCentresPtr_;
}
//- Return stored patch areas information
const scalarField& topoMapper::patchAreas(const label i) const
{
if (!patchAreasPtr_.set(i))
{
FatalErrorIn
(
"const scalarField& topoMapper::patchAreas"
"(const label i) const"
) << nl << " Pointer has not been set at index: " << i
<< abort(FatalError);
}
return patchAreasPtr_[i];
}
//- Return stored patch centre information
const vectorField& topoMapper::patchCentres(const label i) const
{
if (!patchCentresPtr_.set(i))
{
FatalErrorIn
(
"const vectorField& topoMapper::patchCentres"
"(const label i) const"
) << nl << " Pointer has not been set at index: " << i
<< abort(FatalError);
}
return patchCentresPtr_[i];
}
// Conservatively map all volFields in the registry
template <class Type>
void topoMapper::conservativeMapVolFields() const
{
// Define a few typedefs for convenience
typedef typename outerProduct<vector, Type>::type gCmptType;
typedef GeometricField<Type, fvPatchField, volMesh> volType;
typedef GeometricField<gCmptType, fvPatchField, volMesh> gradVolType;
HashTable<const volType*> fields(mesh_.lookupClass<volType>());
// Store old-times before mapping
forAllIter(typename HashTable<const volType*>, fields, fIter)
{
volType& field = const_cast<volType&>(*fIter());
field.storeOldTimes();
}
// Fetch internal/boundary mappers
const topoCellMapper& fMap = volMap();
const topoBoundaryMeshMapper& bMap = boundaryMap();
// Now map all fields
forAllIter(typename HashTable<const volType*>, fields, fIter)
{
volType& field = const_cast<volType&>(*fIter());
if (fvMesh::debug)
{
Info << "Conservatively mapping "
<< field.typeName
<< ' ' << field.name()
<< endl;
}
// Map the internal field
fMap.mapInternalField
(
field.name(),
gradient<gradVolType>(field.name()).internalField(),
field.internalField()
);
// Map patch fields
forAll(bMap, patchI)
{
bMap[patchI].mapPatchField
(
field.name(),
field.boundaryField()[patchI]
);
}
// Set the field instance
field.instance() = field.mesh().thisDb().time().timeName();
}
}
// Conservatively map all surfaceFields in the registry
template <class Type>
void topoMapper::conservativeMapSurfaceFields() const
{
// Define a few typedefs for convenience
typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfType;
HashTable<const surfType*> fields(mesh_.lookupClass<surfType>());
// Store old-times before mapping
forAllIter(typename HashTable<const surfType*>, fields, fIter)
{
surfType& field = const_cast<surfType&>(*fIter());
field.storeOldTimes();
}
// Fetch internal/boundary mappers
const topoSurfaceMapper& fMap = surfaceMap();
const topoBoundaryMeshMapper& bMap = boundaryMap();
// Now map all fields
forAllIter(typename HashTable<const surfType*>, fields, fIter)
{
surfType& field = const_cast<surfType&>(*fIter());
if (fvMesh::debug)
{
Info << "Conservatively mapping "
<< field.typeName
<< ' ' << field.name()
<< endl;
}
// Map the internal field
fMap.mapInternalField
(
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();
}
}
//- Correct fluxes after topology changes, if required
void topoMapper::correctFluxes() const
{
if (surfaceFluxCorrector().required())
{
// Supply a list of inserted faces for interpolation
surfaceFluxCorrector().interpolateFluxes
(
surfaceMap().insertedObjectLabels()
);
// Update fluxes
surfaceFluxCorrector().updateFluxes();
}
}
//- Return volume mapper
const topoCellMapper& topoMapper::volMap() const
{
if (!cellMap_.valid())
{
FatalErrorIn
(
"const topoCellMapper& topoMapper::volMap()"
) << nl << " Volume mapper has not been set. "
<< abort(FatalError);
}
return cellMap_();
}
//- Return surface mapper
const topoSurfaceMapper& topoMapper::surfaceMap() const
{
if (!surfaceMap_.valid())
{
FatalErrorIn
(
"const topoSurfaceMapper& topoMapper::surfaceMap()"
) << nl << " Surface mapper has not been set. "
<< abort(FatalError);
}
return surfaceMap_();
}
//- Return boundary mapper
const topoBoundaryMeshMapper& topoMapper::boundaryMap() const
{
if (!boundaryMap_.valid())
{
FatalErrorIn
(
"const topoBoundaryMeshMapper& topoMapper::boundaryMap()"
) << nl << " Boundary mapper has not been set. "
<< abort(FatalError);
}
return boundaryMap_();
}
//- Return flux-corrector
const fluxCorrector& topoMapper::surfaceFluxCorrector() const
{
if (!fluxCorrector_.valid())
{
FatalErrorIn
(
"const fluxCorrector& topoMapper::surfaceFluxCorrector()"
) << nl << " fluxCorrector has not been set. "
<< abort(FatalError);
}
return fluxCorrector_();
}
//- Clear out member data
void topoMapper::clear()
{
// Clear out mappers
cellMap_.clear();
surfaceMap_.clear();
boundaryMap_.clear();
// Clear stored gradients
sGrads_.clear();
vGrads_.clear();
// Wipe out geomtry information
deleteDemandDrivenData(cellVolumesPtr_);
deleteDemandDrivenData(cellCentresPtr_);
patchAreasPtr_.clear();
patchCentresPtr_.clear();
// Clear maps
faceWeights_.clear();
cellWeights_.clear();
faceCentres_.clear();
cellCentres_.clear();
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void topoMapper::operator=(const topoMapper& rhs)
{
// Check for assignment to self
if (this == &rhs)
{
FatalErrorIn
(
"topoMapper::operator=(const topoMapper&)"
)
<< "Attempted assignment to self"
<< abort(FatalError);
}
}
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,245 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
topoMapper
Description
Class holds all necessary information for mapping fields associated with
dynamicTopoFvMesh and fvMesh.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
topoMapper.C
\*---------------------------------------------------------------------------*/
#ifndef topoMapper_H
#define topoMapper_H
#include "fluxCorrector.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
class topoCellMapper;
class topoSurfaceMapper;
class topoBoundaryMeshMapper;
class mapPolyMesh;
/*---------------------------------------------------------------------------*\
Class topoMapper Declaration
\*---------------------------------------------------------------------------*/
class topoMapper
{
// Private data
//- Reference to fvMesh
const fvMesh& mesh_;
//- Reference to the options dictionary
const dictionary& dict_;
// Demand-driven private data
//- Cell mapper
mutable autoPtr<topoCellMapper> cellMap_;
//- Surface mapper
mutable autoPtr<topoSurfaceMapper> surfaceMap_;
//- Boundary mapper
mutable autoPtr<topoBoundaryMeshMapper> boundaryMap_;
//- Flux corrector
mutable autoPtr<fluxCorrector> fluxCorrector_;
// Stored gradients for mapping
mutable HashTable<autoPtr<volVectorField> > sGrads_;
mutable HashTable<autoPtr<volTensorField> > vGrads_;
//- Geometric information on the old mesh
mutable scalarField* cellVolumesPtr_;
mutable vectorField* cellCentresPtr_;
mutable PtrList<scalarField> patchAreasPtr_;
mutable PtrList<vectorField> patchCentresPtr_;
//- Intersection volume weights
mutable List<scalarField> faceWeights_;
mutable List<scalarField> cellWeights_;
//- Intersection centre weights
mutable List<vectorField> faceCentres_;
mutable List<vectorField> cellCentres_;
// Private Member Functions
//- Disallow default bitwise copy construct
topoMapper(const topoMapper&);
//- Disallow default bitwise assignment
void operator=(const topoMapper&);
// Store gradients of volFields on the mesh
// prior to topology changes
template <class Type, class gradType>
void storeGradients
(
HashTable<autoPtr<gradType> >& gradTable
) const;
//- Fetch the gradient field
template <class Type>
const Type& gradient(const word& name) const;
//- Store gradients prior to mesh reset
void storeGradients() const;
//- Set geometric information
void storeGeometry() const;
public:
// Constructors
//- Construct from mesh and dictionary
topoMapper(const fvMesh& mesh, const dictionary& dict)
:
mesh_(mesh),
dict_(dict),
cellMap_(NULL),
surfaceMap_(NULL),
boundaryMap_(NULL),
fluxCorrector_(fluxCorrector::New(mesh, dict)),
cellVolumesPtr_(NULL),
cellCentresPtr_(NULL),
patchAreasPtr_(mesh.boundary().size()),
patchCentresPtr_(mesh.boundary().size())
{}
// Destructor
~topoMapper();
// Member Functions
//- Return reference to the mesh
const fvMesh& mesh() const;
//- Return reference to objectRegistry storing fields.
const objectRegistry& thisDb() const;
//- Set mapping information
void setMapper(const mapPolyMesh& mpm) const;
//- Set face weighting information
void setFaceWeights
(
const Xfer<List<scalarField> >& weights,
const Xfer<List<vectorField> >& centres
) const;
//- Set cell weighting information
void setCellWeights
(
const Xfer<List<scalarField> >& weights,
const Xfer<List<vectorField> >& centres
) const;
//- Fetch face weights
const List<scalarField>& faceWeights() const;
//- Fetch cell weights
const List<scalarField>& cellWeights() const;
//- Fetch face centres
const List<vectorField>& faceCentres() const;
//- Fetch cell centres
const List<vectorField>& cellCentres() const;
//- Store mesh information for the mapping stage
void storeMeshInformation() const;
//- Return stored cell volume information
const scalarField& cellVolumes() const;
//- Return stored cell centre information
const vectorField& internalCentres() const;
//- Return stored patch area information
const scalarField& patchAreas(const label i) const;
//- Return stored patch centre information
const vectorField& patchCentres(const label i) const;
// Conservatively map all volFields in the registry
template <class Type>
void conservativeMapVolFields() const;
// Conservatively map all surfaceFields in the registry
template <class Type>
void conservativeMapSurfaceFields() const;
//- Correct fluxes after topology change
void correctFluxes() const;
//- Return volume mapper
const topoCellMapper& volMap() const;
//- Return surface mapper
const topoSurfaceMapper& surfaceMap() const;
//- Return boundary mapper
const topoBoundaryMeshMapper& boundaryMap() const;
//- Return flux-corrector
const fluxCorrector& surfaceFluxCorrector() const;
//- Clear out member data
void clear();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#ifdef NoRepository
#include "topoMapper.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,776 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
topoPatchMapper
Description
Implementation of the topoPatchMapper class
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
#include "IOmanip.H"
#include "topoMapper.H"
#include "mapPolyMesh.H"
#include "topoPatchMapper.H"
#include "demandDrivenData.H"
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
//- Clear out local storage
void topoPatchMapper::clearOut()
{
deleteDemandDrivenData(directAddrPtr_);
deleteDemandDrivenData(interpolationAddrPtr_);
deleteDemandDrivenData(weightsPtr_);
deleteDemandDrivenData(insertedFaceLabelsPtr_);
deleteDemandDrivenData(insertedFaceIndexMapPtr_);
deleteDemandDrivenData(insertedFaceAddressingPtr_);
deleteDemandDrivenData(areasPtr_);
deleteDemandDrivenData(centresPtr_);
}
//- Calculate the insertedFaceLabels list
void topoPatchMapper::calcInsertedFaceAddressing() const
{
if (insertedFaceLabelsPtr_ || insertedFaceAddressingPtr_)
{
FatalErrorIn
(
"void topoPatchMapper::calcInsertedFaceAddressing() const"
) << " Inserted labels has already been calculated."
<< abort(FatalError);
}
// Information from the old patch
const label oldPatchSize = mpm_.oldPatchSizes()[patch_.index()];
const label oldPatchStart = mpm_.oldPatchStarts()[patch_.index()];
const label oldPatchEnd = oldPatchStart + oldPatchSize;
// Allocate for inserted face labels and addressing
label nInsertedFaces = 0;
insertedFaceLabelsPtr_ = new labelList(patchSize(), -1);
labelList& insertedFaces = *insertedFaceLabelsPtr_;
insertedFaceIndexMapPtr_ = new labelList(patchSize(), -1);
labelList& insertedFacesMap = *insertedFaceIndexMapPtr_;
insertedFaceAddressingPtr_ = new labelListList(patchSize(), labelList(0));
labelListList& insertedAddressing = *insertedFaceAddressingPtr_;
// Fetch current patch start
label pStart = patch_.patch().start();
// Fetch the current boundary
const polyBoundaryMesh& boundary = mpm_.mesh().boundaryMesh();
// Loop through the facesFromFaces map, and ensure that
// inserted faces are only mapped from faces on the same patch.
const List<objectMap>& fff = mpm_.facesFromFacesMap();
forAll(fff, objectI)
{
const objectMap& fffI = fff[objectI];
// Only pick boundary faces in this patch
if (boundary.whichPatch(fffI.index()) == patch_.index())
{
if (fffI.masterObjects().empty())
{
FatalErrorIn
(
"void topoPatchMapper::"
"calcInsertedFaceAddressing() const"
) << " Mapping for inserted boundary face is incorrect."
<< " Found an empty masterObjects list."
<< nl << " Face: " << fffI.index()
<< nl << " Patch: " << patch_.name()
<< abort(FatalError);
}
else
{
// Make an entry for the inserted label,
// and renumber addressing to patch.
insertedFaces[nInsertedFaces] = fffI.index() - pStart;
// Add a mapping entry for facesFromFaces indices
insertedFacesMap[nInsertedFaces] = objectI;
// Make an entry for addressing
labelList& addr = insertedAddressing[nInsertedFaces];
// Renumber addressing to patch.
// Also, check mapping for hits into
// other patches / internal faces.
addr = fffI.masterObjects();
forAll(addr, faceI)
{
if
(
addr[faceI] >= oldPatchStart
&& addr[faceI] < oldPatchEnd
)
{
addr[faceI] -= oldPatchStart;
}
else
{
FatalErrorIn
(
"void topoPatchMapper::"
"calcInsertedFaceAddressing() const"
)
<< "Addressing into another patch is not allowed."
<< nl << " Patch face index: " << faceI
<< nl << " addr[faceI]: " << addr[faceI]
<< nl << " oldPatchStart: " << oldPatchStart
<< nl << " oldPatchSize: " << oldPatchSize
<< nl << " oldPatchEnd: " << oldPatchEnd
<< abort(FatalError);
}
}
nInsertedFaces++;
}
}
}
// Shorten inserted faces to actual size
insertedFaces.setSize(nInsertedFaces);
insertedFacesMap.setSize(nInsertedFaces);
insertedAddressing.setSize(nInsertedFaces);
}
//- Calculate addressing for mapping
void topoPatchMapper::calcAddressing() const
{
if (directAddrPtr_ || interpolationAddrPtr_)
{
FatalErrorIn("void topoPatchMapper::calcAddressing() const")
<< "Addressing already calculated."
<< abort(FatalError);
}
// Information from the old patch
const label oldPatchSize = mpm_.oldPatchSizes()[patch_.index()];
const label oldPatchStart = mpm_.oldPatchStarts()[patch_.index()];
const label oldPatchEnd = oldPatchStart + oldPatchSize;
// Assemble the maps: slice to patch
if (direct())
{
// Direct mapping - slice to size
directAddrPtr_ = new labelList(patch_.patch().patchSlice(mpm_.faceMap()));
labelList& addr = *directAddrPtr_;
// Shift to local patch indices.
// Also, check mapping for hits into other patches / internal faces.
forAll (addr, faceI)
{
if
(
addr[faceI] >= oldPatchStart
&& addr[faceI] < oldPatchEnd
)
{
addr[faceI] -= oldPatchStart;
}
else
{
FatalErrorIn
(
"void topoPatchMapper::calcAddressing() const"
)
<< "Addressing into another patch is not allowed."
<< nl << " Patch face index: " << faceI
<< nl << " addr[faceI]: " << addr[faceI]
<< nl << " oldPatchStart: " << oldPatchStart
<< nl << " oldPatchSize: " << oldPatchSize
<< nl << " oldPatchEnd: " << oldPatchEnd
<< abort(FatalError);
}
}
}
else
{
// Interpolative addressing
interpolationAddrPtr_ = new labelListList(patchSize(), labelList(0));
labelListList& addr = *interpolationAddrPtr_;
// Fetch the list of inserted faces / addressing
const labelList& insertedFaces = insertedObjectLabels();
const labelListList& insertedAddressing = insertedFaceAddressing();
// Make entries
forAll(insertedFaces, faceI)
{
addr[insertedFaces[faceI]] = insertedAddressing[faceI];
}
// Do mapped faces. Note that this can already be set by insertedFaces
// so check if addressing size still zero.
const labelList& fm = patch_.patch().patchSlice(mpm_.faceMap());
forAll(fm, faceI)
{
if (fm[faceI] > -1 && addr[faceI].size() == 0)
{
// Mapped from a single face
label oldFace = fm[faceI];
if
(
oldFace >= oldPatchStart
&& oldFace < oldPatchEnd
)
{
oldFace -= oldPatchStart;
}
else
{
FatalErrorIn
(
"void topoPatchMapper::calcAddressing() const"
)
<< "Addressing into another patch is not allowed."
<< nl << " Patch face index: " << faceI
<< nl << " faceMap[faceI]: " << oldFace
<< nl << " oldPatchStart: " << oldPatchStart
<< nl << " oldPatchSize: " << oldPatchSize
<< nl << " oldPatchEnd: " << oldPatchEnd
<< abort(FatalError);
}
addr[faceI] = labelList(1, oldFace);
}
}
// Check if we missed anything
forAll(addr, faceI)
{
if (addr[faceI].empty())
{
FatalErrorIn
(
"void topoPatchMapper::calcAddressing() const"
)
<< "Addressing is missing." << nl
<< " Patch face index: " << faceI << nl
<< " nInsertedFaces: " << insertedFaces.size() << nl
<< " faceMap: " << fm[faceI] << nl
<< " Patch: " << patch_.name() << nl
<< abort(FatalError);
}
}
}
}
//- Calculate inverse-distance weights for interpolative mapping
void topoPatchMapper::calcInverseDistanceWeights() const
{
if (weightsPtr_)
{
FatalErrorIn
(
"void topoPatchMapper::calcInverseDistanceWeights() const"
)
<< "Weights already calculated."
<< abort(FatalError);
}
// Fetch interpolative addressing
const labelListList& addr = addressing();
// Allocate memory
weightsPtr_ = new scalarListList(patchSize());
scalarListList& w = *weightsPtr_;
// Obtain face-centre information from old/new meshes
const vectorField& oldCentres = tMapper_.patchCentres(patch_.index());
const vectorField& newCentres = patch_.patch().faceCentres();
forAll(addr, faceI)
{
const labelList& mo = addr[faceI];
// Do mapped faces
if (mo.size() == 1)
{
w[faceI] = scalarList(1, 1.0);
}
else
{
// Map from masters, inverse-distance weights
scalar totalWeight = 0.0;
w[faceI] = scalarList(mo.size(), 0.0);
forAll (mo, oldFaceI)
{
w[faceI][oldFaceI] =
(
1.0/stabilise
(
magSqr
(
newCentres[faceI]
- oldCentres[mo[oldFaceI]]
),
VSMALL
)
);
totalWeight += w[faceI][oldFaceI];
}
// Normalize weights
scalar normFactor = (1.0/totalWeight);
forAll (mo, oldFaceI)
{
w[faceI][oldFaceI] *= normFactor;
}
}
}
}
//- Calculate intersection weights for conservative mapping
void topoPatchMapper::calcIntersectionWeightsAndCentres() const
{
if (areasPtr_ || centresPtr_)
{
FatalErrorIn
(
"void topoPatchMapper::"
"calcIntersectionWeightsAndCentres() const"
)
<< "Weights already calculated."
<< abort(FatalError);
}
// Fetch interpolative addressing
const labelListList& addr = addressing();
// Allocate memory
areasPtr_ = new List<scalarField>(patchSize(), scalarField(0));
List<scalarField>& a = *areasPtr_;
centresPtr_ = new List<vectorField>(patchSize(), vectorField(0));
List<vectorField>& x = *centresPtr_;
// Obtain stored face-centres
const vectorField& faceCentres = tMapper_.patchCentres(patch_.index());
// Fetch maps
const List<vectorField>& mapFaceCentres = tMapper_.faceCentres();
const List<scalarField>& mapFaceWeights = tMapper_.faceWeights();
// Fill in maps first
const labelList& insertedFaces = insertedObjectLabels();
const labelList& insertedFacesMap = insertedObjectMap();
forAll(insertedFaces, faceI)
{
a[insertedFaces[faceI]] = mapFaceWeights[insertedFacesMap[faceI]];
x[insertedFaces[faceI]] = mapFaceCentres[insertedFacesMap[faceI]];
}
// Now do mapped faces
forAll(addr, faceI)
{
const labelList& mo = addr[faceI];
// Check if this is indeed a mapped face
if (mo.size() == 1 && x[faceI].empty() && a[faceI].empty())
{
x[faceI] = vectorField(1, faceCentres[mo[0]]);
a[faceI] = scalarField(1, 1.0);
}
}
// Final check to ensure everything went okay
forAll(a, faceI)
{
if (a[faceI].empty())
{
FatalErrorIn
(
"void topoPatchMapper::"
"calcIntersectionWeightsAndCentres() const"
)
<< "Weights / centres is missing."
<< nl << " Patch face index: " << faceI
<< abort(FatalError);
}
}
}
const List<scalarField>& topoPatchMapper::intersectionWeights() const
{
if (direct())
{
FatalErrorIn
(
"const List<scalarField>& "
"topoPatchMapper::intersectionWeights() const"
) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError);
}
if (!areasPtr_)
{
calcIntersectionWeightsAndCentres();
}
return *areasPtr_;
}
const List<vectorField>& topoPatchMapper::intersectionCentres() const
{
if (direct())
{
FatalErrorIn
(
"const List<vectorField>& "
"topoPatchMapper::intersectionCentres() const"
) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError);
}
if (!centresPtr_)
{
calcIntersectionWeightsAndCentres();
}
return *centresPtr_;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
topoPatchMapper::topoPatchMapper
(
const fvPatch& patch,
const mapPolyMesh& mpm,
const topoMapper& mapper
)
:
patch_(patch),
mpm_(mpm),
tMapper_(mapper),
direct_(false),
directAddrPtr_(NULL),
interpolationAddrPtr_(NULL),
weightsPtr_(NULL),
insertedFaceLabelsPtr_(NULL),
insertedFaceIndexMapPtr_(NULL),
insertedFaceAddressingPtr_(NULL),
areasPtr_(NULL),
centresPtr_(NULL)
{
// Check for the possibility of direct mapping
if (insertedObjects())
{
direct_ = false;
}
else
{
direct_ = true;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
topoPatchMapper::~topoPatchMapper()
{
clearOut();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Return the polyPatch size
label topoPatchMapper::patchSize() const
{
return patch_.patch().size();
}
//- Return size
label topoPatchMapper::size() const
{
return patch_.size();
}
//- Return size before mapping
label topoPatchMapper::sizeBeforeMapping() const
{
if (patch_.type() == "empty")
{
return 0;
}
return mpm_.oldPatchSizes()[patch_.index()];
}
//- Is the mapping direct
bool topoPatchMapper::direct() const
{
return direct_;
}
//- Return direct addressing
const unallocLabelList& topoPatchMapper::directAddressing() const
{
if (!direct())
{
FatalErrorIn
(
"const unallocLabelList& "
"topoPatchMapper::directAddressing() const"
) << "Requested direct addressing for an interpolative mapper."
<< abort(FatalError);
}
if (!directAddrPtr_)
{
calcAddressing();
}
return *directAddrPtr_;
}
//- Return interpolation addressing
const labelListList& topoPatchMapper::addressing() const
{
if (direct())
{
FatalErrorIn
(
"const labelListList& "
"topoPatchMapper::addressing() const"
) << "Requested interpolative addressing for a direct mapper."
<< abort(FatalError);
}
if (!interpolationAddrPtr_)
{
calcAddressing();
}
return *interpolationAddrPtr_;
}
//- Return weights
const scalarListList& topoPatchMapper::weights() const
{
if (direct())
{
FatalErrorIn
(
"const scalarListList& "
"topoPatchMapper::weights() const"
) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError);
}
if (!weightsPtr_)
{
calcInverseDistanceWeights();
}
return *weightsPtr_;
}
//- Are there any inserted faces
bool topoPatchMapper::insertedObjects() const
{
return insertedObjectLabels().size();
}
//- Return list of inserted faces
const labelList& topoPatchMapper::insertedObjectLabels() const
{
if (!insertedFaceLabelsPtr_)
{
calcInsertedFaceAddressing();
}
return *insertedFaceLabelsPtr_;
}
//- Return addressing map for inserted faces
const labelList& topoPatchMapper::insertedObjectMap() const
{
if (!insertedFaceIndexMapPtr_)
{
calcInsertedFaceAddressing();
}
return *insertedFaceIndexMapPtr_;
}
//- Return addressing for inserted faces
const labelListList& topoPatchMapper::insertedFaceAddressing() const
{
if (!insertedFaceAddressingPtr_)
{
calcInsertedFaceAddressing();
}
return *insertedFaceAddressingPtr_;
}
//- Map the patch field
template <class Type>
void topoPatchMapper::mapPatchField
(
const word& fieldName,
Field<Type>& pF
) const
{
// To invoke inverse-distance weighting, use this:
// pF.autoMap(*this);
// Check for possibility of direct mapping
if (direct())
{
pF.autoMap(*this);
}
else
{
if (pF.size() != sizeBeforeMapping())
{
FatalErrorIn
(
"\n\n"
"void topoCellMapper::mapPatchField<Type>\n"
"(\n"
" const word& fieldName,\n"
" Field<Type>& iF\n"
") const\n"
) << "Incompatible size before mapping." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << pF.size() << nl
<< " map size: " << sizeBeforeMapping() << nl
<< abort(FatalError);
}
// Fetch addressing
const labelListList& pAddressing = addressing();
const List<scalarField>& wF = intersectionWeights();
// Compute the integral of the source field
Type intSource = sum(pF * tMapper_.patchAreas(patch_.index()));
// Copy the original field
Field<Type> fieldCpy(pF);
// Resize to current dimensions
pF.setSize(size());
// Map the patch field
forAll(pF, faceI)
{
const labelList& addr = pAddressing[faceI];
pF[faceI] = pTraits<Type>::zero;
// Accumulate area-weighted interpolate
forAll(addr, faceJ)
{
pF[faceI] +=
(
wF[faceI][faceJ] * fieldCpy[addr[faceJ]]
);
}
}
// Compute the integral of the target field
const polyPatch& ppI = mpm_.mesh().boundaryMesh()[patch_.index()];
Type intTarget = sum(pF * mag(ppI.faceAreas()));
if (polyMesh::debug)
{
int oldP = Info().precision();
// Compare the global integral
Info << " Field : " << fieldName
<< " Patch : " << ppI.name()
<< " integral errors : "
<< setprecision(15)
<< " source : " << mag(intSource)
<< " target : " << mag(intTarget)
<< " norm : "
<< (mag(intTarget - intSource) / (mag(intSource) + VSMALL))
<< setprecision(oldP)
<< endl;
}
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void topoPatchMapper::operator=(const topoPatchMapper& rhs)
{
// Check for assignment to self
if (this == &rhs)
{
FatalErrorIn("void topoPatchMapper::operator=")
<< "Attempted assignment to self"
<< abort(FatalError);
}
}
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,207 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
topoPatchMapper
Description
This object provides mapping and fill-in information for patch data
between the two meshes after the topological change. It is
constructed from mapPolyMesh.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
topoPatchMapper.C
\*---------------------------------------------------------------------------*/
#ifndef topoPatchMapper_H
#define topoPatchMapper_H
#include "morphFieldMapper.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
class topoMapper;
class mapPolyMesh;
/*---------------------------------------------------------------------------*\
Class topoPatchMapper Declaration
\*---------------------------------------------------------------------------*/
class topoPatchMapper
:
public fvPatchFieldMapper
{
// Private data
//- Reference to patch
const fvPatch& patch_;
//- Reference to mapPolyMesh
const mapPolyMesh& mpm_;
//- Reference to the topoMapper
const topoMapper& tMapper_;
//- Is the mapping direct
bool direct_;
// Demand-driven private data
//- Direct addressing
mutable labelList* directAddrPtr_;
//- Interpolated addressing
mutable labelListList* interpolationAddrPtr_;
//- Inverse-distance weights
mutable scalarListList* weightsPtr_;
//- Inserted faces
mutable labelList* insertedFaceLabelsPtr_;
//- Inserted face index map
mutable labelList* insertedFaceIndexMapPtr_;
//- Inserted face addressing
mutable labelListList* insertedFaceAddressingPtr_;
//- Interpolation areas
mutable List<scalarField>* areasPtr_;
//- Interpolation centres
mutable List<vectorField>* centresPtr_;
// Private Member Functions
//- Disallow default bitwise copy construct
topoPatchMapper(const topoPatchMapper&);
//- Disallow default bitwise assignment
void operator=(const topoPatchMapper&);
//- Calculate the insertedFace addressing list
void calcInsertedFaceAddressing() const;
//- Calculate addressing for mapping
void calcAddressing() const;
//- Calculate inverse-distance weights for interpolative mapping
void calcInverseDistanceWeights() const;
//- Calculate intersection weights for conservative mapping
void calcIntersectionWeightsAndCentres() const;
//- Return intersection area weights
const List<scalarField>& intersectionWeights() const;
//- Return intersection area centres
const List<vectorField>& intersectionCentres() const;
//- Clear out local storage
void clearOut();
public:
// Constructors
//- Construct from mapPolyMesh
topoPatchMapper
(
const fvPatch& patch,
const mapPolyMesh& mpm,
const topoMapper& mapper
);
// Destructor
virtual ~topoPatchMapper();
// Member Functions
//- Return the polyPatch size
label patchSize() const;
//- Return size
virtual label size() const;
//- Return size of field before mapping
virtual label sizeBeforeMapping() const;
//- Is the mapping direct
virtual bool direct() const;
//- Return direct addressing
virtual const unallocLabelList& directAddressing() const;
//- Return interpolated addressing
virtual const labelListList& addressing() const;
//- Return interpolaion weights
virtual const scalarListList& weights() const;
//- Are there any inserted faces
virtual bool insertedObjects() const;
//- Return list of inserted faces
virtual const labelList& insertedObjectLabels() const;
//- Return addressing map for inserted faces
const labelList& insertedObjectMap() const;
//- Return addressing for inserted faces
const labelListList& insertedFaceAddressing() const;
//- Map the patch field
template <class Type>
void mapPatchField
(
const word& fieldName,
Field<Type>& pF
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#ifdef NoRepository
#include "topoPatchMapper.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,364 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
topoSurfaceMapper
Description
Implementation of the topoSurfaceMapper class
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
#include "topoMapper.H"
#include "mapPolyMesh.H"
#include "topoSurfaceMapper.H"
#include "demandDrivenData.H"
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
//- Clear out local storage
void topoSurfaceMapper::clearOut()
{
deleteDemandDrivenData(directAddrPtr_);
deleteDemandDrivenData(interpolationAddrPtr_);
deleteDemandDrivenData(weightsPtr_);
deleteDemandDrivenData(insertedFaceLabelsPtr_);
}
//- Calculate the insertedFaceLabels list
void topoSurfaceMapper::calcInsertedFaceLabels() const
{
if (insertedFaceLabelsPtr_)
{
FatalErrorIn
(
"void topoSurfaceMapper::calcInsertedFaceLabels() const"
) << " Inserted labels has already been calculated."
<< abort(FatalError);
}
// Allocate for inserted face labels
label nInsertedFaces = 0;
insertedFaceLabelsPtr_ = new labelList(size(), -1);
labelList& insertedFaces = *insertedFaceLabelsPtr_;
label nIntFaces = size();
// Loop through the facesFromFaces map, and ensure that
// inserted internal faces are not mapped from any parent faces.
const List<objectMap>& fff = mpm_.facesFromFacesMap();
forAll(fff, objectI)
{
const objectMap& fffI = fff[objectI];
// Only pick internal faces
if (fffI.index() < nIntFaces)
{
insertedFaces[nInsertedFaces++] = fffI.index();
}
}
// Shorten inserted faces to actual size
insertedFaces.setSize(nInsertedFaces);
}
//- Calculate addressing for mapping
void topoSurfaceMapper::calcAddressing() const
{
if
(
directAddrPtr_
|| interpolationAddrPtr_
)
{
FatalErrorIn("void topoSurfaceMapper::calcAddressing() const")
<< "Addressing already calculated."
<< abort(FatalError);
}
if (direct())
{
// Direct addressing, no weights - slice to size
directAddrPtr_ =
(
new labelList
(
labelList::subList(mpm_.faceMap(), size())
)
);
labelList& addr = *directAddrPtr_;
// Loop through the addressing list,
// and set dummy addressing for newly created faces.
forAll(addr, faceI)
{
if (addr[faceI] < 0)
{
addr[faceI] = 0;
}
}
}
else
{
FatalErrorIn("void topoSurfaceMapper::calcAddressing() const")
<< " Request for interpolative addressing"
<< " on internal surface fields. The surfaceMapper"
<< " is not designed to do this. Interpolate volFields"
<< " to surfaceFields instead."
<< abort(FatalError);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
topoSurfaceMapper::topoSurfaceMapper
(
const mapPolyMesh& mpm,
const topoMapper& mapper
)
:
mesh_(mpm.mesh()),
mpm_(mpm),
tMapper_(mapper),
direct_(false),
directAddrPtr_(NULL),
interpolationAddrPtr_(NULL),
weightsPtr_(NULL),
insertedFaceLabelsPtr_(NULL)
{
// Calculate the insertedFaces list
calcInsertedFaceLabels();
// Set to direct mapping.
direct_ = true;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
topoSurfaceMapper::~topoSurfaceMapper()
{
clearOut();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Return size
label topoSurfaceMapper::size() const
{
return mesh_.nInternalFaces();
}
//- Return size before mapping
label topoSurfaceMapper::sizeBeforeMapping() const
{
return mpm_.nOldInternalFaces();
}
//- Is the mapping direct
bool topoSurfaceMapper::direct() const
{
return direct_;
}
//- Return direct addressing
const unallocLabelList& topoSurfaceMapper::directAddressing() const
{
if (!direct())
{
FatalErrorIn
(
"const unallocLabelList& "
"topoSurfaceMapper::directAddressing() const"
) << "Requested direct addressing for an interpolative mapper."
<< abort(FatalError);
}
if (!directAddrPtr_)
{
calcAddressing();
}
return *directAddrPtr_;
}
//- Return interpolation addressing
const labelListList& topoSurfaceMapper::addressing() const
{
if (direct())
{
FatalErrorIn
(
"const labelListList& "
"topoSurfaceMapper::addressing() const"
) << "Requested interpolative addressing for a direct mapper."
<< abort(FatalError);
}
if (!interpolationAddrPtr_)
{
calcAddressing();
}
return *interpolationAddrPtr_;
}
//- Return weights
const scalarListList& topoSurfaceMapper::weights() const
{
if (direct())
{
FatalErrorIn
(
"const scalarListList& "
"topoSurfaceMapper::weights() const"
) << "Requested interpolative weights for a direct mapper."
<< abort(FatalError);
}
if (!weightsPtr_)
{
calcAddressing();
}
return *weightsPtr_;
}
//- Are there any inserted faces
bool topoSurfaceMapper::insertedObjects() const
{
return insertedObjectLabels().size();
}
//- Return list of inserted faces
const labelList& topoSurfaceMapper::insertedObjectLabels() const
{
if (!insertedFaceLabelsPtr_)
{
calcInsertedFaceLabels();
}
return *insertedFaceLabelsPtr_;
}
//- Return flux flip map
const labelHashSet& topoSurfaceMapper::flipFaceFlux() const
{
return mpm_.flipFaceFlux();
}
//- Map the internal field
template <class Type>
void topoSurfaceMapper::mapInternalField
(
const word& fieldName,
Field<Type>& iF
) const
{
if (iF.size() != sizeBeforeMapping())
{
FatalErrorIn
(
"\n\n"
"void topoSurfaceMapper::mapInternalField<Type>\n"
"(\n"
" Field<Type>& iF\n"
") const\n"
) << "Incompatible size before mapping." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << iF.size() << nl
<< " map size: " << sizeBeforeMapping() << nl
<< abort(FatalError);
}
// Map the internal field
iF.autoMap(*this);
// Flip the flux
const labelList flipFaces = flipFaceFlux().toc();
forAll (flipFaces, i)
{
if (flipFaces[i] < iF.size())
{
iF[flipFaces[i]] *= -1.0;
}
else
{
FatalErrorIn
(
"\n\n"
"void topoSurfaceMapper::mapInternalField<Type>\n"
"(\n"
" Field<Type>& iF\n"
") const\n"
) << "Cannot flip boundary face fluxes." << nl
<< " Field: " << fieldName << nl
<< " Field size: " << iF.size() << nl
<< " Face flip index: " << flipFaces[i] << nl
<< abort(FatalError);
}
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void topoSurfaceMapper::operator=(const topoSurfaceMapper& rhs)
{
// Check for assignment to self
if (this == &rhs)
{
FatalErrorIn("void topoSurfaceMapper::operator=")
<< "Attempted assignment to self"
<< abort(FatalError);
}
}
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,176 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
topoSurfaceMapper
Description
This object provides mapping and fill-in information for internal
face data between the two meshes after topological changes.
It is constructed from mapPolyMesh.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
topoSurfaceMapper.C
\*---------------------------------------------------------------------------*/
#ifndef topoSurfaceMapper_H
#define topoSurfaceMapper_H
#include "morphFieldMapper.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
class topoMapper;
class mapPolyMesh;
/*---------------------------------------------------------------------------*\
Class topoSurfaceMapper Declaration
\*---------------------------------------------------------------------------*/
class topoSurfaceMapper
:
public morphFieldMapper
{
// Private data
//- Reference to polyMesh
const polyMesh& mesh_;
//- Reference to mapPolyMesh
const mapPolyMesh& mpm_;
//- Reference to the topoMapper
const topoMapper& tMapper_;
//- Is the mapping direct
bool direct_;
// Demand-driven private data
//- Direct addressing
mutable labelList* directAddrPtr_;
//- Interpolated addressing
mutable labelListList* interpolationAddrPtr_;
//- Inverse-distance weights
mutable scalarListList* weightsPtr_;
//- Inserted faces
mutable labelList* insertedFaceLabelsPtr_;
// Private Member Functions
//- Disallow default bitwise copy construct
topoSurfaceMapper(const topoSurfaceMapper&);
//- Disallow default bitwise assignment
void operator=(const topoSurfaceMapper&);
//- Calculate the insertedFaceLabels list
void calcInsertedFaceLabels() const;
//- Calculate addressing for mapping
void calcAddressing() const;
//- Clear out local storage
void clearOut();
public:
// Constructors
//- Construct from mapPolyMesh
topoSurfaceMapper
(
const mapPolyMesh& mpm,
const topoMapper& mapper
);
// Destructor
virtual ~topoSurfaceMapper();
// Member Functions
//- Return size
virtual label size() const;
//- Return size before mapping
virtual label sizeBeforeMapping() const;
//- Is the mapping direct
virtual bool direct() const;
//- Return direct addressing
virtual const unallocLabelList& directAddressing() const;
//- Return interpolation addressing
virtual const labelListList& addressing() const;
//- Return weights
virtual const scalarListList& weights() const;
//- Are there any inserted faces
virtual bool insertedObjects() const;
//- Return list of inserted faces
virtual const labelList& insertedObjectLabels() const;
//- Return flux flip map
const labelHashSet& flipFaceFlux() const;
//- Map the internal field
template <class Type>
void mapInternalField
(
const word& fieldName,
Field<Type>& iF
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#ifdef NoRepository
#include "topoSurfaceMapper.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,930 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
meshOps
Description
Various utility functions that perform mesh-related operations.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*----------------------------------------------------------------------------*/
#include "Time.H"
#include "meshOps.H"
#include "ListOps.H"
#include "Pstream.H"
#include "triFace.H"
#include "IOmanip.H"
#include "HashSet.H"
#include "polyMesh.H"
#include "triPointRef.H"
#include "tetPointRef.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
namespace Foam
{
namespace meshOps
{
// Utility method to build a hull of cells
// connected to the edge (for 2D simplical meshes)
void constructPrismHull
(
const label eIndex,
const UList<face>& faces,
const UList<cell>& cells,
const UList<label>& owner,
const UList<label>& neighbour,
const UList<labelList>& edgeFaces,
labelList& hullTriFaces,
labelList& hullCells
)
{
labelHashSet cellSet, triFaceSet;
// Obtain references
const labelList& eFaces = edgeFaces[eIndex];
// Loop through edgeFaces and add cells
forAll(eFaces, faceI)
{
label c0 = owner[eFaces[faceI]];
label c1 = neighbour[eFaces[faceI]];
if (!cellSet.found(c0))
{
// Add this cell
cellSet.insert(c0);
// Find associated triFaces and add them too
const cell& cC = cells[c0];
forAll(cC, faceJ)
{
const face& cF = faces[cC[faceJ]];
if ((cF.size() == 3) && !(triFaceSet.found(cC[faceJ])))
{
triFaceSet.insert(cC[faceJ]);
}
}
}
if (!cellSet.found(c1) && (c1 != -1))
{
// Add this cell
cellSet.insert(c1);
// Find associated triFaces and add them too
const cell& cC = cells[c1];
forAll(cC, faceJ)
{
const face& cF = faces[cC[faceJ]];
if ((cF.size() == 3) && !(triFaceSet.found(cC[faceJ])))
{
triFaceSet.insert(cC[faceJ]);
}
}
}
}
// Obtain lists from hashSets
hullCells = cellSet.toc();
hullTriFaces = triFaceSet.toc();
}
// Utility method to build a hull of cells (and faces)
// around an edge (for 3D simplical meshes)
void constructHull
(
const label eIndex,
const UList<face>& faces,
const UList<edge>& edges,
const UList<cell>& cells,
const UList<label>& owner,
const UList<label>& neighbour,
const UList<labelList>& faceEdges,
const UList<labelList>& edgeFaces,
const UList<labelList>& edgePoints,
labelList& hullEdges,
labelList& hullFaces,
labelList& hullCells,
labelListList& ringEntities
)
{
// [1] hullEdges is an ordered list of edge-labels around eIndex,
// but not connected to it.
// - Ordering is in the same manner as edgePoints.
// [2] hullFaces is an ordered list of face-labels connected to eIndex.
// - Ordering is in the same manner as edgePoints.
// [3] hullCells is an ordered list of cell-labels connected to eIndex.
// - For boundary hulls, the last cell label is -1
// [4] ringEntities are edges and faces connected to eIndex[0] and eIndex[1]
// - ringEntities[0]: edges connected to eIndex[0]
// - ringEntities[1]: faces connected to eIndex[0]
// - ringEntities[2]: edges connected to eIndex[1]
// - ringEntities[3]: faces connected to eIndex[1]
bool found;
label otherPoint = -1, nextPoint = -1;
// Obtain a reference to this edge, and its edgeFaces
const edge& edgeToCheck = edges[eIndex];
const labelList& eFaces = edgeFaces[eIndex];
const labelList& hullVertices = edgePoints[eIndex];
// Loop through all faces of this edge and add them to hullFaces
forAll(eFaces, faceI)
{
const face& faceToCheck = faces[eFaces[faceI]];
// Find the isolated point on this face,
// and compare it with hullVertices
meshOps::findIsolatedPoint
(
faceToCheck,
edgeToCheck,
otherPoint,
nextPoint
);
found = false;
forAll(hullVertices, indexI)
{
if (hullVertices[indexI] == otherPoint)
{
// Fill in the position of this face on the hull
hullFaces[indexI] = eFaces[faceI];
// Obtain edges connected to top and bottom
// vertices of edgeToCheck
const labelList& fEdges = faceEdges[hullFaces[indexI]];
forAll(fEdges, edgeI)
{
if
(
edges[fEdges[edgeI]]
== edge(edgeToCheck[0], otherPoint)
)
{
ringEntities[0][indexI] = fEdges[edgeI];
}
if
(
edges[fEdges[edgeI]]
== edge(edgeToCheck[1], otherPoint)
)
{
ringEntities[2][indexI] = fEdges[edgeI];
}
}
// Depending on the orientation of this face,
// fill in hull cell indices as well
if (nextPoint == edgeToCheck[0])
{
hullCells[indexI] = owner[eFaces[faceI]];
}
else
if (nextPoint == edgeToCheck[1])
{
hullCells[indexI] = neighbour[eFaces[faceI]];
}
else
{
// Something's terribly wrong
FatalErrorIn("void meshOps::constructHull()")
<< nl << " Failed to construct hull. "
<< nl << " Possibly not a tetrahedral mesh. "
<< abort(FatalError);
}
if (hullCells[indexI] != -1)
{
label nextI = hullVertices.fcIndex(indexI);
label nextHullPoint = hullVertices[nextI];
const cell& currCell = cells[hullCells[indexI]];
// Look for the ring-faces
forAll(currCell, faceI)
{
const face& cFace = faces[currCell[faceI]];
// Check if this face contains edgeToCheck[0]
if
(
triFace::compare
(
triFace(cFace),
triFace
(
edgeToCheck[0],
otherPoint,
nextHullPoint
)
)
)
{
ringEntities[1][indexI] = currCell[faceI];
}
// Check if this face contains edgeToCheck[1]
if
(
triFace::compare
(
triFace(cFace),
triFace
(
edgeToCheck[1],
nextHullPoint,
otherPoint
)
)
)
{
ringEntities[3][indexI] = currCell[faceI];
}
}
// Scan one the faces for the ring-edge
const labelList& rFaceEdges =
(
faceEdges[ringEntities[1][indexI]]
);
forAll(rFaceEdges, edgeI)
{
if
(
edges[rFaceEdges[edgeI]]
== edge(otherPoint,nextHullPoint)
)
{
hullEdges[indexI] = rFaceEdges[edgeI];
break;
}
}
}
// Done with this index. Break out.
found = true;
break;
}
}
// Throw an error if the point wasn't found
if (!found)
{
// Something's terribly wrong
FatalErrorIn("void meshOps::constructHull()")
<< " Failed to construct hull. " << nl
<< " edgeFaces connectivity is inconsistent. " << nl
<< " Edge: " << eIndex << ":: " << edgeToCheck << nl
<< " edgeFaces: " << eFaces << nl
<< " edgePoints: " << hullVertices
<< abort(FatalError);
}
}
}
// Given a set of points and edges, find the shortest path
// between the start and end point, using Dijkstra's algorithm.
// - Takes a Map of points and edges that use those points.
// - Edge weights are currently edge-lengths, but can easily be adapted.
// - Returns true if the endPoint was found by the algorithm.
// - The Map 'pi' returns a preceding point for every point in 'points'.
//
// Algorithm is inspired by:
// Renaud Waldura
// Dijkstra's Shortest Path Algorithm in Java
// http://renaud.waldura.com/
bool Dijkstra
(
const Map<point>& points,
const Map<edge>& edges,
const label startPoint,
const label endPoint,
Map<label>& pi
)
{
bool foundEndPoint = false;
// Set of unvisited (Q) / visited (S) points and distances (d)
labelHashSet Q, S;
Map<scalar> d;
// Initialize distances to large values
forAllConstIter(Map<point>, points, pIter)
{
d.insert(pIter.key(), GREAT);
}
// Invert edges to make a local pointEdges list
Map<labelList> localPointEdges;
forAllConstIter(Map<edge>, edges, eIter)
{
const edge& edgeToCheck = eIter();
forAll(edgeToCheck, pointI)
{
if (!localPointEdges.found(edgeToCheck[pointI]))
{
localPointEdges.insert(edgeToCheck[pointI], labelList(0));
}
meshOps::sizeUpList
(
eIter.key(),
localPointEdges[edgeToCheck[pointI]]
);
}
}
// Mark the startPoint as having the smallest distance
d[startPoint] = 0.0;
// Add the startPoint to the list of unvisited points
Q.insert(startPoint);
while (Q.size())
{
// Step 1: Find the node with the smallest distance from the start.
labelHashSet::iterator smallest = Q.begin();
for
(
labelHashSet::iterator iter = ++Q.begin();
iter != Q.end();
iter++
)
{
if (d[iter.key()] < d[smallest.key()])
{
smallest = iter;
}
}
label pointIndex = smallest.key();
scalar smallestDistance = d[pointIndex];
// Move to the visited points list
S.insert(pointIndex);
Q.erase(pointIndex);
// Step 2: Build a list of points adjacent to pointIndex
// but not in the visited list
DynamicList<label> adjacentPoints(10);
const labelList& pEdges = localPointEdges[pointIndex];
forAll(pEdges, edgeI)
{
const edge& edgeToCheck = edges[pEdges[edgeI]];
label otherPoint = edgeToCheck.otherVertex(pointIndex);
if (!S.found(otherPoint))
{
adjacentPoints.append(otherPoint);
}
}
// Step 3: Perform distance-based checks for adjacent points
forAll(adjacentPoints, pointI)
{
label adjPoint = adjacentPoints[pointI];
scalar distance =
(
mag(points[adjPoint] - points[pointIndex])
+ smallestDistance
);
// Check if the end-point has been touched.
if (adjPoint == endPoint)
{
foundEndPoint = true;
}
if (distance < d[adjPoint])
{
// Update to the shorter distance
d[adjPoint] = distance;
// Update the predecessor
if (pi.found(adjPoint))
{
pi[adjPoint] = pointIndex;
}
else
{
pi.insert(adjPoint, pointIndex);
}
// Add to the list of unvisited points
Q.insert(adjPoint);
}
}
}
return foundEndPoint;
}
// Select a list of elements from connectivity,
// and output to a VTK format
void writeVTK
(
const polyMesh& mesh,
const word& name,
const labelList& cList,
const label primitiveType,
const UList<point>& meshPoints,
const UList<edge>& edges,
const UList<face>& faces,
const UList<cell>& cells,
const UList<label>& owner
)
{
label nTotalCells = 0;
label nPoints = 0, nCells = 0;
// Estimate a size for points and cellPoints
pointField points(6*cList.size());
// Connectivity lists
labelListList cpList(cList.size());
// Create a map for local points
Map<label> pointMap, reversePointMap, reverseCellMap;
forAll(cList, cellI)
{
if (cList[cellI] < 0)
{
continue;
}
// Are we looking at points?
if (primitiveType == 0)
{
// Size the list
cpList[nCells].setSize(1);
cpList[nCells] = cList[cellI];
nTotalCells++;
}
// Are we looking at edges?
if (primitiveType == 1)
{
// Size the list
cpList[nCells].setSize(2);
const edge& tEdge = edges[cList[cellI]];
cpList[nCells][0] = tEdge[0];
cpList[nCells][1] = tEdge[1];
nTotalCells += 2;
}
// Are we looking at faces?
if (primitiveType == 2)
{
const face& tFace = faces[cList[cellI]];
if (tFace.size() == 3)
{
// Size the list
cpList[nCells].setSize(3);
// Write out in order
cpList[nCells][0] = tFace[0];
cpList[nCells][1] = tFace[1];
cpList[nCells][2] = tFace[2];
nTotalCells += 3;
}
else
if (tFace.size() == 4)
{
// Size the list
cpList[nCells].setSize(4);
// Write out in order
cpList[nCells][0] = tFace[0];
cpList[nCells][1] = tFace[1];
cpList[nCells][2] = tFace[2];
cpList[nCells][3] = tFace[3];
nTotalCells += 4;
}
}
// Are we looking at cells?
if (primitiveType == 3)
{
const cell& tCell = cells[cList[cellI]];
if (tCell.size() == 4)
{
// Point-ordering for tetrahedra
const face& baseFace = faces[tCell[0]];
const face& checkFace = faces[tCell[1]];
// Size the list
cpList[nCells].setSize(4);
// Get the fourth point
label apexPoint =
(
meshOps::findIsolatedPoint(baseFace, checkFace)
);
// Something's wrong with connectivity.
if (apexPoint == -1)
{
FatalErrorIn
(
"void writeVTK\n"
"(\n"
" const polyMesh& mesh,\n"
" const word& name,\n"
" const labelList& cList,\n"
" const label primitiveType,\n"
" const UList<point>& points,\n"
" const UList<edge>& edges,\n"
" const UList<face>& faces,\n"
" const UList<cell>& cells,\n"
" const UList<label>& owner\n"
") const\n"
)
<< "Cell: " << cList[cellI]
<< ":: " << tCell
<< " has inconsistent connectivity."
<< abort(FatalError);
}
// Write-out in order
label ownCell = owner[tCell[0]];
if (ownCell == cList[cellI])
{
cpList[nCells][0] = baseFace[2];
cpList[nCells][1] = baseFace[1];
cpList[nCells][2] = baseFace[0];
cpList[nCells][3] = apexPoint;
}
else
{
cpList[nCells][0] = baseFace[0];
cpList[nCells][1] = baseFace[1];
cpList[nCells][2] = baseFace[2];
cpList[nCells][3] = apexPoint;
}
nTotalCells += 4;
}
else
if (tCell.size() == 5)
{
// Point-ordering for wedge cells
label firstTriFace = -1;
// Size the list
cpList[nCells].setSize(6);
// Figure out triangle faces
forAll(tCell, faceI)
{
const face& cFace = faces[tCell[faceI]];
if (cFace.size() == 3)
{
if (firstTriFace == -1)
{
firstTriFace = tCell[faceI];
// Right-handedness is assumed here.
// Tri-faces are always on the boundary.
cpList[nCells][0] = cFace[0];
cpList[nCells][1] = cFace[1];
cpList[nCells][2] = cFace[2];
}
else
{
// Detect the three other points.
forAll(tCell, faceJ)
{
const face& nFace = faces[tCell[faceJ]];
if (nFace.size() == 4)
{
// Search for vertices on cFace
// in this face.
forAll(cFace, I)
{
label i = nFace.which(cFace[I]);
if (i != -1)
{
label p = nFace.prevLabel(i);
label n = nFace.nextLabel(i);
if (p == cpList[nCells][0])
{
cpList[nCells][3] = cFace[I];
}
if (p == cpList[nCells][1])
{
cpList[nCells][4] = cFace[I];
}
if (p == cpList[nCells][2])
{
cpList[nCells][5] = cFace[I];
}
if (n == cpList[nCells][0])
{
cpList[nCells][3] = cFace[I];
}
if (n == cpList[nCells][1])
{
cpList[nCells][4] = cFace[I];
}
if (n == cpList[nCells][2])
{
cpList[nCells][5] = cFace[I];
}
}
}
}
}
break;
}
}
}
nTotalCells += 6;
}
}
// Renumber to local ordering
forAll(cpList[nCells], pointI)
{
// Check if this point was added to the map
if (!pointMap.found(cpList[nCells][pointI]))
{
// Point was not found, so add it
points[nPoints] = meshPoints[cpList[nCells][pointI]];
// Update the map
pointMap.insert(cpList[nCells][pointI], nPoints);
reversePointMap.insert(nPoints, cpList[nCells][pointI]);
// Increment the number of points
nPoints++;
}
// Renumber it.
cpList[nCells][pointI] = pointMap[cpList[nCells][pointI]];
}
// Update the cell map.
reverseCellMap.insert(nCells, cList[cellI]);
nCells++;
}
// Finally write it out
meshOps::writeVTK
(
mesh,
name,
nPoints,
nCells,
nTotalCells,
points,
cpList,
primitiveType,
reversePointMap,
reverseCellMap
);
}
// Actual routine to write out the VTK file
void writeVTK
(
const polyMesh& mesh,
const word& name,
const label nPoints,
const label nCells,
const label nTotalCells,
const vectorField& points,
const labelListList& cpList,
const label primitiveType,
const Map<label>& reversePointMap,
const Map<label>& reverseCellMap
)
{
// Make the directory
fileName dirName(mesh.time().path()/"VTK"/mesh.time().timeName());
mkDir(dirName);
// Open stream for output
OFstream file(dirName/name+".vtk");
// Write out the header
file << "# vtk DataFile Version 2.0" << nl
<< name << ".vtk" << nl
<< "ASCII" << nl
<< "DATASET UNSTRUCTURED_GRID" << nl
<< "POINTS " << nPoints << " double" << nl;
for (label i = 0; i < nPoints; i++)
{
file << setprecision(15)
<< points[i].x() << ' '
<< points[i].y() << ' '
<< points[i].z() << ' '
<< nl;
}
file << "CELLS " << nCells << " " << nTotalCells + nCells << endl;
if (cpList.size())
{
forAll(cpList, i)
{
if (cpList[i].size())
{
file << cpList[i].size() << ' ';
forAll(cpList[i], j)
{
file << cpList[i][j] << ' ';
}
file << nl;
}
}
}
else
{
// List of points
for (label i = 0; i < nPoints; i++)
{
file << 1 << ' ' << i << nl;
}
}
file << "CELL_TYPES " << nCells << endl;
if (cpList.size())
{
forAll(cpList, i)
{
if (cpList[i].size() == 1)
{
// Vertex
file << "1" << nl;
}
if (cpList[i].size() == 2)
{
// Edge
file << "3" << nl;
}
if (cpList[i].size() == 3)
{
// Triangle face
file << "5" << nl;
}
if
(
(cpList[i].size() == 4) &&
(primitiveType == 2)
)
{
// Quad face
file << "9" << nl;
}
if
(
(cpList[i].size() == 4) &&
(primitiveType == 3)
)
{
// Tetrahedron
file << "10" << nl;
}
if (cpList[i].size() == 6)
{
// Wedge
file << "13" << nl;
}
}
}
else
{
// List of points
for (label i = 0; i < nPoints; i++)
{
// Vertex
file << '1' << nl;
}
}
// Write out indices for visualization.
if (reverseCellMap.size())
{
file << "CELL_DATA " << nCells << endl;
file << "FIELD CellFields 1" << endl;
file << "CellIds 1 " << nCells << " int" << endl;
for (label i = 0; i < nCells; i++)
{
file << reverseCellMap[i] << ' ';
}
file << endl;
}
// Write out indices for visualization.
if (reversePointMap.size())
{
file << "POINT_DATA " << nPoints << endl;
file << "FIELD PointFields 1" << endl;
file << "PointIds 1 " << nPoints << " int" << endl;
for (label i = 0; i < nPoints; i++)
{
file << reversePointMap[i] << ' ';
}
file << endl;
}
}
} // End namespace meshOps
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,311 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
meshOps
Description
Various utility functions that perform mesh-related operations.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
meshOpsI.H
meshOps.C
\*---------------------------------------------------------------------------*/
#ifndef meshOps_H
#define meshOps_H
#include "Map.H"
#include "line.H"
#include "point.H"
#include "label.H"
#include "scalar.H"
#include "HashSet.H"
#include "cellList.H"
#include "edgeList.H"
#include "faceList.H"
#include "triangle.H"
#include "vectorField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class polyMesh;
/*---------------------------------------------------------------------------*\
Namespace meshOps Declaration
\*---------------------------------------------------------------------------*/
namespace meshOps
{
// Compute the centre for a given face, using UList
inline void faceCentre
(
const face& faceToCheck,
const UList<vector>& points,
vector& centre
);
// Compute the normal for a given face, using UList
inline void faceNormal
(
const face& faceToCheck,
const UList<vector>& points,
vector& normal
);
// Method to find the common-edge between two faces.
inline bool findCommonEdge
(
const label first,
const label second,
const UList<labelList>& faceEdges,
label& common
);
// Method to find the interior/boundary faces
// for an input quad-face and adjacent triangle-prism cell.
inline void findPrismFaces
(
const label fIndex,
const label cIndex,
const UList<face>& faces,
const UList<cell>& cells,
const UList<label>& neighbour,
FixedList<face,2>& bdyf,
FixedList<label,2>& bidx,
FixedList<face,2>& intf,
FixedList<label,2>& iidx
);
// Utility method to build a hull of cells
// connected to the edge (for 2D simplical meshes)
void constructPrismHull
(
const label eIndex,
const UList<face>& faces,
const UList<cell>& cells,
const UList<label>& owner,
const UList<label>& neighbour,
const UList<labelList>& edgeFaces,
labelList& hullTriFaces,
labelList& hullCells
);
// Utility method to build a hull of cells (and faces)
// around an edge (for 3D simplical meshes)
void constructHull
(
const label eIndex,
const UList<face>& faces,
const UList<edge>& edges,
const UList<cell>& cells,
const UList<label>& owner,
const UList<label>& neighbour,
const UList<labelList>& faceEdges,
const UList<labelList>& edgeFaces,
const UList<labelList>& edgePoints,
labelList& hullEdges,
labelList& hullFaces,
labelList& hullCells,
labelListList& ringEntities
);
// Utility method to find the isolated point
// given two triangular faces.
inline label findIsolatedPoint
(
const face& baseFace,
const face& checkFace
);
// Method to find the isolated point on a triangular face
// that doesn't lie on the specified edge
inline void findIsolatedPoint
(
const face& f,
const edge& e,
label& ptIndex,
label& nextPtIndex
);
// Given a face and cell index, find the apex point for a tet cell.
inline label tetApexPoint
(
const label cIndex,
const label fIndex,
const UList<face> faces,
const UList<cell> cells
);
// Given a cell index, find the centroid / volume of a cell.
template <class T1, class T2>
inline bool cellCentreAndVolume
(
const label cIndex,
const UList<Vector<T1> >& points,
const UList<face>& faces,
const UList<cell>& cells,
const UList<label>& owner,
Vector<T2>& centre,
T2& volume,
bool checkClosed = false
);
// Determine whether a segment intersects a triangular face
template <class T>
inline bool segmentTriFaceIntersection
(
const triangle<Vector<T>, const Vector<T>&>& faceToCheck,
const line<Vector<T>, const Vector<T>&>& edgeToCheck,
const T& matchTol,
Vector<T>& intPoint
);
// Determine whether the particular point lies
// inside the given triangular face
template <class T>
inline bool pointInTriFace
(
const triangle<Vector<T>, const Vector<T>&>& faceToCheck,
const Vector<T>& checkPoint,
const T& matchTol,
bool testCoplanarity
);
// Dijkstra's algorithm for the shortest path problem
bool Dijkstra
(
const Map<point>& points,
const Map<edge>& edges,
const label startPoint,
const label endPoint,
Map<label>& pi
);
// Method to insert labels in a face, so that
// right-handedness is preserved.
template <class T>
inline void insertPointLabels
(
const Vector<T>& refNorm,
const Field<Vector<T> >& points,
const labelHashSet& pLabels,
face& modFace
);
// Method to insert a label between two labels in a list
inline void insertLabel
(
const label newLabel,
const label labelA,
const label labelB,
labelList& list
);
// Utility method to replace a label in a given list
inline void replaceLabel
(
const label original,
const label replacement,
labelList& list
);
// Utility method to size-up the list to include an item
template <class Type>
inline void sizeUpList
(
const Type& item,
List<Type>& list
);
// Utility method to size-down the list to remove an item
template <class Type>
inline void sizeDownList
(
const Type& item,
List<Type>& list
);
// Remove an item at a particular index in the list
template <class Type>
inline void removeIndex
(
const label index,
List<Type>& list
);
// Select a list of elements from connectivity,
// and output to a VTK format
void writeVTK
(
const polyMesh& mesh,
const word& name,
const labelList& cList,
const label primitiveType,
const UList<point>& meshPoints,
const UList<edge>& edges,
const UList<face>& faces,
const UList<cell>& cells,
const UList<label>& owner
);
// Actual routine to write out the VTK file
void writeVTK
(
const polyMesh& mesh,
const word& name,
const label nPoints,
const label nCells,
const label nTotalCells,
const vectorField& points,
const labelListList& cpList = labelListList(0),
const label primitiveType = 3,
const Map<label>& reversePointMap = Map<label>(),
const Map<label>& reverseCellMap = Map<label>()
);
} // End namespace meshOps
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "meshOpsI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,843 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
meshOps
Description
Various utility functions that perform mesh-related operations.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
\*---------------------------------------------------------------------------*/
#include "Pair.H"
#include "meshOps.H"
#include "ListOps.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace meshOps
{
// Compute the centroid for a given face, using UList
inline void faceCentre
(
const face& faceToCheck,
const UList<vector>& points,
vector& centre
)
{
// Reset to zero
centre = vector::zero;
vector a(vector::zero);
vector b(vector::zero);
vector c(vector::zero);
// If the face is a triangle, do a direct calculation
// to avoid round-off error-related problems
if (faceToCheck.size() == 3)
{
a = points[faceToCheck[0]];
b = points[faceToCheck[1]];
c = points[faceToCheck[2]];
centre = ((1.0 / 3.0) * (a + b + c));
return;
}
label nPoints = faceToCheck.size();
register label pI;
// Store the centre point in c
for (pI = 0; pI < nPoints; pI++)
{
c += points[faceToCheck[pI]];
}
c /= nPoints;
scalar ta = 0.0;
scalar sumA = 0.0;
vector ttc = vector::zero;
vector sumAc = vector::zero;
for (pI = 0; pI < nPoints; pI++)
{
a = points[faceToCheck[pI]];
b = points[faceToCheck[(pI + 1) % nPoints]];
// Calculate 3*triangle centre
ttc = (a + b + c);
// Calculate 2*triangle area
ta = mag((a - c)^(b - c));
sumA += ta;
sumAc += ta*ttc;
}
if (sumA > VSMALL)
{
centre = (sumAc/(3.0 * sumA));
}
else
{
centre = c;
}
}
// Compute the normal for a given face, using UList
inline void faceNormal
(
const face& faceToCheck,
const UList<vector>& points,
vector& normal
)
{
// Reset to zero
normal = vector::zero;
vector a(vector::zero);
vector b(vector::zero);
vector c(vector::zero);
// If the face is a triangle, do a direct calculation
// to avoid round-off error-related problems
if (faceToCheck.size() == 3)
{
a = points[faceToCheck[0]];
b = points[faceToCheck[1]];
c = points[faceToCheck[2]];
normal = (0.5 * ((b - a)^(c - a)));
return;
}
label nPoints = faceToCheck.size();
register label pI;
// Store the centre point in c
for (pI = 0; pI < nPoints; pI++)
{
c += points[faceToCheck[pI]];
}
c /= nPoints;
for (pI = 0; pI < nPoints; pI++)
{
a = points[faceToCheck[pI]];
b = points[faceToCheck[(pI + 1) % nPoints]];
normal += (0.5 * ((b - a)^(c - a)));
}
}
// Utility method to find the common edge between two faces.
inline bool findCommonEdge
(
const label first,
const label second,
const UList<labelList>& faceEdges,
label& common
)
{
bool found = false;
const labelList& fEi = faceEdges[first];
const labelList& fEj = faceEdges[second];
forAll(fEi, edgeI)
{
forAll(fEj, edgeJ)
{
if (fEi[edgeI] == fEj[edgeJ])
{
common = fEi[edgeI];
found = true;
break;
}
}
if (found)
{
break;
}
}
return found;
}
// Utility method to find the interior (quad) / boundary (tri) faces
// for an input quad-face and adjacent triangle-prism cell.
inline void findPrismFaces
(
const label fIndex,
const label cIndex,
const UList<face>& faces,
const UList<cell>& cells,
const UList<label>& neighbour,
FixedList<face,2>& bdyf,
FixedList<label,2>& bidx,
FixedList<face,2>& intf,
FixedList<label,2>& iidx
)
{
label indexO = 0, indexI = 0;
const cell& c = cells[cIndex];
forAll(c, i)
{
label faceIndex = c[i];
// Don't count the face under consideration
if (faceIndex != fIndex)
{
const face& fi = faces[faceIndex];
if (neighbour[faceIndex] == -1)
{
if (fi.size() == 3)
{
// Triangular face on the boundary
bidx[indexO] = faceIndex;
bdyf[indexO++] = fi;
}
else
{
// This seems to be a non-triangular face on the boundary
// Consider this as "interior" and move on
iidx[indexI] = faceIndex;
intf[indexI++] = fi;
}
}
else
{
// Face on the interior
iidx[indexI] = faceIndex;
intf[indexI++] = fi;
}
}
}
}
// Utility method to find the isolated point given two triangular faces.
// - Returns the point on checkFace that does not belong to baseFace.
inline label findIsolatedPoint
(
const face& baseFace,
const face& checkFace
)
{
// Get the fourth point
forAll(checkFace, pointI)
{
if
(
checkFace[pointI] != baseFace[0] &&
checkFace[pointI] != baseFace[1] &&
checkFace[pointI] != baseFace[2]
)
{
return checkFace[pointI];
}
}
return -1;
}
// Utility method to find the isolated point on a triangular face
// that doesn't lie on the specified edge. Also returns the point next to it.
inline void findIsolatedPoint
(
const face& f,
const edge& e,
label& ptIndex,
label& nextPtIndex
)
{
// Check the first point
if ( f[0] != e.start() && f[0] != e.end() )
{
ptIndex = f[0];
nextPtIndex = f[1];
return;
}
// Check the second point
if ( f[1] != e.start() && f[1] != e.end() )
{
ptIndex = f[1];
nextPtIndex = f[2];
return;
}
// Check the third point
if ( f[2] != e.start() && f[2] != e.end() )
{
ptIndex = f[2];
nextPtIndex = f[0];
return;
}
// This bit should never happen.
FatalErrorIn
(
"inline void meshOps::findIsolatedPoint"
"(const face&, const edge&, label&, label&)"
)
<< "Cannot find isolated point in face " << f << endl
<< " Using edge: " << e
<< abort(FatalError);
}
// Given a face and cell index, find the apex point for a tet cell.
inline label tetApexPoint
(
const label cIndex,
const label fIndex,
const UList<face> faces,
const UList<cell> cells
)
{
label apexPoint = -1;
bool foundApex = false;
const cell& cellToCheck = cells[cIndex];
const face& baseFace = faces[fIndex];
forAll(cellToCheck, faceI)
{
const face& faceToCheck = faces[cellToCheck[faceI]];
apexPoint = findIsolatedPoint(baseFace, faceToCheck);
if (foundApex > -1)
{
foundApex = true;
break;
}
}
if (!foundApex)
{
Info << "fIndex: " << fIndex << ":: " << faces[fIndex] << endl;
Info << "cIndex: " << cIndex << ":: " << cellToCheck << endl;
forAll(cellToCheck, faceI)
{
Info << '\t' << cellToCheck[faceI] << ":: "
<< faces[cellToCheck[faceI]] << endl;
}
FatalErrorIn
(
"\n\n"
"inline label meshOps::tetApexPoint\n"
"(\n"
" const label cIndex,\n"
" const label fIndex,\n"
" const UList<face> faces,\n"
" const UList<cell> cells\n"
")\n"
)
<< "Could not find an apex point in cell: " << cIndex
<< " given face: " << fIndex
<< abort(FatalError);
}
return apexPoint;
}
// Given a cell index, find the centroid / volume of a cell.
// - If checking is enabled, return whether the cell is closed
template <class T1, class T2>
inline bool cellCentreAndVolume
(
const label cIndex,
const UList<Vector<T1> >& points,
const UList<face>& faces,
const UList<cell>& cells,
const UList<label>& owner,
Vector<T2>& centre,
T2& volume,
bool checkClosed
)
{
// Reset inputs
volume = pTraits<T2>::zero;
centre = Vector<T2>::zero;
const cell& cellToCheck = cells[cIndex];
// Average face-centres to get an estimate centroid
Vector<T2> cEst(Vector<T2>::zero);
Vector<T2> fCentre(Vector<T2>::zero);
Vector<T2> fArea(Vector<T2>::zero);
Vector<T2> sumClosed(Vector<T2>::zero), sumMagClosed(Vector<T2>::zero);
const T2 three(3.0), four(4.0);
const T2 oneThird = (pTraits<T2>::one / three);
const T2 oneFourth = (pTraits<T2>::one / four);
const T2 threeFourths = (three / four);
forAll(cellToCheck, faceI)
{
const face& checkFace = faces[cellToCheck[faceI]];
if (checkFace.empty())
{
continue;
}
meshOps::faceCentre
(
checkFace,
points,
fCentre
);
cEst += fCentre;
}
cEst /= T2(cellToCheck.size());
forAll(cellToCheck, faceI)
{
const face& checkFace = faces[cellToCheck[faceI]];
if (checkFace.empty())
{
continue;
}
meshOps::faceNormal
(
checkFace,
points,
fArea
);
meshOps::faceCentre
(
checkFace,
points,
fCentre
);
// Flip if necessary
if (owner[cellToCheck[faceI]] != cIndex)
{
fArea *= -1.0;
}
// Calculate 3*face-pyramid volume
T2 pyr3Vol = fArea & (fCentre - cEst);
// Calculate face-pyramid centre
Vector<T2> pc = (threeFourths*fCentre) + (oneFourth*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 *= oneThird;
if (checkClosed)
{
// Check the sum across components
T2 closed = pTraits<T2>::zero;
for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
{
closed =
(
mag(sumClosed[cmpt]) / (sumMagClosed[cmpt] + VSMALL)
);
if (closed > 1e-06)
{
// Volume is invalid
return false;
}
}
}
// Volume is valid
return true;
}
// Determine whether a segment intersects a triangular face
template <class T>
inline bool segmentTriFaceIntersection
(
const triangle<Vector<T>, const Vector<T>&>& faceToCheck,
const line<Vector<T>, const Vector<T>&>& edgeToCheck,
const T& matchTol,
Vector<T>& intPoint
)
{
// Fetch references
const Vector<T>& p1 = edgeToCheck.start();
const Vector<T>& p2 = edgeToCheck.end();
const Vector<T>& a = faceToCheck.a();
const Vector<T>& b = faceToCheck.b();
const Vector<T>& c = faceToCheck.c();
// Define custom face-normal
Vector<T> n = (T(0.5) * ((b - a)^(c - a)));
n /= mag(n) + VSMALL;
// Compute uValue
T numerator = n & (a - p1);
T denominator = n & (p2 - p1);
// Check if the edge is parallel to the face
if (mag(denominator) < VSMALL)
{
return false;
}
T u = (numerator / denominator);
T tolerance = (matchTol * mag(p2 - p1));
// Check for intersection along line.
if ((u > tolerance) && (u < (pTraits<T>::one - tolerance)))
{
// Compute point of intersection
intPoint = p1 + u*(p2 - p1);
// Also make sure that intPoint lies within the face
if (pointInTriFace(faceToCheck, intPoint, matchTol, false))
{
return true;
}
}
// Failed to fall within edge-bounds, or within face
return false;
}
// Determine whether the particular point lies
// inside the given triangular face
template <class T>
inline bool pointInTriFace
(
const triangle<Vector<T>, const Vector<T>&>& faceToCheck,
const Vector<T>& cP,
const T& matchTol,
bool testCoplanarity
)
{
// Copy inputs
const Vector<T>& a = faceToCheck.a();
const Vector<T>& b = faceToCheck.b();
const Vector<T>& c = faceToCheck.c();
const T pointFive(0.5);
// Compute the normal
Vector<T> nf = (pointFive * ((b - a)^(c - a)));
if ( ((pointFive * ((b - a)^(cP - a))) & nf) < pTraits<T>::zero)
{
return false;
}
if ( ((pointFive * ((c - b)^(cP - b))) & nf) < pTraits<T>::zero)
{
return false;
}
if ( ((pointFive * ((a - c)^(cP - c))) & nf) < pTraits<T>::zero)
{
return false;
}
if (testCoplanarity)
{
Vector<T> ftp(a - cP);
// Normalize vectors
nf /= mag(nf) + VSMALL;
ftp /= mag(ftp) + VSMALL;
if (mag(ftp & nf) > matchTol)
{
return false;
}
}
// Passed test with all edges
return true;
}
// 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
)
{
// Need to configure a new face.
face newFace(modFace);
const T pointFive(0.5);
forAllConstIter(labelHashSet, pLabels, pIter)
{
forAll(newFace, pI)
{
label nI = newFace.fcIndex(pI);
const Vector<T>& a = points[newFace[pI]];
const Vector<T>& b = points[pIter.key()];
const Vector<T>& c = points[newFace[nI]];
// Compute the normal.
Vector<T> newNorm = (pointFive * ((b - a)^(c - a)));
if ((refNorm & newNorm) > pTraits<T>::zero)
{
// Insert the point.
meshOps::insertLabel
(
pIter.key(),
newFace[pI],
newFace[nI],
newFace
);
break;
}
}
}
// Take over storage
modFace.transfer(newFace);
}
// Method to insert a label between two labels in a list
// Assumes that all labels are unique.
inline void insertLabel
(
const label newLabel,
const label labelA,
const label labelB,
labelList& list
)
{
// Create a new list
bool found = false;
label origSize = list.size();
labelList newList(origSize + 1);
label index = 0, nextI = -1;
// Start a linear search
forAll(list, itemI)
{
newList[index++] = list[itemI];
nextI = list.fcIndex(itemI);
if
(
(
(list[itemI] == labelA && list[nextI] == labelB) ||
(list[itemI] == labelB && list[nextI] == labelA)
) &&
!found
)
{
found = true;
newList[index++] = newLabel;
}
}
if (!found)
{
FatalErrorIn
(
"inline void meshOps::insertLabel"
"(const label, const label, const label, labelList&)"
) << nl << "Cannot insert " << newLabel
<< " in list: " << list << nl
<< " Labels: "
<< labelA << " and " << labelB
<< " were not found in sequence."
<< abort(FatalError);
}
// Transfer the list
list.transfer(newList);
}
// Utility method to replace a label in a given list
inline void replaceLabel
(
const label original,
const label replacement,
labelList& list
)
{
label index = -1;
if ((index = findIndex(list, original)) > -1)
{
list[index] = replacement;
}
else
{
FatalErrorIn
(
"inline void label meshOps::replaceLabel"
"(const label, const label, labelList&)"
) << nl << "Cannot find " << original
<< " in list: " << list << nl
<< " Label: " << replacement
<< " was not used in replacement."
<< abort(FatalError);
}
}
// Utility method to size-up the list to include an item
template <class Type>
inline void sizeUpList
(
const Type& item,
List<Type>& list
)
{
list.setSize(list.size() + 1, item);
}
// Utility method to size-down the list to remove an item
template <class Type>
inline void sizeDownList
(
const Type& item,
List<Type>& list
)
{
label index = -1;
if ((index = findIndex(list, item)) > -1)
{
meshOps::removeIndex(index, list);
}
else
{
FatalErrorIn
(
"inline void meshOps::sizeDownList"
"(const Type& item, List<Type>& list)"
)
<< nl << "Item: " << item
<< " was not found in list. " << nl
<< " List: " << nl << list
<< abort(FatalError);
}
}
// Remove an item at a particular index in the list
template <class Type>
inline void removeIndex
(
const label index,
List<Type>& list
)
{
// Create a new list
List<Type> newList(list.size() - 1);
// Copy individual items
label n = 0;
forAll(list, itemI)
{
if (itemI == index)
{
continue;
}
newList[n++] = list[itemI];
}
// Overwrite
list.transfer(newList);
}
} // End namespace meshOps
} // End namespace Foam
// ************************************************************************* //

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,278 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
lengthScaleEstimator
Description
Utility class used to provide length-scale estimates at various
mesh locations. These estimates are based on specified boundary
conditions provided through dictionary entries.
Author
Sandeep Menon
University of Massachusetts Amherst
All rights reserved
SourceFiles
lengthScaleEstimator.C
lengthScaleEstimatorI.H
\*---------------------------------------------------------------------------*/
#ifndef lengthScaleEstimator_H
#define lengthScaleEstimator_H
#include "polyMesh.H"
#include "dictionary.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class lengthScaleEstimator Declaration
\*---------------------------------------------------------------------------*/
class lengthScaleEstimator
{
// Private data
//- Const reference to polyMesh
const polyMesh& mesh_;
//- Edge bisection/collapse criteria
scalar ratioMin_;
scalar ratioMax_;
scalar growthFactor_;
scalar minLengthScale_;
scalar maxLengthScale_;
scalar curvatureDeviation_;
//- Specific to proximity-based refinement
label spatialRes_;
boundBox proxBoundBox_;
labelListList proximityBins_;
//- Specific to mesh-slicing operations
scalar sliceThreshold_;
label sliceHoldOff_;
List<boundBox> sliceBoxes_;
// Buffers for parallel length-scale calculations
labelListList sendLblBuffer_;
labelListList recvLblBuffer_;
scalarListList sendSclBuffer_;
scalarListList recvSclBuffer_;
//- Sub-dictionary which specifies
// fixed length-scales for patches
dictionary fixedPatches_;
//- Sub-dictionary which specifies
// floating length-scales for patches
dictionary freePatches_;
//- Sub-dictionary which specifies
// patches with curvature-based length-scale
dictionary curvaturePatches_;
//- Sub-dictionary which specifies
// patches with proximity-based length-scale
dictionary proximityPatches_;
//- Patches for which edge-refinements
// are to be avoided
labelList noModPatchIDs_;
//- Field-based refinement
word field_;
scalar fieldLength_;
scalar lowerRefineLevel_;
scalar upperRefineLevel_;
//- Specify limits for refinement criteria
scalar meanScale_;
label maxRefineLevel_;
// Private Member Functions
// Check for legitimacy of patches
void checkPatches(const wordList& patchList) const;
// Prepare for proximity-based refinement, if necessary
void prepareProximityPatches();
// Perform spatial hashing on a set of points
static void spatialHash
(
const pointField& pointLocations,
const labelList& pointIndices,
const boundBox& box,
const label resolution,
labelListList& bins
);
// Send length-scale info across processors
void writeLengthScaleInfo
(
const labelList& cellLevels,
const scalarList& lengthScale
);
// Receive length-scale info across processors
void readLengthScaleInfo
(
const label level,
label& visitedCells,
labelList& cellLevels,
UList<scalar>& lengthScale,
labelHashSet& levelCells
) const;
public:
// Declare the name of the class and its debug switch
TypeName("lengthScaleEstimator");
// Constructors
//- Construct from polyMesh and dictionary
explicit lengthScaleEstimator
(
const polyMesh& mesh
);
// Destructor
virtual ~lengthScaleEstimator();
// Member Functions
//- Read edge refinement options from the dictionary
void readRefinementOptions
(
const dictionary& refineDict,
bool reRead = false,
bool mandatory = false
);
//- Set explicitly coupled patch information
void setCoupledPatches
(
const dictionary& coupledPatches
);
//- Calculate the length scale field
void calculateLengthScale(UList<scalar>& lengthScale);
//- Return refinement criteria
inline scalar ratioMin() const;
inline scalar ratioMax() const;
inline scalar growthFactor() const;
//- Limit length scale for surface-edges
inline void limitScale(scalar& scale) const;
//- Check if a particular patch is free-floating
// (i.e., gets its length-scale from the interior)
inline bool isFreePatch(const label pIndex) const;
//- Check if a particular patch is flagged
// for proximity-based refinement
inline bool isProximityPatch(const label pIndex) const;
//- Check if a particular patch is flagged
// for curvature-based refinement
inline bool isCurvaturePatch(const label pIndex) const;
//- Return reference curvature deviation
inline scalar curvatureDeviation() const
{
return curvatureDeviation_;
}
//- Check whether a particular point is too close
// to a previous mesh slice location
inline bool checkOldSlices(const vector& gCentre) const;
//- Add a boundBox to the existing set of sliceBoxes
inline void appendBox(const boundBox& bBox);
//- Clear the list of sliceBoxes
inline void clearBoxes();
//- Set an initial hold-off value
inline void setHoldOff(const label hVal)
{
sliceHoldOff_ = hVal;
}
//- Return the current holdOff value
inline label holdOff() const
{
return sliceHoldOff_;
}
//- Decrement the current holdOff value
inline void decrementHoldOff()
{
sliceHoldOff_--;
}
//- Check whether a particular patch permits refinement
inline bool checkRefinementPatch
(
const label pIndex
) const;
//- Return the appropriate length-scale for boundary face
inline scalar fixedLengthScale
(
const label fIndex,
const label pIndex,
bool usePolyMesh = false
) const;
//- Test for proximity to patch faces
inline bool testProximity
(
const vector& gCentre,
const vector& gNormal,
const scalar testStep,
label& proxFace,
scalar& proxDistance
) const;
};
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "lengthScaleEstimatorI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

View file

@ -0,0 +1,329 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline scalar lengthScaleEstimator::ratioMin() const
{
return ratioMin_;
}
inline scalar lengthScaleEstimator::ratioMax() const
{
return ratioMax_;
}
inline scalar lengthScaleEstimator::growthFactor() const
{
return growthFactor_;
}
//- Limit length scale for surface-edges
inline void lengthScaleEstimator::limitScale(scalar& scale) const
{
// Need to account for multiplication by ratioMin / ratioMax as well
scale = Foam::max(scale, minLengthScale_ / ratioMin_);
scale = Foam::min(scale, maxLengthScale_ / ratioMax_);
}
//- Check if a particular patch is free-floating
// (i.e., gets its length-scale from the interior)
inline bool lengthScaleEstimator::isFreePatch
(
const label pIndex
) const
{
const polyBoundaryMesh& boundary = mesh_.boundaryMesh();
if (freePatches_.found(boundary[pIndex].name()))
{
return true;
}
// Certain boundary types are also considered floating
const word& pType = boundary[pIndex].type();
if
(
(pType == "processor") ||
(pType == "cyclic") ||
(pType == "wedge") ||
(pType == "empty") ||
(pType == "symmetryPlane")
)
{
return true;
}
return false;
}
//- Check if a particular patch is flagged
// for proximity-based refinement
inline bool lengthScaleEstimator::isProximityPatch
(
const label pIndex
) const
{
if (proximityPatches_.found(mesh_.boundaryMesh()[pIndex].name()))
{
return true;
}
return false;
}
//- Check if a particular patch is flagged
// for curvature-based refinement
inline bool lengthScaleEstimator::isCurvaturePatch
(
const label pIndex
) const
{
if (curvaturePatches_.found(mesh_.boundaryMesh()[pIndex].name()))
{
return true;
}
return false;
}
//- Check whether a particular point is too close
// to a previous mesh slice location
inline bool lengthScaleEstimator::checkOldSlices
(
const vector& gCentre
) const
{
forAll(sliceBoxes_, boxI)
{
if (sliceBoxes_[boxI].contains(gCentre))
{
// Too close to another slice-point. Bail out.
return true;
}
}
return false;
}
//- Add a boundBox to the existing set of sliceBoxes
inline void lengthScaleEstimator::appendBox
(
const boundBox& bBox
)
{
label currentSize = sliceBoxes_.size();
sliceBoxes_.setSize(currentSize + 1);
sliceBoxes_[currentSize] = bBox;
}
//- Clear the list of sliceBoxes
inline void lengthScaleEstimator::clearBoxes()
{
sliceBoxes_.clear();
}
//- Check whether a particular patch permits refinement
// - Return true if not permitted
inline bool lengthScaleEstimator::checkRefinementPatch
(
const label pIndex
) const
{
if (pIndex < 0)
{
return false;
}
if (findIndex(noModPatchIDs_, pIndex) > -1)
{
return true;
}
// Refinement is allowed
return false;
}
// Return the fixed length-scale value for a boundary face
inline scalar lengthScaleEstimator::fixedLengthScale
(
const label fIndex,
const label pIndex,
bool usePolyMesh
) const
{
scalar scale = 0.0;
const polyBoundaryMesh& boundary = mesh_.boundaryMesh();
// Check fixed length-scale patches
// If the value is negative, average face length-scales.
if (fixedPatches_.found(boundary[pIndex].name()))
{
scalar dictValue =
(
fixedPatches_[boundary[pIndex].name()][0].scalarToken()
);
if (dictValue > 0.0)
{
return dictValue;
}
}
// Approximate a length-scale from face area
if (usePolyMesh)
{
scale = Foam::sqrt(2.0 * mag(mesh_.faceAreas()[fIndex]));
}
return scale;
}
// Test an edge / face for proximity with other faces on proximity patches
// and return the scalar distance to an oppositely-oriented face.
inline bool lengthScaleEstimator::testProximity
(
const vector& gCentre,
const vector& gNormal,
const scalar testStep,
label& proxFace,
scalar& proxDistance
) const
{
// Reset input
proxDistance = GREAT;
DynamicList<label> posIndices(20);
scalar minDeviation = -0.9;
label nD = spatialRes_, binSize = proximityBins_.size();
// Obtain min/max extents
const point& bMin = proxBoundBox_.min();
const point& bMax = proxBoundBox_.max();
// Extend bounding-box dimensions a bit to avoid edge-effects.
scalar ext = 0.02*(mag(bMax - bMin));
// Define an inverse grid-cell size.
scalar xL = nD/(bMax.x() - bMin.x() + ext);
scalar yL = nD/(bMax.y() - bMin.y() + ext);
scalar zL = nD/(bMax.z() - bMin.z() + ext);
// Reset the proximity face
proxFace = -1;
// Now take multiple steps in both normal directions,
// and add to the list of boxes to be checked.
for (scalar dir = -1.0; dir < 2.0; dir += 2.0)
{
for (scalar step = 0.0; step < 5.0*testStep; step += testStep)
{
// Hash the point-location
point p = (gCentre + (dir*step*gNormal)) - bMin;
label i = label(mag(::floor(p.x()*xL)));
label j = label(mag(::floor(p.y()*yL)));
label k = label(mag(::floor(p.z()*zL)));
label pos = mag(((k*nD*nD)+(j*nD)+i) % binSize);
if (findIndex(posIndices, pos) == -1)
{
posIndices.append(pos);
}
}
}
// Obtain old-mesh face geometry for reference.
const vectorField& faceAreas = mesh_.faceAreas();
const vectorField& faceCentres = mesh_.faceCentres();
forAll(posIndices, indexI)
{
const labelList& posBin = proximityBins_[posIndices[indexI]];
forAll(posBin, faceI)
{
// Step 1: Measure the distance to the face.
vector rFace = (faceCentres[posBin[faceI]] - gCentre);
scalar distance = mag(rFace);
// Step 2: Check if this face is oriented away from face / edge.
const vector& fNormal = faceAreas[posBin[faceI]];
scalar deviation = (gNormal & (fNormal/(mag(fNormal) + VSMALL)));
if
(
(deviation < minDeviation) &&
(distance < proxDistance)
)
{
// Update statistics
proxFace = posBin[faceI];
proxDistance = distance;
// minDeviation = deviation;
}
}
}
// Check for fall below threshold
if ((proxFace != -1) && (proxDistance < sliceThreshold_))
{
return true;
}
// Does not fall below threshold
return false;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -3,7 +3,7 @@ EXE_INC = \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/tetDecompositionMotionSolver/lnInclude \
-I$(LIB_SRC)/dynamicMesh/meshMotion/tetDecompositionMotionSolver/lnInclude \
-I$(LIB_SRC)/tetDecompositionFiniteElement/lnInclude \
$(WM_DECOMP_INC)