ENH: cvMesh: split refinement and smoothing of background mesh

This commit is contained in:
laurence 2013-01-28 18:01:22 +00:00
parent 07cb342aeb
commit 53533a2da2
26 changed files with 1676 additions and 1017 deletions

View File

@ -29,6 +29,10 @@ $(cellSizeAndAlignmentControl)/searchableSurfaceControl/searchableSurfaceControl
cellShapeControl/cellAspectRatioControl/cellAspectRatioControl.C
cellShapeControl/smoothAlignmentSolver/smoothAlignmentSolver.C
cellShapeControl/controlMeshRefinement/controlMeshRefinement.C
/*cellSizeControlSurfaces/cellSizeControlSurfaces.C*/
cellSizeFunctions = cellSizeControlSurfaces/cellSizeFunction

View File

@ -31,213 +31,16 @@ License
#include "searchableSurfaceControl.H"
#include "cellSizeFunction.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(cellShapeControl, 0);
}
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template <class Triangulation, class Type>
Foam::tmp<Foam::Field<Type> > Foam::cellShapeControl::filterFarPoints
(
const Triangulation& mesh,
const Field<Type>& field
)
{
tmp<Field<Type> > tNewField(new Field<Type>(field.size()));
Field<Type>& newField = tNewField();
label added = 0;
label count = 0;
for
(
typename Triangulation::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
vit != mesh.finite_vertices_end();
++vit
)
{
if (vit->real())
{
newField[added++] = field[count];
}
count++;
}
newField.resize(added);
return tNewField;
}
template <class Triangulation>
Foam::autoPtr<Foam::mapDistribute> Foam::cellShapeControl::buildReferredMap
(
const Triangulation& mesh,
labelList& indices
)
{
globalIndex globalIndexing(mesh.vertexCount());
DynamicList<label> dynIndices(mesh.vertexCount()/10);
for
(
typename Triangulation::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
vit != mesh.finite_vertices_end();
++vit
)
{
if (vit->referred())
{
dynIndices.append
(
globalIndexing.toGlobal(vit->procIndex(), vit->index())
);
}
}
indices.transfer(dynIndices);
List<Map<label> > compactMap;
return autoPtr<mapDistribute>
(
new mapDistribute
(
globalIndexing,
indices,
compactMap
)
);
}
template <class Triangulation>
Foam::autoPtr<Foam::mapDistribute> Foam::cellShapeControl::buildMap
(
const Triangulation& mesh,
labelListList& pointPoints
)
{
pointPoints.setSize(mesh.vertexCount());
globalIndex globalIndexing(mesh.vertexCount());
for
(
typename Triangulation::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
vit != mesh.finite_vertices_end();
++vit
)
{
if (!vit->real())
{
continue;
}
std::list<typename Triangulation::Vertex_handle> adjVerts;
mesh.finite_adjacent_vertices(vit, std::back_inserter(adjVerts));
DynamicList<label> indices(adjVerts.size());
for
(
typename std::list<typename Triangulation::Vertex_handle>::
const_iterator adjVertI = adjVerts.begin();
adjVertI != adjVerts.end();
++adjVertI
)
{
typename Triangulation::Vertex_handle vh = *adjVertI;
if (!vh->farPoint())
{
indices.append
(
globalIndexing.toGlobal(vh->procIndex(), vh->index())
);
}
}
pointPoints[vit->index()].transfer(indices);
}
List<Map<label> > compactMap;
return autoPtr<mapDistribute>
(
new mapDistribute
(
globalIndexing,
pointPoints,
compactMap
)
);
}
template <class Triangulation>
Foam::tmp<Foam::triadField> Foam::cellShapeControl::buildAlignmentField
(
const Triangulation& mesh
)
{
tmp<triadField> tAlignments
(
new triadField(mesh.vertexCount(), triad::unset)
);
triadField& alignments = tAlignments();
for
(
typename Triangulation::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
vit != mesh.finite_vertices_end();
++vit
)
{
if (!vit->real())
{
continue;
}
alignments[vit->index()] = vit->alignment();
}
return tAlignments;
}
template <class Triangulation>
Foam::tmp<Foam::pointField> Foam::cellShapeControl::buildPointField
(
const Triangulation& mesh
)
{
tmp<pointField> tPoints
(
new pointField(mesh.vertexCount(), point(GREAT, GREAT, GREAT))
);
pointField& points = tPoints();
for
(
typename Triangulation::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
vit != mesh.finite_vertices_end();
++vit
)
{
if (!vit->real())
{
continue;
}
points[vit->index()] = topoint(vit->point());
}
return tPoints;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -261,7 +64,7 @@ Foam::cellShapeControl::cellShapeControl
(
runTime,
motionDict.subDict("shapeControlFunctions"),
geometryToConformTo,
geometryToConformTo_,
defaultCellSize_
)
{}
@ -300,16 +103,7 @@ Foam::scalar Foam::cellShapeControl::cellSize(const point& pt) const
scalar size = 0;
label nFarPoints = 0;
for (label pI = 0; pI < 4; ++pI)
{
if (ch->vertex(pI)->farPoint())
{
nFarPoints++;
}
}
if (shapeControlMesh_.is_infinite(ch))
if (shapeControlMesh_.dimension() < 3 || shapeControlMesh_.is_infinite(ch))
{
// if (nFarPoints != 0)
// {
@ -323,22 +117,25 @@ Foam::scalar Foam::cellShapeControl::cellSize(const point& pt) const
// }
// }
// Look up nearest point
cellShapeControlMesh::Vertex_handle nearV =
shapeControlMesh_.nearest_vertex
(
toPoint<cellShapeControlMesh::Point>(pt)
);
size = nearV->targetCellSize();
// Find nearest surface
size = sizeAndAlignment_.cellSize(pt);
}
else
{
label nFarPoints = 0;
for (label pI = 0; pI < 4; ++pI)
{
if (ch->vertex(pI)->farPoint())
{
nFarPoints++;
}
}
if (nFarPoints != 0)
{
for (label pI = 0; pI < 4; ++pI)
{
if (!ch->vertex(pI)->farPoint())
if (!ch->vertex(pI)->uninitialised())
{
size = ch->vertex(pI)->targetCellSize();
return size;
@ -368,54 +165,58 @@ Foam::tensor Foam::cellShapeControl::cellAlignment(const point& pt) const
tensor alignment = tensor::zero;
label nFarPoints = 0;
for (label pI = 0; pI < 4; ++pI)
if (shapeControlMesh_.dimension() < 3 || shapeControlMesh_.is_infinite(ch))
{
if (ch->vertex(pI)->farPoint())
{
nFarPoints++;
}
}
if (shapeControlMesh_.is_infinite(ch) || nFarPoints == 4)
{
Pout<< "At Infinite vertex" << endl;
if (nFarPoints != 0)
{
for (label pI = 0; pI < 4; ++pI)
{
if (!ch->vertex(pI)->farPoint())
{
alignment = ch->vertex(pI)->alignment();
return alignment;
}
}
}
// cellShapeControlMesh::Vertex_handle nearV =
// shapeControlMesh_.nearest_vertex
// (
// toPoint<cellShapeControlMesh::Point>(pt)
// );
//
// alignment = nearV->alignment();
alignment = tensor::I;
}
else
{
// forAll(bary, pI)
label nFarPoints = 0;
for (label pI = 0; pI < 4; ++pI)
{
if (ch->vertex(pI)->farPoint())
{
nFarPoints++;
}
}
// if (nFarPoints != 0)
// {
// alignment += bary[pI]*ch->vertex(pI)->alignment();
// for (label pI = 0; pI < 4; ++pI)
// {
// if (!ch->vertex(pI)->farPoint())
// {
// alignment = ch->vertex(pI)->alignment();
// }
// }
// }
// else
{
triad tri;
cellShapeControlMesh::Vertex_handle nearV =
shapeControlMesh_.nearest_vertex_in_cell
(
toPoint<cellShapeControlMesh::Point>(pt),
ch
);
for (label pI = 0; pI < 4; ++pI)
{
if (bary[pI] > SMALL)
{
tri += triad(bary[pI]*ch->vertex(pI)->alignment());
}
}
alignment = nearV->alignment();
tri.normalize();
tri.orthogonalize();
// tri = tri.sortxyz();
alignment = tri;
}
// cellShapeControlMesh::Vertex_handle nearV =
// shapeControlMesh_.nearest_vertex_in_cell
// (
// toPoint<cellShapeControlMesh::Point>(pt),
// ch
// );
//
// alignment = nearV->alignment();
}
return alignment;
@ -437,491 +238,81 @@ void Foam::cellShapeControl::cellSizeAndAlignment
alignment = tensor::zero;
size = 0;
label nFarPoints = 0;
for (label pI = 0; pI < 4; ++pI)
if (shapeControlMesh_.dimension() < 3 || shapeControlMesh_.is_infinite(ch))
{
if (ch->vertex(pI)->farPoint())
{
nFarPoints++;
}
}
if (shapeControlMesh_.is_infinite(ch))
{
if (nFarPoints != 0)
{
for (label pI = 0; pI < 4; ++pI)
{
if (!ch->vertex(pI)->farPoint())
{
size = ch->vertex(pI)->targetCellSize();
alignment = ch->vertex(pI)->alignment();
return;
}
}
}
// cellShapeControlMesh::Vertex_handle nearV =
// shapeControlMesh_.nearest_vertex
// (
// toPoint<cellShapeControlMesh::Point>(pt)
// );
//
// size = nearV->targetCellSize();
// alignment = nearV->alignment();
// Find nearest surface
size = sizeAndAlignment_.cellSize(pt);
alignment = tensor::I;
}
else
{
label nFarPoints = 0;
for (label pI = 0; pI < 4; ++pI)
{
if (ch->vertex(pI)->farPoint())
{
nFarPoints++;
}
}
if (nFarPoints != 0)
{
for (label pI = 0; pI < 4; ++pI)
{
if (!ch->vertex(pI)->farPoint())
if (!ch->vertex(pI)->uninitialised())
{
size = ch->vertex(pI)->targetCellSize();
alignment = ch->vertex(pI)->alignment();
return;
}
}
}
else
{
// triad tri;
triad tri;
forAll(bary, pI)
for (label pI = 0; pI < 4; ++pI)
{
size += bary[pI]*ch->vertex(pI)->targetCellSize();
// triad triTmp2 = ch->vertex(pI)->alignment();
// tri += triTmp2*bary[pI];
if (bary[pI] > SMALL)
{
tri += triad(bary[pI]*ch->vertex(pI)->alignment());
}
}
// tri.normalize();
// tri.orthogonalize();
tri.normalize();
tri.orthogonalize();
// tri = tri.sortxyz();
alignment = tri;
// cellShapeControlMesh::Vertex_handle nearV =
// shapeControlMesh_.nearest_vertex
// (
// toPoint<cellShapeControlMesh::Point>(pt)
// );
//
// alignment = tri;
cellShapeControlMesh::Vertex_handle nearV =
shapeControlMesh_.nearest_vertex
(
toPoint<cellShapeControlMesh::Point>(pt)
);
alignment = nearV->alignment();
// alignment = nearV->alignment();
}
}
}
void Foam::cellShapeControl::initialMeshPopulation
(
const autoPtr<backgroundMeshDecomposition>& decomposition
)
{
// Need to pass in the background mesh decomposition so that can test if
// a point to insert is on the processor.
if (Pstream::parRun())
for (label dir = 0; dir < 3; dir++)
{
shapeControlMesh_.insertBoundingPoints(decomposition().procBounds());
}
else
{
shapeControlMesh_.insertBoundingPoints(allGeometry_.bounds());
}
const triad v = alignment;
const PtrList<cellSizeAndAlignmentControl>& controlFunctions =
sizeAndAlignment_.controlFunctions();
forAll(controlFunctions, fI)
{
const cellSizeAndAlignmentControl& controlFunction =
controlFunctions[fI];
Info<< "Inserting points from " << controlFunction.name()
<< " (" << controlFunction.type() << ")" << endl;
pointField pts;
scalarField sizes;
triadField alignments;
controlFunction.initialVertices(pts, sizes, alignments);
List<Vb> vertices(pts.size());
// Clip the minimum size
forAll(vertices, vI)
if (!v.set(dir) || size == 0)
{
vertices[vI] = Vb(pts[vI], Vb::vtInternal);
vertices[vI].targetCellSize() = max(sizes[vI], minimumCellSize_);
vertices[vI].alignment() = alignments[vI];
}
pts.clear();
sizes.clear();
alignments.clear();
label nRejected = 0;
PackedBoolList keepVertex(vertices.size(), true);
if (Pstream::parRun())
{
forAll(vertices, vI)
{
const bool onProc = decomposition().positionOnThisProcessor
(
topoint(vertices[vI].point())
);
if (!onProc)
{
keepVertex[vI] = false;
}
}
}
inplaceSubset(keepVertex, vertices);
const label preInsertedSize = shapeControlMesh_.number_of_vertices();
shapeControlMesh_.rangeInsertWithInfo(vertices.begin(), vertices.end());
Info<< " Inserted "
<< returnReduce
(
shapeControlMesh_.number_of_vertices()
- preInsertedSize, sumOp<label>()
)
<< "/" << vertices.size()
<< endl;
}
}
Foam::label Foam::cellShapeControl::refineMesh
(
const autoPtr<backgroundMeshDecomposition>& decomposition
)
{
const pointField cellCentres(shapeControlMesh_.cellCentres());
Info<< " Created cell centres" << endl;
const PtrList<cellSizeAndAlignmentControl>& controlFunctions =
sizeAndAlignment_.controlFunctions();
DynamicList<Vb> verts(shapeControlMesh_.number_of_vertices());
forAll(cellCentres, cellI)
{
const Foam::point& pt = cellCentres[cellI];
if (!geometryToConformTo_.inside(pt))
{
continue;
}
scalarList bary;
cellShapeControlMesh::Cell_handle ch;
shapeControlMesh_.barycentricCoords(pt, bary, ch);
if (shapeControlMesh_.is_infinite(ch))
{
continue;
}
scalar interpolatedSize = 0;
forAll(bary, pI)
{
interpolatedSize += bary[pI]*ch->vertex(pI)->targetCellSize();
}
label lastPriority = labelMax;
scalar lastCellSize = GREAT;
forAll(controlFunctions, fI)
{
const cellSizeAndAlignmentControl& controlFunction =
controlFunctions[fI];
if (controlFunction.priority() > lastPriority)
{
continue;
}
if (isA<searchableSurfaceControl>(controlFunction))
{
const cellSizeFunction& sizeFunction =
dynamicCast<const searchableSurfaceControl>
(
controlFunction
).sizeFunction();
scalar cellSize = 0;
if (sizeFunction.cellSize(pt, cellSize))
{
if (controlFunction.priority() == lastPriority)
{
if (cellSize < lastCellSize)
{
lastCellSize = cellSize;
}
}
else
{
lastCellSize = cellSize;
}
lastPriority = controlFunction.priority();
}
}
}
if
(
lastCellSize < GREAT
//&& mag(interpolatedSize - lastCellSize)/lastCellSize > 0.2
)
{
if (Pstream::parRun())
{
if (!decomposition().positionOnThisProcessor(pt))
{
continue;
}
}
verts.append
FatalErrorIn
(
Vb
(
toPoint<cellShapeControlMesh::Point>(pt),
Vb::vtInternal
)
);
verts.last().targetCellSize() = lastCellSize;
verts.last().alignment() = triad::unset;
}
}
shapeControlMesh_.insertPoints(verts);
return verts.size();
}
void Foam::cellShapeControl::smoothMesh()
{
label maxSmoothingIterations = readLabel(lookup("maxSmoothingIterations"));
scalar minResidual = 0;
labelListList pointPoints;
autoPtr<mapDistribute> meshDistributor = buildMap
(
shapeControlMesh_,
pointPoints
);
triadField alignments(buildAlignmentField(shapeControlMesh_));
pointField points(buildPointField(shapeControlMesh_));
// Setup the sizes and alignments on each point
triadField fixedAlignments(shapeControlMesh_.vertexCount(), triad::unset);
for
(
CellSizeDelaunay::Finite_vertices_iterator vit =
shapeControlMesh_.finite_vertices_begin();
vit != shapeControlMesh_.finite_vertices_end();
++vit
)
{
if (vit->real())
{
fixedAlignments[vit->index()] = vit->alignment();
}
}
Info<< nl << "Smoothing alignments" << endl;
for (label iter = 0; iter < maxSmoothingIterations; iter++)
{
Info<< "Iteration " << iter;
meshDistributor().distribute(points);
meshDistributor().distribute(alignments);
scalar residual = 0;
triadField triadAv(alignments.size(), triad::unset);
forAll(pointPoints, pI)
{
const labelList& pPoints = pointPoints[pI];
if (pPoints.empty())
{
continue;
}
triad& newTriad = triadAv[pI];
forAll(pPoints, adjPointI)
{
const label adjPointIndex = pPoints[adjPointI];
scalar dist = mag(points[pI] - points[adjPointIndex]);
triad tmpTriad = alignments[adjPointIndex];
for (direction dir = 0; dir < 3; dir++)
{
if (tmpTriad.set(dir))
{
tmpTriad[dir] *= (1.0/dist);
}
}
newTriad += tmpTriad;
}
newTriad.normalize();
newTriad.orthogonalize();
newTriad = newTriad.sortxyz();
// Enforce the boundary conditions
const triad& fixedAlignment = fixedAlignments[pI];
label nFixed = 0;
forAll(fixedAlignment, dirI)
{
if (fixedAlignment.set(dirI))
{
nFixed++;
}
}
if (nFixed == 1)
{
forAll(fixedAlignment, dirI)
{
if (fixedAlignment.set(dirI))
{
newTriad.align(fixedAlignment[dirI]);
}
}
}
else if (nFixed == 2)
{
forAll(fixedAlignment, dirI)
{
if (fixedAlignment.set(dirI))
{
newTriad[dirI] = fixedAlignment[dirI];
}
else
{
newTriad[dirI] = triad::unset[dirI];
}
}
newTriad.orthogonalize();
}
else if (nFixed == 3)
{
forAll(fixedAlignment, dirI)
{
if (fixedAlignment.set(dirI))
{
newTriad[dirI] = fixedAlignment[dirI];
}
}
}
const triad& oldTriad = alignments[pI];
for (direction dir = 0; dir < 3; ++dir)
{
if
(
newTriad.set(dir)
&& oldTriad.set(dir)
//&& !fixedAlignment.set(dir)
)
{
scalar dotProd = (oldTriad[dir] & newTriad[dir]);
scalar diff = mag(dotProd) - 1.0;
residual += mag(diff);
}
}
// if (iter == 198 || iter == 199)
// {
// Info<< "Triad" << nl
// << " Fixed (" << nFixed << ") = " << fixedAlignment
// << nl
// << " Old = " << oldTriad << nl
// << " Pre-Align= " << triadAv[pI] << nl
// << " New = " << newTriad << nl
// << " Residual = " << residual << endl;
// }
}
forAll(alignments, pI)
{
alignments[pI] = triadAv[pI].sortxyz();
}
reduce(residual, sumOp<scalar>());
Info<< ", Residual = " << residual << endl;
if (iter > 0 && residual <= minResidual)
{
break;
}
}
meshDistributor().distribute(alignments);
for
(
CellSizeDelaunay::Finite_vertices_iterator vit =
shapeControlMesh_.finite_vertices_begin();
vit != shapeControlMesh_.finite_vertices_end();
++vit
)
{
if (vit->real())
{
vit->alignment() = alignments[vit->index()];
}
}
labelList referredPoints;
autoPtr<mapDistribute> referredDistributor = buildReferredMap
(
shapeControlMesh_,
referredPoints
);
alignments.setSize(shapeControlMesh_.vertexCount());
referredDistributor().distribute(alignments);
label referredI = 0;
for
(
CellSizeDelaunay::Finite_vertices_iterator vit =
shapeControlMesh_.finite_vertices_begin();
vit != shapeControlMesh_.finite_vertices_end();
++vit
)
{
if (vit->referred())
{
vit->alignment() = alignments[referredPoints[referredI++]];
"Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()"
) << "Point has bad alignment! "
<< pt << " " << size << " " << alignment << nl
<< "Bary Coords = " << bary << nl
<< ch->vertex(0)->info() << nl
<< ch->vertex(1)->info() << nl
<< ch->vertex(2)->info() << nl
<< ch->vertex(3)->info()
<< abort(FatalError);
}
}
}

View File

@ -84,33 +84,6 @@ class cellShapeControl
// Private Member Functions
template <class Triangulation, class Type>
tmp<Field<Type> > filterFarPoints
(
const Triangulation& mesh,
const Field<Type>& field
);
template <class Triangulation>
autoPtr<mapDistribute> buildMap
(
const Triangulation& mesh,
labelListList& pointPoints
);
template <class Triangulation>
autoPtr<mapDistribute> buildReferredMap
(
const Triangulation& mesh,
labelList& indices
);
template <class Triangulation>
tmp<triadField> buildAlignmentField(const Triangulation& mesh);
template <class Triangulation>
tmp<pointField> buildPointField(const Triangulation& mesh);
//- Disallow default bitwise copy construct
cellShapeControl(const cellShapeControl&);
@ -155,6 +128,8 @@ public:
inline const cellSizeAndAlignmentControls& sizeAndAlignment() const;
inline const scalar& minimumCellSize() const;
// Query
@ -172,29 +147,6 @@ public:
scalar& size,
tensor& alignment
) const;
// Edit
void initialMeshPopulation
(
const autoPtr<backgroundMeshDecomposition>& decomposition
);
label refineMesh
(
const autoPtr<backgroundMeshDecomposition>& decomposition
);
void smoothMesh();
//- Add a control point with a specified size and alignment
// virtual void addControlPoint
// (
// const point& pt,
// const scalar& size,
// const tensor& alignment
// );
};

View File

@ -52,4 +52,17 @@ Foam::cellShapeControl::aspectRatio() const
}
inline const Foam::cellSizeAndAlignmentControls&
Foam::cellShapeControl::sizeAndAlignment() const
{
return sizeAndAlignment_;
}
inline const Foam::scalar& Foam::cellShapeControl::minimumCellSize() const
{
return minimumCellSize_;
}
// ************************************************************************* //

