zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: foamyHexMesh: Enable reading and writing of the background meshes

Browse Source
This commit is contained in:
laurence
2013-06-14 16:41:14 +01:00
parent 9d9bc02121
commit dd502d1e05
19 changed files with 723 additions and 295 deletions
Download Patch File Download Diff File Expand all files Collapse all files

115
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C
View File

@ -26,6 +26,10 @@ License
#include "DelaunayMesh.H"
#include "labelPair.H"
#include "PrintTable.H"
#include "pointIOField.H"
#include "scalarIOField.H"
#include "labelIOField.H"
#include "pointConversion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -36,14 +40,121 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Triangulation>
Foam::DelaunayMesh<Triangulation>::DelaunayMesh()
Foam::DelaunayMesh<Triangulation>::DelaunayMesh(const Time& runTime)
:
Triangulation(),
vertexCount_(0),
cellCount_(0)
cellCount_(0),
runTime_(runTime)
{}
template<class Triangulation>
Foam::DelaunayMesh<Triangulation>::DelaunayMesh
(
const Time& runTime,
const word& meshName
)
:
Triangulation(),
vertexCount_(0),
cellCount_(0),
runTime_(runTime)
{
pointIOField pts
(
IOobject
(
"points",
runTime.timeName(),
meshName/polyMesh::meshSubDir,
runTime,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
)
);
labelIOField types
(
IOobject
(
"types",
runTime.timeName(),
meshName,
runTime,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
)
);
labelIOField indices
(
IOobject
(
"indices",
runTime.timeName(),
meshName,
runTime,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
)
);
labelIOField processorIndices
(
IOobject
(
"processorIndices",
runTime.timeName(),
meshName,
runTime,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
)
);
if (pts.headerOk())
{
forAll(pts, ptI)
{
Vertex_handle vh = this->insert(toPoint<Point>(pts[ptI]));
if (indices.headerOk())
{
vh->index() = indices[ptI];
vertexCount()++;
}
else
{
vh->index() = getNewVertexIndex();
}
if (processorIndices.headerOk())
{
vh->procIndex() = processorIndices[ptI];
}
else
{
vh->procIndex() = Pstream::myProcNo();
}
if (types.headerOk())
{
vh->type() =
static_cast<Foam::indexedVertexEnum::vertexType>
(
types[ptI]
);
}
else
{
vh->type() = Vb::vtUnassigned;
}
}
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Triangulation>

35
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H
View File

@ -43,6 +43,7 @@ SourceFiles
#include "boundBox.H"
#include "indexedVertex.H"
#include "CGALTriangulation3Ddefs.H"
#include "Time.H"
#include "autoPtr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -82,6 +83,12 @@ public:
FixedList<label, 2>::Hash<>
> labelPairHashSet;
typedef HashTable
<
label,
labelPair,
FixedList<label, 2>::Hash<>
> labelTolabelPairHashTable;
private:
@ -95,6 +102,9 @@ private:
// This allows a unique index to be assigned to each cell.
mutable label cellCount_;
//- Reference to Time
const Time& runTime_;
//- Spatial sort traits to use with a pair of point pointers and an int.
// Taken from a post on the CGAL lists: 2010-01/msg00004.html by
// Sebastien Loriot (Geometry Factory).
@ -159,7 +169,13 @@ public:
// Constructors
//- Construct from components
DelaunayMesh();
explicit DelaunayMesh(const Time& runTime);
DelaunayMesh
(
const Time& runTime,
const word& meshName
);
//- Destructor
@ -168,6 +184,14 @@ public:
// Member Functions
inline const Time& time() const;
inline void timeCheck
(
const string& description,
const bool check = true
) const;
inline label getNewVertexIndex() const;
inline label getNewCellIndex() const;
@ -177,6 +201,7 @@ public:
inline void resetCellCount();
inline label vertexCount() const;
inline label& vertexCount();
inline void resetVertexCount();
@ -209,12 +234,12 @@ public:
//- Create an fvMesh from the triangulation.
// The mesh is not parallel consistent - only used for viewing
autoPtr<fvMesh> createMesh
autoPtr<polyMesh> createMesh
(
const fileName& name,
const Time& runTime,
labelList& vertexMap,
labelList& cellMap
labelTolabelPairHashTable& vertexMap,
labelList& cellMap,
const bool writeDelaunayData = true
) const;
};

40
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshI.H
View File

@ -36,6 +36,40 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Triangulation>
inline const Foam::Time& Foam::DelaunayMesh<Triangulation>::time() const
{
return runTime_;
}
template<class Triangulation>
void Foam::DelaunayMesh<Triangulation>::timeCheck
(
const string& description,
const bool check
) const
{
if (check)
{
Info<< nl << "--- [ cpuTime "
<< time().elapsedCpuTime() << " s, "
<< "delta " << time().cpuTimeIncrement()<< " s";
if (description != word::null)
{
Info<< ", " << description << " ";
}
else
{
Info<< " ";
}
Info<< "] --- " << endl;
}
}
template<class Triangulation>
inline Foam::label Foam::DelaunayMesh<Triangulation>::getNewVertexIndex() const
{
@ -90,6 +124,12 @@ Foam::label Foam::DelaunayMesh<Triangulation>::vertexCount() const
return vertexCount_;
}
template<class Triangulation>
Foam::label& Foam::DelaunayMesh<Triangulation>::vertexCount()
{
return vertexCount_;
}
template<class Triangulation>
void Foam::DelaunayMesh<Triangulation>::resetVertexCount()

107
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshIO.C
View File

@ -28,6 +28,7 @@ License
#include "pointConversion.H"
#include "wallPolyPatch.H"
#include "processorPolyPatch.H"
#include "labelIOField.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -331,13 +332,13 @@ void Foam::DelaunayMesh<Triangulation>::printVertexInfo(Ostream& os) const
template<class Triangulation>
Foam::autoPtr<Foam::fvMesh>
Foam::autoPtr<Foam::polyMesh>
Foam::DelaunayMesh<Triangulation>::createMesh
(
const fileName& name,
const Time& runTime,
labelList& vertexMap,
labelList& cellMap
labelTolabelPairHashTable& vertexMap,
labelList& cellMap,
const bool writeDelaunayData
) const
{
pointField points(Triangulation::number_of_vertices());
@ -354,12 +355,54 @@ Foam::DelaunayMesh<Triangulation>::createMesh
List<DynamicList<face> > patchFaces(1, DynamicList<face>());
List<DynamicList<label> > patchOwners(1, DynamicList<label>());
vertexMap.setSize(vertexCount(), -1);
vertexMap.resize(vertexCount());
cellMap.setSize(Triangulation::number_of_finite_cells(), -1);
// Calculate pts and a map of point index to location in pts.
label vertI = 0;
labelIOField indices
(
IOobject
(
"indices",
time().timeName(),
name,
time(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Triangulation::number_of_vertices()
);
labelIOField types
(
IOobject
(
"types",
time().timeName(),
name,
time(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Triangulation::number_of_vertices()
);
labelIOField processorIndices
(
IOobject
(
"processorIndices",
time().timeName(),
name,
time(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Triangulation::number_of_vertices()
);
for
(
Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
@ -369,13 +412,20 @@ Foam::DelaunayMesh<Triangulation>::createMesh
{
if (!vit->farPoint())
{
vertexMap[vit->index()] = vertI;
vertexMap(labelPair(vit->index(), vit->procIndex())) = vertI;
points[vertI] = topoint(vit->point());
indices[vertI] = vit->index();
types[vertI] = static_cast<label>(vit->type());
processorIndices[vertI] = vit->procIndex();
vertI++;
}
}
points.setSize(vertI);
indices.setSize(vertI);
types.setSize(vertI);
processorIndices.setSize(vertI);
// Index the cells
label cellI = 0;
@ -391,6 +441,7 @@ Foam::DelaunayMesh<Triangulation>::createMesh
(
!cit->hasFarPoint()
&& !Triangulation::is_infinite(cit)
&& cit->real()
)
{
cellMap[cit->cellIndex()] = cellI++;
@ -424,7 +475,12 @@ Foam::DelaunayMesh<Triangulation>::createMesh
label c1I = Cb::ctFar;
bool c1Real = false;
if (!c1->hasFarPoint() && !Triangulation::is_infinite(c1))
if
(
!c1->hasFarPoint()
&& !Triangulation::is_infinite(c1)
&& c1->real()
)
{
c1I = cellMap[c1->cellIndex()];
c1Real = true;
@ -432,7 +488,12 @@ Foam::DelaunayMesh<Triangulation>::createMesh
label c2I = Cb::ctFar;
bool c2Real = false;
if (!c2->hasFarPoint() && !Triangulation::is_infinite(c2))
if
(
!c2->hasFarPoint()
&& !Triangulation::is_infinite(c2)
&& c2->real()
)
{
c2I = cellMap[c2->cellIndex()];
c2Real = true;
@ -451,10 +512,17 @@ Foam::DelaunayMesh<Triangulation>::createMesh
{
verticesOnTriFace[i] = vertexMap
[
c1->vertex
labelPair
(
Triangulation::vertex_triple_index(oppositeVertex, i)
)->index()
c1->vertex
(
Triangulation::vertex_triple_index(oppositeVertex, i)
)->index(),
c1->vertex
(
Triangulation::vertex_triple_index(oppositeVertex, i)
)->procIndex()
)
];
}
@ -524,15 +592,15 @@ Foam::DelaunayMesh<Triangulation>::createMesh
Info<< "Creating mesh" << endl;
autoPtr<fvMesh> meshPtr
autoPtr<polyMesh> meshPtr
(
new fvMesh
new polyMesh
(
IOobject
(
name,
runTime.timeName(),
runTime,
time().timeName(),
time(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
@ -565,7 +633,14 @@ Foam::DelaunayMesh<Triangulation>::createMesh
patches.setSize(nValidPatches);
meshPtr().addFvPatches(patches);
meshPtr().addPatches(patches);
if (writeDelaunayData)
{
indices.write();
types.write();
processorIndices.write();
}
Info<< "Mesh created" << endl;

19
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
View File

@ -126,9 +126,24 @@ Foam::DistributedDelaunayMesh<Triangulation>::buildMap
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Triangulation>
Foam::DistributedDelaunayMesh<Triangulation>::DistributedDelaunayMesh()
Foam::DistributedDelaunayMesh<Triangulation>::DistributedDelaunayMesh
(
const Time& runTime
)
:
DelaunayMesh<Triangulation>(),
DelaunayMesh<Triangulation>(runTime),
allBackgroundMeshBounds_()
{}
template<class Triangulation>
Foam::DistributedDelaunayMesh<Triangulation>::DistributedDelaunayMesh
(
const Time& runTime,
const word& meshName
)
:
DelaunayMesh<Triangulation>(runTime, meshName),
allBackgroundMeshBounds_()
{}

14
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.H
View File

@ -135,13 +135,25 @@ public:
// Constructors
//- Construct from components
DistributedDelaunayMesh();
explicit DistributedDelaunayMesh(const Time& runTime);
DistributedDelaunayMesh
(
const Time& runTime,
const word& meshName
);
//- Destructor
~DistributedDelaunayMesh();
// Queries
//- Use DelaunayMesh timeCheck function
using DelaunayMesh<Triangulation>::timeCheck;
// Member Functions
//- Build a mapDistribute for the supplied destination processor data

5
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
View File

@ -791,10 +791,11 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
(
IOobject
(
fvMesh::defaultRegion,
"backgroundMeshDecomposition",
runTime_.timeName(),
runTime_,
IOobject::MUST_READ
IOobject::MUST_READ,
IOobject::AUTO_WRITE
)
),
meshCutter_

171
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
View File

@ -27,7 +27,7 @@ License
#include "cellSizeAndAlignmentControls.H"
#include "pointIOField.H"
#include "scalarIOField.H"
#include "tensorIOField.H"
#include "triadIOField.H"
#include "tetrahedron.H"
#include "plane.H"
#include "transform.H"
@ -38,6 +38,8 @@ License
namespace Foam
{
defineTypeNameAndDebug(cellShapeControlMesh, 0);
word cellShapeControlMesh::meshSubDir = "cellShapeControlMesh";
}
@ -366,9 +368,89 @@ void Foam::cellShapeControlMesh::writeTriangulation()
Foam::cellShapeControlMesh::cellShapeControlMesh(const Time& runTime)
:
DistributedDelaunayMesh<CellSizeDelaunay>
(
runTime,
meshSubDir
),
runTime_(runTime),
defaultCellSize_(0.0)
{}
{
if (this->vertexCount())
{
fvMesh mesh
(
IOobject
(
meshSubDir,
runTime.timeName(),
runTime,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
)
);
if (mesh.nPoints() == this->vertexCount())
{
pointScalarField sizes
(
IOobject
(
"sizes",
runTime.timeName(),
meshSubDir,
runTime,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
pointMesh::New(mesh)
);
triadIOField alignments
(
IOobject
(
"alignments",
mesh.time().timeName(),
meshSubDir,
mesh.time(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
)
);
if
(
sizes.size() == this->vertexCount()
&& alignments.size() == this->vertexCount()
)
{
label count = 0;
for
(
Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
++vit
)
{
vit->targetCellSize() = sizes[count];
vit->alignment() = alignments[count];
count++;
}
}
else
{
FatalErrorIn
(
"Foam::cellShapeControlMesh::cellShapeControlMesh"
"(const Time&)"
) << "Cell size point field is not the same size as the "
<< "mesh."
<< abort(FatalError);
}
}
}
}
//Foam::triangulatedMesh::triangulatedMesh
@ -744,9 +826,7 @@ void Foam::cellShapeControlMesh::insertBoundingPoints
void Foam::cellShapeControlMesh::write() const
{
Info<< "Writing cell size and alignment mesh" << endl;
const fileName name("cellSizeAndAlignmentMesh");
Info<< "Writing " << meshSubDir << endl;
// Reindex the cells
label cellCount = 0;
@ -763,17 +843,16 @@ void Foam::cellShapeControlMesh::write() const
}
}
labelList vertexMap;
DelaunayMesh<CellSizeDelaunay>::labelTolabelPairHashTable vertexMap;
labelList cellMap;
autoPtr<fvMesh> meshPtr = DelaunayMesh<CellSizeDelaunay>::createMesh
autoPtr<polyMesh> meshPtr = DelaunayMesh<CellSizeDelaunay>::createMesh
(
name,
runTime_,
meshSubDir,
vertexMap,
cellMap
);
const fvMesh& mesh = meshPtr();
const polyMesh& mesh = meshPtr();
pointScalarField sizes
(
@ -781,7 +860,8 @@ void Foam::cellShapeControlMesh::write() const
(
"sizes",
mesh.time().timeName(),
mesh,
meshSubDir,
mesh.time(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
@ -789,7 +869,22 @@ void Foam::cellShapeControlMesh::write() const
scalar(0)
);
OFstream str(runTime_.path()/"alignments.obj");
triadIOField alignments
(
IOobject
(
"alignments",
mesh.time().timeName(),
meshSubDir,
mesh.time(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sizes.size()
);
// Write alignments
// OFstream str(runTime_.path()/"alignments.obj");
for
(
@ -801,35 +896,41 @@ void Foam::cellShapeControlMesh::write() const
if (!vit->farPoint())
{
// Populate sizes
sizes[vertexMap[vit->index()]] = vit->targetCellSize();
sizes[vertexMap[labelPair(vit->index(), vit->procIndex())]] =
vit->targetCellSize();
// Write alignments
const tensor& alignment = vit->alignment();
pointFromPoint pt = topoint(vit->point());
alignments[vertexMap[labelPair(vit->index(), vit->procIndex())]] =
vit->alignment();
if
(
alignment.x() == triad::unset[0]
|| alignment.y() == triad::unset[0]
|| alignment.z() == triad::unset[0]
)
{
Info<< "Bad alignment = " << vit->info();
vit->alignment() = tensor::I;
Info<< "New alignment = " << vit->info();
continue;
}
meshTools::writeOBJ(str, pt, alignment.x() + pt);
meshTools::writeOBJ(str, pt, alignment.y() + pt);
meshTools::writeOBJ(str, pt, alignment.z() + pt);
// // Write alignments
// const tensor& alignment = vit->alignment();
// pointFromPoint pt = topoint(vit->point());
//
// if
// (
// alignment.x() == triad::unset[0]
// || alignment.y() == triad::unset[0]
// || alignment.z() == triad::unset[0]
// )
// {
// Info<< "Bad alignment = " << vit->info();
//
// vit->alignment() = tensor::I;
//
// Info<< "New alignment = " << vit->info();
//
// continue;
// }
//
// meshTools::writeOBJ(str, pt, alignment.x() + pt);
// meshTools::writeOBJ(str, pt, alignment.y() + pt);
// meshTools::writeOBJ(str, pt, alignment.z() + pt);
}
}
mesh.write();
sizes.write();
alignments.write();
}

3
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.H
View File

@ -95,6 +95,9 @@ public:
//- Runtime type information
ClassName("cellShapeControlMesh");
//- Return the mesh sub-directory name (usually "cellShapeControlMesh")
static word meshSubDir;
// Constructors

9
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C
View File

@ -223,6 +223,13 @@ void Foam::controlMeshRefinement::initialMeshPopulation
const autoPtr<backgroundMeshDecomposition>& decomposition
)
{
if (shapeController_.shapeControlMesh().vertexCount() > 0)
{
// Mesh already populated.
Info<< "Cell size and alignment mesh already populated." << endl;
return;
}
autoPtr<boundBox> overallBoundBox;
// Need to pass in the background mesh decomposition so that can test if
@ -268,7 +275,7 @@ void Foam::controlMeshRefinement::initialMeshPopulation
controlFunction.initialVertices(pts, sizes, alignments);
Info<< " Got initial vertices list" << endl;
Info<< " Got initial vertices list of size " << pts.size() << endl;
List<Vb> vertices(pts.size());

30
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
View File

@ -981,7 +981,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
const dictionary& foamyHexMeshDict
)
:
DistributedDelaunayMesh<Delaunay>(),
DistributedDelaunayMesh<Delaunay>(runTime),
runTime_(runTime),
rndGen_(64293*Pstream::myProcNo()),
foamyHexMeshControls_(foamyHexMeshDict),
@ -1135,8 +1135,6 @@ void Foam::conformalVoronoiMesh::initialiseForMotion()
Foam::indexedVertexEnum::vtExternalFeaturePoint
);
}
//writeFixedPoints("fixedPointsStart.obj");
}
@ -1814,32 +1812,6 @@ void Foam::conformalVoronoiMesh::move()
if (time().outputTime())
{
writeMesh(time().timeName());
// label cellI = 0;
// for
// (
// Finite_cells_iterator cit = finite_cells_begin();
// cit != finite_cells_end();
// ++cit
// )
// {
// if
// (
// !cit->hasFarPoint()
// && !is_infinite(cit)
// )
// {
// cit->cellIndex() = cellI++;
// }
// }
//
// labelList vertexMap;
// labelList cellMap;
// autoPtr<fvMesh> tetMesh =
// createMesh("tetMesh", runTime_, vertexMap, cellMap);
//
// tetMesh().write();
//writeFixedPoints("fixedPointsStart_" + runTime_.timeName() + ".obj");
}
updateSizesAndAlignments(pointsToInsert);

2
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
View File

@ -2808,7 +2808,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
patchPointPairSlaves[patchI].transfer(patchPPSlaves[patchI]);
}
// if (foamyHexMeshControls().objOutput())
if (foamyHexMeshControls().objOutput())
{
Info<< "Writing processor interfaces" << endl;

16
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
View File

@ -337,11 +337,10 @@ inline void Foam::conformalVoronoiMesh::createPointPair
const Foam::point internalPt = surfPt - ppDistn;
const Foam::point externalPt = surfPt + ppDistn;
if
(
geometryToConformTo_.inside(internalPt)
&& geometryToConformTo_.outside(externalPt)
)
bool internalInside = geometryToConformTo_.inside(internalPt);
bool externalOutside = geometryToConformTo_.outside(externalPt);
if (internalInside && externalOutside)
{
pts.append
(
@ -369,11 +368,8 @@ inline void Foam::conformalVoronoiMesh::createPointPair
{
Info<< "Warning: point pair not inside/outside" << nl
<< " surfPt = " << surfPt << nl
<< " internal = "
<< internalPt << " " << geometryToConformTo_.inside(internalPt)
<< nl
<< " external = "
<< externalPt << " " << geometryToConformTo_.outside(externalPt)
<< " internal = " << internalPt << " " << internalInside << nl
<< " external = " << externalPt << " " << externalOutside
<< endl;
}
}

411
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
View File

@ -421,202 +421,241 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
}
}
if (foamyHexMeshControls().writeCellShapeControlMesh())
{
cellShapeControls().shapeControlMesh().write();
}
if (foamyHexMeshControls().writeBackgroundMeshDecomposition())
{
Info<< nl << "Writing " << "backgroundMeshDecomposition" << endl;
decomposition_().mesh().write();
}
if (foamyHexMeshControls().writeTetDualMesh())
{
// Determine map from Delaunay vertex to Dual mesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// From all Delaunay vertices to cell (positive index)
// or patch face (negative index)
labelList vertexToDualAddressing(number_of_vertices(), 0);
forAll(cellToDelaunayVertex, cellI)
label cellI = 0;
for
(
Finite_cells_iterator cit = finite_cells_begin();
cit != finite_cells_end();
++cit
)
{
label vertI = cellToDelaunayVertex[cellI];
if (vertexToDualAddressing[vertI] != 0)
if
(
!cit->hasFarPoint()
&& !is_infinite(cit)
)
{
FatalErrorIn("conformalVoronoiMesh::writeMesh(..)")
<< "Delaunay vertex " << vertI
<< " from cell " << cellI
<< " is already mapped to "
<< vertexToDualAddressing[vertI]
<< exit(FatalError);
}
vertexToDualAddressing[vertI] = cellI+1;
}
forAll(patchToDelaunayVertex, patchI)
{
const labelList& patchVertices = patchToDelaunayVertex[patchI];
forAll(patchVertices, i)
{
label vertI = patchVertices[i];
if (vertexToDualAddressing[vertI] > 0)
{
FatalErrorIn("conformalVoronoiMesh::writeMesh(..)")
<< "Delaunay vertex " << vertI
<< " from patch " << patchI
<< " local index " << i
<< " is already mapped to cell "
<< vertexToDualAddressing[vertI]-1
<< exit(FatalError);
}
// Vertex might be used by multiple faces. Which one to
// use? For now last one wins.
label dualFaceI = dualPatchStarts[patchI]+i;
vertexToDualAddressing[vertI] = -dualFaceI-1;
cit->cellIndex() = cellI++;
}
}
Info<< nl << "Writing " << "tetDualMesh" << endl;
// Calculate tet mesh addressing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DistributedDelaunayMesh<Delaunay>::labelTolabelPairHashTable vertexMap;
labelList cellMap;
autoPtr<polyMesh> tetMesh =
createMesh("tetDualMesh", vertexMap, cellMap);
pointField points;
labelList boundaryPts(number_of_finite_cells(), -1);
// From tet point back to Delaunay vertex index
labelList pointToDelaunayVertex;
faceList faces;
labelList owner;
labelList neighbour;
wordList patchTypes;
wordList patchNames;
PtrList<dictionary> patchDicts;
pointField cellCentres;
tetMesh().write();
calcTetMesh
(
points,
pointToDelaunayVertex,
faces,
owner,
neighbour,
patchTypes,
patchNames,
patchDicts
);
// Calculate map from tet points to dual mesh cells/patch faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelIOList pointDualAddressing
(
IOobject
(
"pointDualAddressing",
instance,
"tetDualMesh"/polyMesh::meshSubDir,
runTime_,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
UIndirectList<label>
(
vertexToDualAddressing,
pointToDelaunayVertex
)()
);
label pointI = findIndex(pointDualAddressing, -1);
if (pointI != -1)
{
WarningIn
(
"conformalVoronoiMesh::writeMesh\n"
"(\n"
" const fileName& instance,\n"
" bool filterFaces\n"
")\n"
) << "Delaunay vertex " << pointI
<< " does not have a corresponding dual cell." << endl;
}
Info<< "Writing map from tetDualMesh points to Voronoi mesh to "
<< pointDualAddressing.objectPath() << endl;
pointDualAddressing.write();
// Write tet points corresponding to the Voronoi cell/face centre
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
{
// Read Voronoi mesh
fvMesh mesh
(
IOobject
(
Foam::polyMesh::defaultRegion,
instance,
runTime_,
IOobject::MUST_READ
)
);
pointIOField dualPoints
(
IOobject
(
"dualPoints",
instance,
"tetDualMesh"/polyMesh::meshSubDir,
runTime_,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
points
);
forAll(pointDualAddressing, pointI)
{
label index = pointDualAddressing[pointI];
if (index > 0)
{
label cellI = index-1;
dualPoints[pointI] = mesh.cellCentres()[cellI];
}
else if (index < 0)
{
label faceI = -index-1;
if (faceI >= mesh.nInternalFaces())
{
dualPoints[pointI] = mesh.faceCentres()[faceI];
}
}
}
Info<< "Writing new tetDualMesh points mapped onto Voronoi mesh to "
<< dualPoints.objectPath() << endl
<< "Replace the polyMesh/points with these." << endl;
dualPoints.write();
}
Info<< nl << "Writing tetDualMesh to " << instance << endl;
PackedBoolList boundaryFacesToRemove;
writeMesh
(
"tetDualMesh",
instance,
points,
boundaryPts,
faces,
owner,
neighbour,
patchTypes,
patchNames,
patchDicts,
cellCentres,
boundaryFacesToRemove
);
// // Determine map from Delaunay vertex to Dual mesh
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// // From all Delaunay vertices to cell (positive index)
// // or patch face (negative index)
// labelList vertexToDualAddressing(number_of_vertices(), 0);
//
// forAll(cellToDelaunayVertex, cellI)
// {
// label vertI = cellToDelaunayVertex[cellI];
//
// if (vertexToDualAddressing[vertI] != 0)
// {
// FatalErrorIn("conformalVoronoiMesh::writeMesh(..)")
// << "Delaunay vertex " << vertI
// << " from cell " << cellI
// << " is already mapped to "
// << vertexToDualAddressing[vertI]
// << exit(FatalError);
// }
// vertexToDualAddressing[vertI] = cellI+1;
// }
//
// forAll(patchToDelaunayVertex, patchI)
// {
// const labelList& patchVertices = patchToDelaunayVertex[patchI];
//
// forAll(patchVertices, i)
// {
// label vertI = patchVertices[i];
//
// if (vertexToDualAddressing[vertI] > 0)
// {
// FatalErrorIn("conformalVoronoiMesh::writeMesh(..)")
// << "Delaunay vertex " << vertI
// << " from patch " << patchI
// << " local index " << i
// << " is already mapped to cell "
// << vertexToDualAddressing[vertI]-1
// << exit(FatalError);
// }
//
// // Vertex might be used by multiple faces. Which one to
// // use? For now last one wins.
// label dualFaceI = dualPatchStarts[patchI]+i;
// vertexToDualAddressing[vertI] = -dualFaceI-1;
// }
// }
//
//
// // Calculate tet mesh addressing
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// pointField points;
// labelList boundaryPts(number_of_finite_cells(), -1);
// // From tet point back to Delaunay vertex index
// labelList pointToDelaunayVertex;
// faceList faces;
// labelList owner;
// labelList neighbour;
// wordList patchTypes;
// wordList patchNames;
// PtrList<dictionary> patchDicts;
// pointField cellCentres;
//
// calcTetMesh
// (
// points,
// pointToDelaunayVertex,
// faces,
// owner,
// neighbour,
// patchTypes,
// patchNames,
// patchDicts
// );
//
//
//
// // Calculate map from tet points to dual mesh cells/patch faces
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// labelIOList pointDualAddressing
// (
// IOobject
// (
// "pointDualAddressing",
// instance,
// "tetDualMesh"/polyMesh::meshSubDir,
// runTime_,
// IOobject::NO_READ,
// IOobject::AUTO_WRITE,
// false
// ),
// UIndirectList<label>
// (
// vertexToDualAddressing,
// pointToDelaunayVertex
// )()
// );
//
// label pointI = findIndex(pointDualAddressing, -1);
// if (pointI != -1)
// {
// WarningIn
// (
// "conformalVoronoiMesh::writeMesh\n"
// "(\n"
// " const fileName& instance,\n"
// " bool filterFaces\n"
// ")\n"
// ) << "Delaunay vertex " << pointI
// << " does not have a corresponding dual cell." << endl;
// }
//
// Info<< "Writing map from tetDualMesh points to Voronoi mesh to "
// << pointDualAddressing.objectPath() << endl;
// pointDualAddressing.write();
//
//
//
// // Write tet points corresponding to the Voronoi cell/face centre
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// {
// // Read Voronoi mesh
// fvMesh mesh
// (
// IOobject
// (
// Foam::polyMesh::defaultRegion,
// instance,
// runTime_,
// IOobject::MUST_READ
// )
// );
// pointIOField dualPoints
// (
// IOobject
// (
// "dualPoints",
// instance,
// "tetDualMesh"/polyMesh::meshSubDir,
// runTime_,
// IOobject::NO_READ,
// IOobject::AUTO_WRITE,
// false
// ),
// points
// );
//
// forAll(pointDualAddressing, pointI)
// {
// label index = pointDualAddressing[pointI];
//
// if (index > 0)
// {
// label cellI = index-1;
// dualPoints[pointI] = mesh.cellCentres()[cellI];
// }
// else if (index < 0)
// {
// label faceI = -index-1;
// if (faceI >= mesh.nInternalFaces())
// {
// dualPoints[pointI] = mesh.faceCentres()[faceI];
// }
// }
// }
//
// Info<< "Writing tetDualMesh points mapped onto Voronoi mesh to "
// << dualPoints.objectPath() << endl
// << "Replace the polyMesh/points with these." << endl;
// dualPoints.write();
// }
//
//
// Info<< nl << "Writing tetDualMesh to " << instance << endl;
//
// PackedBoolList boundaryFacesToRemove;
// writeMesh
// (
// "tetDualMesh",
// instance,
// points,
// boundaryPts,
// faces,
// owner,
// neighbour,
// patchTypes,
// patchNames,
// patchDicts,
// cellCentres,
// boundaryFacesToRemove
// );
}
}

9
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cvControls/cvControls.C
View File

@ -247,6 +247,15 @@ Foam::cvControls::cvControls
}
writeTetDualMesh_ = Switch(filteringDict.lookup("writeTetDualMesh"));
writeCellShapeControlMesh_ =
Switch(filteringDict.lookup("writeCellShapeControlMesh"));
writeBackgroundMeshDecomposition_ =
(
Switch(filteringDict.lookup("writeBackgroundMeshDecomposition"))
&& Pstream::parRun()
);
}

10
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cvControls/cvControls.H
View File

@ -202,6 +202,10 @@ class cvControls
//- Write tet mesh at output time (it always writes the Voronoi)
Switch writeTetDualMesh_;
Switch writeCellShapeControlMesh_;
Switch writeBackgroundMeshDecomposition_;
// Private Member Functions
@ -335,6 +339,12 @@ public:
//- Write tetMesh at output time
inline Switch writeTetDualMesh() const;
//- Write cellShapeControlMesh at output time
inline Switch writeCellShapeControlMesh() const;
//- Write backgroundMeshDecomposition at output time
inline Switch writeBackgroundMeshDecomposition() const;
};

10
applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
View File

@ -215,5 +215,15 @@ inline Foam::Switch Foam::cvControls::writeTetDualMesh() const
return writeTetDualMesh_;
}
inline Foam::Switch Foam::cvControls::writeCellShapeControlMesh() const
{
return writeCellShapeControlMesh_;
}
inline Foam::Switch Foam::cvControls::writeBackgroundMeshDecomposition() const
{
return writeBackgroundMeshDecomposition_;
}
// ************************************************************************* //