View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "cellShapeControlMesh.H"
#include "cellSizeAndAlignmentControls.H"
#include "pointIOField.H"
#include "scalarIOField.H"
#include "tensorIOField.H"
@ -36,9 +37,7 @@ License
namespace Foam
{
defineTypeNameAndDebug(cellShapeControlMesh, 0);
}
@ -537,12 +536,9 @@ void Foam::cellShapeControlMesh::barycentricCoords
) const
{
// Use the previous cell handle as a hint on where to start searching
ch = locate
(
Point(pt.x(), pt.y(), pt.z())
);
ch = locate(toPoint<Point>(pt));
if (!is_infinite(ch))
if (dimension() > 2 && !is_infinite(ch))
{
oldCellHandle_ = ch;

View File

@ -53,6 +53,8 @@ SourceFiles
namespace Foam
{
class cellSizeAndAlignmentControls;
/*---------------------------------------------------------------------------*\
Class cellShapeControlMesh Declaration
\*---------------------------------------------------------------------------*/

View File

@ -48,7 +48,7 @@ Foam::cellSizeAndAlignmentControl::cellSizeAndAlignmentControl
const Time& runTime,
const word& name,
const dictionary& controlFunctionDict,
const conformationSurfaces& allGeometry
const conformationSurfaces& geometryToConformTo
)
:
runTime_(runTime),
@ -64,7 +64,7 @@ Foam::cellSizeAndAlignmentControl::New
const Time& runTime,
const word& name,
const dictionary& controlFunctionDict,
const conformationSurfaces& allGeometry
const conformationSurfaces& geometryToConformTo
)
{
word cellSizeAndAlignmentControlTypeName
@ -72,7 +72,7 @@ Foam::cellSizeAndAlignmentControl::New
controlFunctionDict.lookup("type")
);
Info<< nl << "Selecting cellSizeAndAlignmentControl "
Info<< " Selecting cellSizeAndAlignmentControl "
<< cellSizeAndAlignmentControlTypeName << endl;
dictionaryConstructorTable::iterator cstrIter =
@ -101,7 +101,7 @@ Foam::cellSizeAndAlignmentControl::New
runTime,
name,
controlFunctionDict,
allGeometry
geometryToConformTo
)
);
}

View File

@ -91,9 +91,9 @@ public:
const Time& runTime,
const word& name,
const dictionary& controlFunctionDict,
const conformationSurfaces& allGeometry
const conformationSurfaces& geometryToConformTo
),
(runTime, name, controlFunctionDict, allGeometry)
(runTime, name, controlFunctionDict, geometryToConformTo)
);
@ -106,7 +106,7 @@ public:
const Time& runTime,
const word& name,
const dictionary& controlFunctionDict,
const conformationSurfaces& allGeometry
const conformationSurfaces& geometryToConformTo
);
@ -118,7 +118,7 @@ public:
const Time& runTime,
const word& name,
const dictionary& controlFunctionDict,
const conformationSurfaces& allGeometry
const conformationSurfaces& geometryToConformTo
);

View File

@ -30,9 +30,7 @@ License
namespace Foam
{
defineTypeNameAndDebug(cellSizeAndAlignmentControls, 0);
}
@ -116,12 +114,12 @@ Foam::cellSizeAndAlignmentControls::cellSizeAndAlignmentControls
(
const Time& runTime,
const dictionary& shapeControlDict,
const conformationSurfaces& allGeometry,
const conformationSurfaces& geometryToConformTo,
const scalar defaultCellSize
)
:
shapeControlDict_(shapeControlDict),
allGeometry_(allGeometry),
geometryToConformTo_(geometryToConformTo),
controlFunctions_(shapeControlDict_.size()),
defaultCellSize_(defaultCellSize)
{
@ -136,6 +134,8 @@ Foam::cellSizeAndAlignmentControls::cellSizeAndAlignmentControls
shapeControlDict_.subDict(shapeControlEntryName)
);
Info<< nl << "Shape Control : " << shapeControlEntryName << endl;
controlFunctions_.set
(
functionI,
@ -144,7 +144,7 @@ Foam::cellSizeAndAlignmentControls::cellSizeAndAlignmentControls
runTime,
shapeControlEntryName,
controlFunctionDict,
allGeometry
geometryToConformTo
)
);

View File

@ -52,7 +52,7 @@ class cellSizeAndAlignmentControls
const dictionary& shapeControlDict_;
const conformationSurfaces& allGeometry_;
const conformationSurfaces& geometryToConformTo_;
PtrList<cellSizeAndAlignmentControl> controlFunctions_;
@ -83,7 +83,7 @@ public:
(
const Time& runTime,
const dictionary& shapeControlDict,
const conformationSurfaces& allGeometry,
const conformationSurfaces& geometryToConformTo,
const scalar defaultCellSize
);
@ -96,11 +96,17 @@ public:
// Access
const PtrList<cellSizeAndAlignmentControl>& controlFunctions() const
inline const PtrList<cellSizeAndAlignmentControl>&
controlFunctions() const
{
return controlFunctions_;
}
inline const conformationSurfaces& geometryToConformTo() const
{
return geometryToConformTo_;
}
// Query

View File

@ -58,7 +58,7 @@ Foam::fileControl::fileControl
const Time& runTime,
const word& name,
const dictionary& controlFunctionDict,
const conformationSurfaces& allGeometry
const conformationSurfaces& geometryToConformTo
)
:
cellSizeAndAlignmentControl
@ -66,7 +66,7 @@ Foam::fileControl::fileControl
runTime,
name,
controlFunctionDict,
allGeometry
geometryToConformTo
),
pointsFile_(controlFunctionDict.lookup("pointsFile")),
sizesFile_(controlFunctionDict.lookup("sizesFile")),

View File

@ -85,7 +85,7 @@ public:
const Time& runTime,
const word& name,
const dictionary& controlFunctionDict,
const conformationSurfaces& allGeometry
const conformationSurfaces& geometryToConformTo
);
//- Destructor

View File

@ -158,7 +158,7 @@ Foam::searchableSurfaceControl::searchableSurfaceControl
const Time& runTime,
const word& name,
const dictionary& controlFunctionDict,
const conformationSurfaces& allGeometry
const conformationSurfaces& geometryToConformTo
)
:
cellSizeAndAlignmentControl
@ -166,11 +166,11 @@ Foam::searchableSurfaceControl::searchableSurfaceControl
runTime,
name,
controlFunctionDict,
allGeometry
geometryToConformTo
),
surfaceName_(controlFunctionDict.lookupOrDefault<word>("surface", name)),
searchableSurface_(allGeometry.geometry()[surfaceName_]),
allGeometry_(allGeometry),
searchableSurface_(geometryToConformTo.geometry()[surfaceName_]),
geometryToConformTo_(geometryToConformTo),
cellSizeFunction_
(
cellSizeFunction::New(controlFunctionDict, searchableSurface_)
@ -579,7 +579,7 @@ void Foam::searchableSurfaceControl::initialVertices
pointIndexHit info;
label infoFeature;
allGeometry_.findFeaturePointNearest
geometryToConformTo_.findFeaturePointNearest
(
pts[pI],
nearFeatDistSqrCoeff,
@ -592,7 +592,7 @@ void Foam::searchableSurfaceControl::initialVertices
if (info.hit())
{
const extendedFeatureEdgeMesh& features =
allGeometry_.features()[infoFeature];
geometryToConformTo_.features()[infoFeature];
vectorField norms = features.featurePointNormals(info.index());
@ -608,7 +608,7 @@ void Foam::searchableSurfaceControl::initialVertices
}
else
{
allGeometry_.findEdgeNearest
geometryToConformTo_.findEdgeNearest
(
pts[pI],
nearFeatDistSqrCoeff,
@ -619,7 +619,7 @@ void Foam::searchableSurfaceControl::initialVertices
if (info.hit())
{
const extendedFeatureEdgeMesh& features =
allGeometry_.features()[infoFeature];
geometryToConformTo_.features()[infoFeature];
vectorField norms = features.edgeNormals(info.index());

View File

@ -59,7 +59,7 @@ class searchableSurfaceControl
//- Reference to the searchableSurface object holding the geometry data
const searchableSurface& searchableSurface_;
const conformationSurfaces& allGeometry_;
const conformationSurfaces& geometryToConformTo_;
autoPtr<cellSizeFunction> cellSizeFunction_;
@ -106,7 +106,7 @@ public:
const Time& runTime,
const word& name,
const dictionary& controlFunctionDict,
const conformationSurfaces& allGeometry
const conformationSurfaces& geometryToConformTo
);
//- Destructor

View File

@ -0,0 +1,683 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "controlMeshRefinement.H"
#include "cellSizeAndAlignmentControl.H"
#include "OFstream.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(controlMeshRefinement, 0);
}
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::controlMeshRefinement::detectEdge
(
point a,
point b,
DynamicList<Vb>& pointsFound
) const
{
const Foam::point midpoint = 0.5*(a + b);
const scalar h = mag(a - b);
scalar cellSizeA = sizeControls_.cellSize(a);
scalar cellSizeMid = sizeControls_.cellSize(midpoint);
scalar cellSizeB = sizeControls_.cellSize(b);
scalar firstDerivative = (cellSizeA - cellSizeB)/h;
scalar secondDerivative =
(
cellSizeB
- 2*cellSizeMid
+ cellSizeA
)
/magSqr(h/2.0);
// Info<< a << " " << midpoint << " " << b << endl;
// Info<< " length = " << h << endl;
// Info<< "size(a) = " << cellSizeA << endl;
// Info<< "size(m) = " << cellSizeMid << endl;
// Info<< "size(b) = " << cellSizeB << endl;
// Info<< " f' = " << firstDerivative << endl;
// Info<< " f'' = " << secondDerivative << endl;
if (mag(secondDerivative) < 1e-1)
{
return false;
}
if (mag(a - b) < 1e-1*cellSizeMid)
{
if (geometryToConformTo_.outside(midpoint))
{
return false;
}
else
{
// Info<< " Point Added" << endl;
pointsFound.append
(
Vb
(
toPoint<cellShapeControlMesh::Point>(midpoint),
Vb::vtInternal
)
);
pointsFound.last().targetCellSize() =
sizeControls_.cellSize(midpoint);
pointsFound.last().alignment() = triad::unset;
return true;
}
}
detectEdge(a, midpoint, pointsFound);
detectEdge(midpoint, b, pointsFound);
}
Foam::DynamicList<Vb>
Foam::controlMeshRefinement::findDiscontinuities(const linePointRef& l) const
{
DynamicList<Vb> verts(2);
// Divide up the line into a reasonable number of intervals
// Call detectEdge on each interval
detectEdge(l.start(), l.end(), verts);
if (verts.size())
{
Info<< l << " has " << verts.size() << " discontinuities" << endl;
OFstream str("plotSizes_" + name(verts.size()));
forAll(verts, vI)
{
scalar x =
mag(topoint(verts[vI].point()) - l.start())
/mag(l.end() - l.start());
str << x << " " << verts[vI].targetCellSize() << endl;
}
// forAll(verts, vI)
// {
// Info<< verts[vI].info();
// }
}
return verts;
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::controlMeshRefinement::controlMeshRefinement
(
cellShapeControl& shapeController
)
:
shapeController_(shapeController),
mesh_(shapeController.shapeControlMesh()),
sizeControls_(shapeController.sizeAndAlignment()),
geometryToConformTo_(sizeControls_.geometryToConformTo())
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::controlMeshRefinement::~controlMeshRefinement()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::controlMeshRefinement::initialMeshPopulation
(
const autoPtr<backgroundMeshDecomposition>& decomposition
)
{
autoPtr<boundBox> overallBoundBox;
// Need to pass in the background mesh decomposition so that can test if
// a point to insert is on the processor.
if (Pstream::parRun())
{
overallBoundBox.set(new boundBox(decomposition().procBounds()));
}
else
{
overallBoundBox.set
(
new boundBox(geometryToConformTo_.geometry().bounds())
);
}
mesh_.insertBoundingPoints
(
overallBoundBox(),
sizeControls_
);
const PtrList<cellSizeAndAlignmentControl>& controlFunctions =
sizeControls_.controlFunctions();
forAll(controlFunctions, fI)
{
const cellSizeAndAlignmentControl& controlFunction =
controlFunctions[fI];
Info<< "Inserting points from " << controlFunction.name()
<< " (" << controlFunction.type() << ")" << endl;
pointField pts;
scalarField sizes;
triadField alignments;
controlFunction.initialVertices(pts, sizes, alignments);
List<Vb> vertices(pts.size());
// Clip the minimum size
forAll(vertices, vI)
{
vertices[vI] = Vb(pts[vI], Vb::vtInternalNearBoundary);
vertices[vI].targetCellSize() = max
(
sizes[vI],
shapeController_.minimumCellSize()
);
vertices[vI].alignment() = alignments[vI];
}
pts.clear();
sizes.clear();
alignments.clear();
label nRejected = 0;
PackedBoolList keepVertex(vertices.size(), true);
if (Pstream::parRun())
{
forAll(vertices, vI)
{
const bool onProc = decomposition().positionOnThisProcessor
(
topoint(vertices[vI].point())
);
if (!onProc)
{
keepVertex[vI] = false;
}
}
}
inplaceSubset(keepVertex, vertices);
const label preInsertedSize = mesh_.number_of_vertices();
forAll(vertices, vI)
{
bool insertPoint = false;
pointFromPoint pt(topoint(vertices[vI].point()));
if
(
mesh_.dimension() < 3
|| mesh_.is_infinite
(
mesh_.locate(vertices[vI].point())
)
)
{
insertPoint = true;
}
const scalar interpolatedCellSize = shapeController_.cellSize(pt);
const triad interpolatedAlignment =
shapeController_.cellAlignment(pt);
const scalar calculatedCellSize = vertices[vI].targetCellSize();
const triad calculatedAlignment = vertices[vI].alignment();
if (debug)
{
Info<< "Point = " << pt << nl
<< " Size(interp) = " << interpolatedCellSize << nl
<< " Size(calc) = " << calculatedCellSize << nl
<< " Align(interp) = " << interpolatedAlignment << nl
<< " Align(calc) = " << calculatedAlignment << nl
<< endl;
}
const scalar sizeDiff =
mag(interpolatedCellSize - calculatedCellSize);
const scalar alignmentDiff =
diff(interpolatedAlignment, calculatedAlignment);
if (debug)
{
Info<< " size difference = " << sizeDiff
<< ", alignment difference = " << alignmentDiff << endl;
}
// @todo Also need to base it on the alignments
if
(
sizeDiff/interpolatedCellSize > 0.1
|| alignmentDiff > 0.15
)
{
insertPoint = true;
}
if (insertPoint)
{
mesh_.insert
(
pt,
calculatedCellSize,
vertices[vI].alignment(),
Vb::vtInternalNearBoundary
);
}
}
//mesh_.rangeInsertWithInfo(vertices.begin(), vertices.end());
Info<< " Inserted "
<< returnReduce
(
label(mesh_.number_of_vertices()) - preInsertedSize,
sumOp<label>()
)
<< "/" << returnReduce(vertices.size(), sumOp<label>())
<< endl;
}
}
Foam::label Foam::controlMeshRefinement::refineMesh
(
const autoPtr<backgroundMeshDecomposition>& decomposition
)
{
Info<< "Iterate over cell size mesh edges" << endl;
DynamicList<Vb> verts(mesh_.number_of_vertices());
for
(
CellSizeDelaunay::Finite_edges_iterator eit =
mesh_.finite_edges_begin();
eit != mesh_.finite_edges_end();
++eit
)
{
CellSizeDelaunay::Cell_handle c = eit->first;
CellSizeDelaunay::Vertex_handle vA = c->vertex(eit->second);
CellSizeDelaunay::Vertex_handle vB = c->vertex(eit->third);
if
(
mesh_.is_infinite(vA)
|| mesh_.is_infinite(vB)
|| (vA->referred() && vB->referred())
|| (vA->referred() && vA->procIndex() > vB->procIndex())
|| (vB->referred() && vB->procIndex() > vA->procIndex())
)
{
continue;
}
pointFromPoint ptA(topoint(vA->point()));
pointFromPoint ptB(topoint(vB->point()));
linePointRef l(ptA, ptB);
// Info<< "Edge " << l << endl;
//
// Info<< vA->info() << vB->info();
const Foam::point midPoint = l.centre();
const scalar length = l.mag();
scalar evaluatedSize = sizeControls_.cellSize(midPoint);
scalar interpolatedSize =
(vA->targetCellSize() + vB->targetCellSize())/2;
const scalar diff = mag(interpolatedSize - evaluatedSize);
// Info<< " Evaluated size = " << evaluatedSize << nl
// << " Interpolated size = " << interpolatedSize << nl
// << " vA cell size = " << vA->targetCellSize() << nl
// << " vB cell size = " << vB->targetCellSize() << endl;
// verts.append(findDiscontinuities(l));
if (diff/interpolatedSize > 0.1)
{
// Create a new point
// Walk along the edge, placing points where the gradient in cell
// size changes
if (vA->targetCellSize() <= vB->targetCellSize())
{
// Walk from vA
const label nPoints = length/vA->targetCellSize();
if (nPoints == 0)
{
continue;
}
scalar prevSize = vA->targetCellSize();
for (label pI = 1; pI < nPoints; ++pI)
{
const Foam::point testPt =
topoint(vA->point()) + pI*l.vec()/nPoints;
// if (geometryToConformTo_.outside(testPt))
// {
// continue;
// }
scalar testPtSize = sizeControls_.cellSize(testPt);
if (mag(testPtSize - prevSize) < 1e-3*testPtSize)
{
// Info<< "Adding " << testPt << " " << testPtSize
// << endl;
verts.append
(
Vb
(
toPoint<cellShapeControlMesh::Point>(testPt),
Vb::vtInternal
)
);
verts.last().targetCellSize() = testPtSize;
verts.last().alignment() = triad::unset;
break;
}
prevSize = testPtSize;
}
}
else
{
// Walk from vB
const label nPoints = length/vB->targetCellSize();
if (nPoints == 0)
{
continue;
}
scalar prevSize = vA->targetCellSize();
for (label pI = 1; pI < nPoints; ++pI)
{
const Foam::point testPt =
topoint(vB->point()) - pI*l.vec()/nPoints;
// if (geometryToConformTo_.outside(testPt))
// {
// continue;
// }
scalar testPtSize = sizeControls_.cellSize(testPt);
if (mag(testPtSize - prevSize) < 1e-3*testPtSize)
//if (testPtSize > prevSize)// || testPtSize < prevSize)
{
// Info<< "Adding " << testPt << " " << testPtSize
// << endl;
verts.append
(
Vb
(
toPoint<cellShapeControlMesh::Point>(testPt),
Vb::vtInternal
)
);
verts.last().targetCellSize() = testPtSize;
verts.last().alignment() = triad::unset;
break;
}
prevSize = testPtSize;
}
}
}
}
// const pointField cellCentres(mesh_.cellCentres());
//
// Info<< " Created cell centres" << endl;
//
// const PtrList<cellSizeAndAlignmentControl>& controlFunctions =
// sizeControls_.controlFunctions();
//
// DynamicList<Vb> verts(mesh_.number_of_vertices());
//
// forAll(cellCentres, cellI)
// {
// Foam::point pt = cellCentres[cellI];
//
// if (!geometryToConformTo_.inside(pt))
// {
// continue;
// }
//
// scalarList bary;
// cellShapeControlMesh::Cell_handle ch;
//
// mesh_.barycentricCoords(pt, bary, ch);
//
// if (mesh_.is_infinite(ch))
// {
// continue;
// }
//
// scalar interpolatedSize = 0;
// forAll(bary, pI)
// {
// interpolatedSize += bary[pI]*ch->vertex(pI)->targetCellSize();
// }
//
// label lastPriority = labelMin;
// scalar lastCellSize = GREAT;
// forAll(controlFunctions, fI)
// {
// const cellSizeAndAlignmentControl& controlFunction =
// controlFunctions[fI];
//
// if (controlFunction.priority() < lastPriority)
// {
// continue;
// }
//
// if (isA<searchableSurfaceControl>(controlFunction))
// {
// const cellSizeFunction& sizeFunction =
// dynamicCast<const searchableSurfaceControl>
// (
// controlFunction
// ).sizeFunction();
//
// scalar cellSize = 0;
// if (sizeFunction.cellSize(pt, cellSize))
// {
// if (controlFunction.priority() == lastPriority)
// {
// if (cellSize < lastCellSize)
// {
// lastCellSize = cellSize;
// }
// }
// else
// {
// lastCellSize = cellSize;
// }
//
// lastPriority = controlFunction.priority();
// }
// }
// }
//
// const scalar diff = mag(interpolatedSize - lastCellSize);
//
// if (diff/interpolatedSize > 0.1)
// {
// bool pointFound = false;
//
// // Travel along lines to each tet point
// for (label vI = 0; vI < 4; ++vI)
// {
// if (pointFound)
// {
// break;
// }
//
// pointFromPoint tetVertex = topoint(ch->vertex(vI)->point());
// scalar tetVertexSize = sizeControls_.cellSize(tetVertex);
//
// if (tetVertexSize < interpolatedSize)
// {
// linePointRef l(tetVertex, pt);
//
// label nTestPoints = l.mag()/tetVertexSize;
//
// for (label i = 1; i < nTestPoints; ++i)
// {
// const Foam::point testPt
// (
// tetVertex
// + l.vec()*i/nTestPoints
// );
//
// scalar testPtSize = sizeControls_.cellSize(testPt);
//
// if (testPtSize > tetVertexSize)
// {
// pt = testPt;
// lastCellSize = testPtSize;
// pointFound = true;
// break;
// }
// }
// }
// else
// {
// linePointRef l(pt, tetVertex);
//
// label nTestPoints = l.mag()/tetVertexSize;
//
// for (label i = 1; i < nTestPoints; ++i)
// {
// const Foam::point testPt
// (
// tetVertex
// + l.vec()*i/nTestPoints
// );
//
// scalar testPtSize = sizeControls_.cellSize(testPt);
//
// if (testPtSize < interpolatedSize)
// {
// pt = testPt;
// lastCellSize = testPtSize;
// pointFound = true;
// break;
// }
// }
// }
// }
// }
//
// // Decide whether to insert or not
// if (lastCellSize < GREAT)
// {
// if (!geometryToConformTo_.inside(pt))
// {
// continue;
// }
//
// if (Pstream::parRun())
// {
// if (!decomposition().positionOnThisProcessor(pt))
// {
// continue;
// }
// }
//
// verts.append
// (
// Vb
// (
// toPoint<cellShapeControlMesh::Point>(pt),
// Vb::vtInternal
// )
// );
// verts.last().targetCellSize() = lastCellSize;
// verts.last().alignment() = triad::unset;
// }
// }
mesh_.insertPoints(verts);
return verts.size();
}
// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -0,0 +1,133 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::controlMeshRefinement
Description
SourceFiles
controlMeshRefinementI.H
controlMeshRefinement.C
controlMeshRefinementIO.C
\*---------------------------------------------------------------------------*/
#ifndef controlMeshRefinement_H
#define controlMeshRefinement_H
#include "cellShapeControl.H"
#include "cellShapeControlMesh.H"
#include "cellSizeAndAlignmentControls.H"
#include "conformationSurfaces.H"
#include "backgroundMeshDecomposition.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class controlMeshRefinement Declaration
\*---------------------------------------------------------------------------*/
class controlMeshRefinement
{
// Private data
const cellShapeControl& shapeController_;
cellShapeControlMesh& mesh_;
const cellSizeAndAlignmentControls& sizeControls_;
const conformationSurfaces& geometryToConformTo_;
// Private Member Functions
bool detectEdge
(
point a,
point b,
DynamicList<Vb>& pointsFound
) const;
DynamicList<Vb> findDiscontinuities(const linePointRef& l) const;
//- Disallow default bitwise copy construct
controlMeshRefinement(const controlMeshRefinement&);
//- Disallow default bitwise assignment
void operator=(const controlMeshRefinement&);
public:
//- Runtime type information
ClassName("controlMeshRefinement");
// Constructors
//- Construct null
controlMeshRefinement(cellShapeControl& shapeController);
//- Destructor
~controlMeshRefinement();
// Member Functions
// Access
// Check
// Edit
void initialMeshPopulation
(
const autoPtr<backgroundMeshDecomposition>& decomposition
);
label refineMesh
(
const autoPtr<backgroundMeshDecomposition>& decomposition
);
// Write
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,488 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "smoothAlignmentSolver.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template <class Triangulation, class Type>
Foam::tmp<Foam::Field<Type> > Foam::smoothAlignmentSolver::filterFarPoints
(
const Triangulation& mesh,
const Field<Type>& field
)
{
tmp<Field<Type> > tNewField(new Field<Type>(field.size()));
Field<Type>& newField = tNewField();
label added = 0;
label count = 0;
for
(
typename Triangulation::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
vit != mesh.finite_vertices_end();
++vit
)
{
if (vit->real())
{
newField[added++] = field[count];
}
count++;
}
newField.resize(added);
return tNewField;
}
template <class Triangulation>
Foam::autoPtr<Foam::mapDistribute> Foam::smoothAlignmentSolver::buildReferredMap
(
const Triangulation& mesh,
labelList& indices
)
{
globalIndex globalIndexing(mesh.vertexCount());
DynamicList<label> dynIndices(mesh.vertexCount()/10);
for
(
typename Triangulation::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
vit != mesh.finite_vertices_end();
++vit
)
{
if (vit->referred())
{
dynIndices.append
(
globalIndexing.toGlobal(vit->procIndex(), vit->index())
);
}
}
indices.transfer(dynIndices);
List<Map<label> > compactMap;
return autoPtr<mapDistribute>
(
new mapDistribute
(
globalIndexing,
indices,
compactMap
)
);
}
template <class Triangulation>
Foam::autoPtr<Foam::mapDistribute> Foam::smoothAlignmentSolver::buildMap
(
const Triangulation& mesh,
labelListList& pointPoints
)
{
pointPoints.setSize(mesh.vertexCount());
globalIndex globalIndexing(mesh.vertexCount());
for
(
typename Triangulation::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
vit != mesh.finite_vertices_end();
++vit
)
{
if (!vit->real())
{
continue;
}
std::list<typename Triangulation::Vertex_handle> adjVerts;
mesh.finite_adjacent_vertices(vit, std::back_inserter(adjVerts));
DynamicList<label> indices(adjVerts.size());
for
(
typename std::list<typename Triangulation::Vertex_handle>::
const_iterator adjVertI = adjVerts.begin();
adjVertI != adjVerts.end();
++adjVertI
)
{
typename Triangulation::Vertex_handle vh = *adjVertI;
if (!vh->farPoint())
{
indices.append
(
globalIndexing.toGlobal(vh->procIndex(), vh->index())
);
}
}
pointPoints[vit->index()].transfer(indices);
}
List<Map<label> > compactMap;
return autoPtr<mapDistribute>
(
new mapDistribute
(
globalIndexing,
pointPoints,
compactMap
)
);
}
template <class Triangulation>
Foam::tmp<Foam::triadField> Foam::smoothAlignmentSolver::buildAlignmentField
(
const Triangulation& mesh
)
{
tmp<triadField> tAlignments
(
new triadField(mesh.vertexCount(), triad::unset)
);
triadField& alignments = tAlignments();
for
(
typename Triangulation::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
vit != mesh.finite_vertices_end();
++vit
)
{
if (!vit->real())
{
continue;
}
alignments[vit->index()] = vit->alignment();
}
return tAlignments;
}
template <class Triangulation>
Foam::tmp<Foam::pointField> Foam::smoothAlignmentSolver::buildPointField
(
const Triangulation& mesh
)
{
tmp<pointField> tPoints
(
new pointField(mesh.vertexCount(), point(GREAT, GREAT, GREAT))
);
pointField& points = tPoints();
for
(
typename Triangulation::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
vit != mesh.finite_vertices_end();
++vit
)
{
if (!vit->real())
{
continue;
}
points[vit->index()] = topoint(vit->point());
}
return tPoints;
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::smoothAlignmentSolver::smoothAlignmentSolver(cellShapeControlMesh& mesh)
:
mesh_(mesh)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::smoothAlignmentSolver::~smoothAlignmentSolver()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::smoothAlignmentSolver::smoothAlignments
(
const label maxSmoothingIterations
)
{
scalar minResidual = 0;
labelListList pointPoints;
autoPtr<mapDistribute> meshDistributor = buildMap
(
mesh_,
pointPoints
);
triadField alignments(buildAlignmentField(mesh_));
pointField points(buildPointField(mesh_));
// Setup the sizes and alignments on each point
triadField fixedAlignments(mesh_.vertexCount(), triad::unset);
for
(
CellSizeDelaunay::Finite_vertices_iterator vit =
mesh_.finite_vertices_begin();
vit != mesh_.finite_vertices_end();
++vit
)
{
if (vit->real())
{
fixedAlignments[vit->index()] = vit->alignment();
}
}
Info<< nl << "Smoothing alignments" << endl;
for (label iter = 0; iter < maxSmoothingIterations; iter++)
{
Info<< "Iteration " << iter;
meshDistributor().distribute(points);
meshDistributor().distribute(alignments);
scalar residual = 0;
triadField triadAv(alignments.size(), triad::unset);
forAll(pointPoints, pI)
{
const labelList& pPoints = pointPoints[pI];
if (pPoints.empty())
{
continue;
}
triad& newTriad = triadAv[pI];
forAll(pPoints, adjPointI)
{
const label adjPointIndex = pPoints[adjPointI];
scalar dist = mag(points[pI] - points[adjPointIndex]);
triad tmpTriad = alignments[adjPointIndex];
for (direction dir = 0; dir < 3; dir++)
{
if (tmpTriad.set(dir))
{
tmpTriad[dir] *= 1.0/(dist + SMALL);
}
}
newTriad += tmpTriad;
}
newTriad.normalize();
newTriad.orthogonalize();
// newTriad = newTriad.sortxyz();
// Enforce the boundary conditions
const triad& fixedAlignment = fixedAlignments[pI];
label nFixed = 0;
forAll(fixedAlignment, dirI)
{
if (fixedAlignment.set(dirI))
{
nFixed++;
}
}
if (nFixed == 1)
{
forAll(fixedAlignment, dirI)
{
if (fixedAlignment.set(dirI))
{
newTriad.align(fixedAlignment[dirI]);
}
}
}
else if (nFixed == 2)
{
forAll(fixedAlignment, dirI)
{
if (fixedAlignment.set(dirI))
{
newTriad[dirI] = fixedAlignment[dirI];
}
else
{
newTriad[dirI] = triad::unset[dirI];
}
}
newTriad.orthogonalize();
}
else if (nFixed == 3)
{
forAll(fixedAlignment, dirI)
{
if (fixedAlignment.set(dirI))
{
newTriad[dirI] = fixedAlignment[dirI];
}
}
}
const triad& oldTriad = alignments[pI];
for (direction dir = 0; dir < 3; ++dir)
{
if
(
newTriad.set(dir)
&& oldTriad.set(dir)
&& !fixedAlignment.set(dir)
)
{
residual += diff(oldTriad, newTriad);
}
}
// if (iter == 198 || iter == 199)
// {
// Info<< "Triad" << nl
// << " Fixed (" << nFixed << ") = " << fixedAlignment
// << nl
// << " Old = " << oldTriad << nl
// << " Pre-Align= " << triadAv[pI] << nl
// << " New = " << newTriad << nl
// << " Residual = " << residual << endl;
// }
}
forAll(alignments, pI)
{
alignments[pI] = triadAv[pI].sortxyz();
}
reduce(residual, sumOp<scalar>());
Info<< ", Residual = "
<< residual
/returnReduce(points.size(), sumOp<label>())
<< endl;
if (iter > 0 && residual <= minResidual)
{
break;
}
}
meshDistributor().distribute(alignments);
for
(
CellSizeDelaunay::Finite_vertices_iterator vit =
mesh_.finite_vertices_begin();
vit != mesh_.finite_vertices_end();
++vit
)
{
if (vit->real())
{
vit->alignment() = alignments[vit->index()];
}
}
labelList referredPoints;
autoPtr<mapDistribute> referredDistributor = buildReferredMap
(
mesh_,
referredPoints
);
alignments.setSize(mesh_.vertexCount());
referredDistributor().distribute(alignments);
label referredI = 0;
for
(
CellSizeDelaunay::Finite_vertices_iterator vit =
mesh_.finite_vertices_begin();
vit != mesh_.finite_vertices_end();
++vit
)
{
if (vit->referred())
{
vit->alignment() = alignments[referredPoints[referredI++]];
}
}
}
// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -0,0 +1,139 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::smoothAlignmentSolver
Description
SourceFiles
smoothAlignmentSolverI.H
smoothAlignmentSolver.C
smoothAlignmentSolverIO.C
\*---------------------------------------------------------------------------*/
#ifndef smoothAlignmentSolver_H
#define smoothAlignmentSolver_H
#include "cellShapeControlMesh.H"
#include "triadField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class smoothAlignmentSolver Declaration
\*---------------------------------------------------------------------------*/
class smoothAlignmentSolver
{
// Private data
cellShapeControlMesh& mesh_;
// Private Member Functions
template <class Triangulation, class Type>
tmp<Field<Type> > filterFarPoints
(
const Triangulation& mesh,
const Field<Type>& field
);
template <class Triangulation>
autoPtr<mapDistribute> buildMap
(
const Triangulation& mesh,
labelListList& pointPoints
);
template <class Triangulation>
autoPtr<mapDistribute> buildReferredMap
(
const Triangulation& mesh,
labelList& indices
);
template <class Triangulation>
tmp<triadField> buildAlignmentField(const Triangulation& mesh);
template <class Triangulation>
tmp<pointField> buildPointField(const Triangulation& mesh);
//- Disallow default bitwise copy construct
smoothAlignmentSolver(const smoothAlignmentSolver&);
//- Disallow default bitwise assignment
void operator=(const smoothAlignmentSolver&);
public:
// Constructors
//- Construct null
smoothAlignmentSolver(cellShapeControlMesh& mesh);
//- Destructor
~smoothAlignmentSolver();
// Member Functions
// Access
// Check
// Edit
//- Smooth the alignments on the mesh
void smoothAlignments(const label maxSmoothingIterations);
// Write
// Member Operators
// Friend Functions
// Friend Operators
// IOstream Operators
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -31,15 +31,15 @@ License
#include "vectorTools.H"
#include "IOmanip.H"
#include "indexedCellChecks.H"
#include "controlMeshRefinement.H"
#include "smoothAlignmentSolver.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(conformalVoronoiMesh, 0);
}
const Foam::scalar Foam::conformalVoronoiMesh::tolParallel = 1e-3;
@ -634,160 +634,6 @@ void Foam::conformalVoronoiMesh::insertInitialPoints()
}
Foam::scalar Foam::conformalVoronoiMesh::calculateLoadUnbalance() const
{
label nRealVertices = 0;
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
++vit
)
{
// Only store real vertices that are not feature vertices
if (vit->real() && !vit->featurePoint())
{
nRealVertices++;
}
}
scalar globalNRealVertices = returnReduce
(
nRealVertices,
sumOp<label>()
);
scalar unbalance = returnReduce
(
mag(1.0 - nRealVertices/(globalNRealVertices/Pstream::nProcs())),
maxOp<scalar>()
);
Info<< " Processor unbalance " << unbalance << endl;
return unbalance;
}
bool Foam::conformalVoronoiMesh::distributeBackground()
{
if (!Pstream::parRun())
{
return false;
}
Info<< nl << "Redistributing points" << endl;
timeCheck("Before distribute");
label iteration = 0;
scalar previousLoadUnbalance = 0;
while (true)
{
scalar maxLoadUnbalance = calculateLoadUnbalance();
if
(
maxLoadUnbalance <= cvMeshControls().maxLoadUnbalance()
|| maxLoadUnbalance <= previousLoadUnbalance
)
{
// If this is the first iteration, return false, if it was a
// subsequent one, return true;
return iteration != 0;
}
previousLoadUnbalance = maxLoadUnbalance;
Info<< " Total number of vertices before redistribution "
<< returnReduce(label(number_of_vertices()), sumOp<label>())
<< endl;
const fvMesh& bMesh = decomposition_().mesh();
volScalarField cellWeights
(
IOobject
(
"cellWeights",
bMesh.time().timeName(),
bMesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
bMesh,
dimensionedScalar("weight", dimless, 1e-2),
zeroGradientFvPatchScalarField::typeName
);
meshSearch cellSearch(bMesh, polyMesh::FACEPLANES);
labelList cellVertices(bMesh.nCells(), 0);
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
++vit
)
{
// Only store real vertices that are not feature vertices
if (vit->real() && !vit->featurePoint())
{
pointFromPoint v = topoint(vit->point());
label cellI = cellSearch.findCell(v);
if (cellI == -1)
{
// Pout<< "findCell conformalVoronoiMesh::distribute "
// << "findCell "
// << vit->type() << " "
// << vit->index() << " "
// << v << " "
// << cellI
// << " find nearest cellI ";
cellI = cellSearch.findNearestCell(v);
}
cellVertices[cellI]++;
}
}
forAll(cellVertices, cI)
{
// Give a small but finite weight for empty cells. Some
// decomposition methods have difficulty with integer overflows in
// the sum of the normalised weight field.
cellWeights.internalField()[cI] = max
(
cellVertices[cI],
1e-2
);
}
autoPtr<mapDistributePolyMesh> mapDist = decomposition_().distribute
(
cellWeights
);
cellShapeControl_.shapeControlMesh().distribute(decomposition_);
distribute();
timeCheck("After distribute");
iteration++;
}
return true;
}
void Foam::conformalVoronoiMesh::distribute()
{
if (!Pstream::parRun())
@ -871,13 +717,24 @@ void Foam::conformalVoronoiMesh::distribute()
void Foam::conformalVoronoiMesh::buildCellSizeAndAlignmentMesh()
{
cellShapeControl_.initialMeshPopulation(decomposition_);
controlMeshRefinement meshRefinement
(
cellShapeControl_
);
smoothAlignmentSolver meshAlignmentSmoother
(
cellShapeControl_.shapeControlMesh()
);
meshRefinement.initialMeshPopulation(decomposition_);
cellShapeControlMesh& cellSizeMesh = cellShapeControl_.shapeControlMesh();
if (Pstream::parRun())
{
cellSizeMesh.distribute(decomposition_);
distributeBackground(cellSizeMesh);
// cellSizeMesh.distribute(decomposition_);
}
label nMaxIter = readLabel
@ -892,40 +749,48 @@ void Foam::conformalVoronoiMesh::buildCellSizeAndAlignmentMesh()
for (label i = 0; i < nMaxIter; ++i)
{
// label nRemoved = cellSizeMesh.removePoints();
label nRemoved = 0;
reduce(nRemoved, sumOp<label>());
label nAdded = cellShapeControl_.refineMesh(decomposition_);
// label nAdded = 0;
label nAdded = meshRefinement.refineMesh(decomposition_);
//cellShapeControl_.refineMesh(decomposition_);
reduce(nAdded, sumOp<label>());
if (Pstream::parRun())
{
cellSizeMesh.distribute(decomposition_);
}
if (nRemoved + nAdded == 0)
{
break;
distributeBackground(cellSizeMesh);
}
Info<< " Iteration " << i
<< " Added = " << nAdded << " points"
<< ", Removed = " << nRemoved << " points"
<< endl;
if (nAdded == 0)
{
break;
}
}
cellShapeControl_.smoothMesh();
label maxSmoothingIterations = readLabel
(
cvMeshControls().cvMeshDict().subDict("motionControl").lookup
(
"maxSmoothingIterations"
)
);
meshAlignmentSmoother.smoothAlignments(maxSmoothingIterations);
Info<< "Background cell size and alignment mesh:" << endl;
cellSizeMesh.printInfo(Info);
if (cvMeshControls().objOutput())
Info<< "Triangulation is "
<< (cellSizeMesh.is_valid() ? "valid" : "not valid!" )
<< endl;
if (!Pstream::parRun())
{
cellSizeMesh.writeTriangulation();
//cellSizeMesh.writeTriangulation();
cellSizeMesh.write();
}
cellSizeMesh.printVertexInfo(Info);
}
@ -1102,11 +967,6 @@ void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
vit->targetCellSize(),
vit->alignment()
);
//vit->alignment() = tensor(1,0,0,0,1,0,0,0,1);
//vit->alignment() = requiredAlignment(pt);
//vit->targetCellSize() = cellShapeControls().cellSize(pt);
}
}
}
@ -1424,13 +1284,13 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
// Improve the guess that the backgroundMeshDecomposition makes with the
// initial positions. Use before building the surface conformation to
// better balance the surface conformation load.
distributeBackground();
distributeBackground(*this);
buildSurfaceConformation();
// The introduction of the surface conformation may have distorted the
// balance of vertices, distribute if necessary.
distributeBackground();
distributeBackground(*this);
if (Pstream::parRun())
{
@ -1701,22 +1561,22 @@ void Foam::conformalVoronoiMesh::move()
delta
);
if (targetFaceArea == 0)
{
Pout<< vA->info() << vB->info();
Cell_handle ch = locate(vA->point());
if (is_infinite(ch))
{
Pout<< "vA " << vA->targetCellSize() << endl;
}
ch = locate(vB->point());
if (is_infinite(ch))
{
Pout<< "vB " << vB->targetCellSize() << endl;
}
}
// if (targetFaceArea == 0)
// {
// Pout<< vA->info() << vB->info();
//
// Cell_handle ch = locate(vA->point());
// if (is_infinite(ch))
// {
// Pout<< "vA " << vA->targetCellSize() << endl;
// }
//
// ch = locate(vB->point());
// if (is_infinite(ch))
// {
// Pout<< "vB " << vB->targetCellSize() << endl;
// }
// }
delta *= faceAreaWeightModel_->faceAreaWeight
(

View File

@ -34,6 +34,7 @@ SourceFiles
conformalVoronoiMeshFeaturePoints.C
conformalVoronoiMeshFeaturePointSpecialisations.C
conformalVoronoiMeshCalcDualMesh.C
conformalVoronoiMeshTemplates.C
\*---------------------------------------------------------------------------*/

View File

@ -69,7 +69,7 @@ void Foam::conformalVoronoiMesh::conformToSurface()
// Rebuild, insert and store new surface conformation
buildSurfaceConformation();
if (distributeBackground())
if (distributeBackground(*this))
{
if (Pstream::parRun())
{

View File

@ -613,17 +613,9 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints()
}
label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
reduce(nFeatureVertices, sumOp<label>());
if (Pstream::parRun())
{
reduce(nFeatureVertices, sumOp<label>());
}
if (nFeatureVertices > 0)
{
Info<< " Inserted " << nFeatureVertices
<< " feature vertices" << endl;
}
Info<< " Inserted " << nFeatureVertices << " feature vertices" << endl;
featureVertices_.clear();
featureVertices_.setSize(pts.size());

View File

@ -92,6 +92,7 @@ surface2.nas
// Output the proximity of feature points and edges to each other
featureProximity no;
// The maximum search distance to use when looking for other feature
// points and edges
maxFeatureProximity 1;

View File

@ -2,7 +2,7 @@
# ========= |
# \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
# \\ / O peration |
# \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
# \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
# \\/ M anipulation |
#------------------------------------------------------------------------------
# License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -75,8 +75,6 @@ const triad triad::unset
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //