ENH: Add forcePointInsertion option to controlMeshRefinement

This commit is contained in:
laurence
2013-04-11 20:28:08 +01:00
parent 26a42b37a0
commit 4af889f860
8 changed files with 171 additions and 88 deletions

View File

@ -103,7 +103,11 @@ Foam::scalar Foam::cellShapeControl::cellSize(const point& pt) const
scalar size = 0;
if (shapeControlMesh_.dimension() < 3 || shapeControlMesh_.is_infinite(ch))
if (shapeControlMesh_.dimension() < 3)
{
size = sizeAndAlignment_.cellSize(pt);
}
else if (shapeControlMesh_.is_infinite(ch))
{
// if (nFarPoints != 0)
// {
@ -204,7 +208,7 @@ Foam::tensor Foam::cellShapeControl::cellAlignment(const point& pt) const
tri.normalize();
tri.orthogonalize();
// tri = tri.sortxyz();
tri = tri.sortxyz();
alignment = tri;
}
@ -282,7 +286,7 @@ void Foam::cellShapeControl::cellSizeAndAlignment
tri.normalize();
tri.orthogonalize();
// tri = tri.sortxyz();
tri = tri.sortxyz();
alignment = tri;
@ -298,21 +302,51 @@ void Foam::cellShapeControl::cellSizeAndAlignment
for (label dir = 0; dir < 3; dir++)
{
const triad v = alignment;
triad v = alignment;
if (!v.set(dir) || size == 0)
{
FatalErrorIn
(
"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);
// Force orthogonalization of triad.
scalar dotProd = GREAT;
if (dir == 0)
{
dotProd = v[1] & v[2];
v[dir] = v[1] ^ v[2];
}
if (dir == 1)
{
dotProd = v[0] & v[2];
v[dir] = v[0] ^ v[2];
}
if (dir == 2)
{
dotProd = v[0] & v[1];
v[dir] = v[0] ^ v[1];
}
v.normalize();
v.orthogonalize();
Pout<< "Dot prod = " << dotProd << endl;
Pout<< "Alignment = " << v << endl;
alignment = v;
// FatalErrorIn
// (
// "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

@ -536,6 +536,7 @@ void Foam::cellShapeControlMesh::barycentricCoords
) const
{
// Use the previous cell handle as a hint on where to start searching
// Giving a hint causes strange errors...
ch = locate(toPoint<Point>(pt));
if (dimension() > 2 && !is_infinite(ch))
@ -672,8 +673,6 @@ void Foam::cellShapeControlMesh::distribute
true
);
Info<< nl << " Inserted distributed background tessellation" << endl;
sync(decomposition.procBounds());
Info<< " Total number of vertices after redistribution "

View File

@ -52,7 +52,15 @@ Foam::cellSizeAndAlignmentControl::cellSizeAndAlignmentControl
:
runTime_(runTime),
defaultCellSize_(defaultCellSize),
name_(name)
name_(name),
forceInitialPointInsertion_
(
controlFunctionDict.lookupOrDefault<Switch>
(
"forceInitialPointInsertion",
"off"
)
)
{}

View File

@ -59,6 +59,8 @@ protected:
const scalar& defaultCellSize_;
Switch forceInitialPointInsertion_;
private:
@ -146,6 +148,11 @@ public:
return name_;
}
const Switch& forceInitialPointInsertion() const
{
return forceInitialPointInsertion_;
}
// Query

View File

@ -50,7 +50,7 @@ bool Foam::cellSizeAndAlignmentControls::evalCellSizeFunctions
{
// Maintain priority of current hit. Initialise so it always goes
// through at least once.
label previousPriority = -1;
label previousPriority = labelMin;
forAll(controlFunctions_, i)
{

View File

@ -186,62 +186,94 @@ Foam::searchableSurfaceControl::searchableSurfaceControl
(
controlFunctionDict,
searchableSurface_,
defaultCellSize_
defaultCellSize_,
labelList()
)
);
Info<< decrIndent;
PtrList<cellSizeFunction> regionCellSizeFunctions
(
searchableSurface_.regions().size()
);
PtrList<cellSizeFunction> regionCellSizeFunctions;
label functionI = 0;
DynamicList<label> defaultCellSizeRegions;
label nRegionCellSizeFunctions = 0;
// Loop over regions - if any entry is not specified they should
// inherit values from the parent surface.
if (controlFunctionDict.found("regions"))
{
const dictionary& regionsDict = controlFunctionDict.subDict("regions");
const wordList& regionNames = geometryToConformTo.patchNames();
const wordList& regionNames = geometryToConformTo_.patchNames();
label nRegions = regionsDict.size();
regionCellSizeFunctions.setSize(nRegions + 1);
defaultCellSizeRegions.setCapacity(nRegions);
forAll(regionNames, regionI)
{
const word& regionName = regionNames[regionI];
label regionID = geometryToConformTo_.geometry().findSurfaceRegionID
(
this->name(),
regionName
);
if (regionsDict.found(regionName))
{
// Get the dictionary for region
const dictionary& regionDict = regionsDict.subDict(regionName);
Info<< indent << "Region " << regionName << " settings:"
Info<< indent << "Region " << regionName
<< " (ID = " << regionID << ")" << " settings:"
<< endl;
Info<< incrIndent;
regionCellSizeFunctions.set
(
functionI,
nRegionCellSizeFunctions,
cellSizeFunction::New
(
regionDict,
searchableSurface_,
defaultCellSize_
defaultCellSize_,
labelList(1, regionID)
)
);
Info<< decrIndent;
functionI++;
nRegionCellSizeFunctions++;
}
else
{
// Add to default list
defaultCellSizeRegions.append(regionID);
}
}
}
regionCellSizeFunctions.setSize(functionI);
if (functionI > 0)
if (defaultCellSizeRegions.empty() && !regionCellSizeFunctions.empty())
{
cellSizeFunctions_.transfer(regionCellSizeFunctions);
}
else if (nRegionCellSizeFunctions > 0)
{
regionCellSizeFunctions.set
(
nRegionCellSizeFunctions,
cellSizeFunction::New
(
controlFunctionDict,
searchableSurface_,
defaultCellSize_,
labelList()
)
);
cellSizeFunctions_.transfer(regionCellSizeFunctions);
}
forAll(cellSizeFunctions_, funcI)
{
@ -252,6 +284,9 @@ Foam::searchableSurfaceControl::searchableSurfaceControl
maxPriority_ = funcPriority;
}
}
Info<< nl << "There are " << cellSizeFunctions_.size()
<< " region control functions" << endl;
}
@ -387,7 +422,7 @@ bool Foam::searchableSurfaceControl::cellSize
continue;
}
scalar sizeI = GREAT;
scalar sizeI = -1;
if (sizeFunc.cellSize(pt, sizeI))
{
@ -403,15 +438,15 @@ bool Foam::searchableSurfaceControl::cellSize
else
{
cellSize = sizeI;
priority = sizeFunc.priority();
}
if (debug)
{
Info<< "sizeI " << sizeI
<<" minSize " << cellSize << endl;
Info<< " sizeI " << sizeI
<<" minSize " << cellSize << endl;
}
priority = sizeFunc.priority();
}
}

View File

@ -77,13 +77,8 @@ bool Foam::controlMeshRefinement::detectEdge
Foam::point a(startPt);
Foam::point b(endPt);
// scalar cellSizeA = sizeControls_.cellSize(a);
// scalar cellSizeB = sizeControls_.cellSize(b);
Foam::point midPoint = (a + b)/2.0;
// scalar cellSizeMid = sizeControls_.cellSize(midPoint);
label nIterations = 0;
while (true)
@ -106,10 +101,10 @@ bool Foam::controlMeshRefinement::detectEdge
scalar cellSizeA = sizeControls_.cellSize(a);
scalar cellSizeB = sizeControls_.cellSize(b);
if (magSqr(cellSizeA - cellSizeB) < 1e-6)
{
return false;
}
// if (magSqr(cellSizeA - cellSizeB) < 1e-6)
// {
// return false;
// }
scalar cellSizeMid = sizeControls_.cellSize(midPoint);
@ -164,16 +159,12 @@ bool Foam::controlMeshRefinement::detectEdge
if (magSqr(secondDerivative1) > magSqr(secondDerivative2))
{
b = midPoint;
// cellSizeB = cellSizeMid;
midPoint = midPoint1;
// cellSizeMid = cellSizeMid1;
}
else
{
a = midPoint;
// cellSizeA = cellSizeMid;
midPoint = midPoint2;
// cellSizeMid = cellSizeMid2;
}
}
}
@ -186,8 +177,8 @@ Foam::pointHit Foam::controlMeshRefinement::findDiscontinuities
{
pointHit p(point::max);
const scalar tolSqr = sqr(1e-1);
const scalar secondDerivTolSqr = sqr(1e-1);
const scalar tolSqr = sqr(1e-3);
const scalar secondDerivTolSqr = sqr(1e-3);
detectEdge
(
@ -242,16 +233,16 @@ void Foam::controlMeshRefinement::initialMeshPopulation
}
else
{
overallBoundBox.set
(
new boundBox(geometryToConformTo_.geometry().bounds())
);
mesh_.insertBoundingPoints
(
overallBoundBox(),
sizeControls_
);
// overallBoundBox.set
// (
// new boundBox(geometryToConformTo_.geometry().bounds())
// );
//
// mesh_.insertBoundingPoints
// (
// overallBoundBox(),
// sizeControls_
// );
}
const PtrList<cellSizeAndAlignmentControl>& controlFunctions =
@ -262,8 +253,12 @@ void Foam::controlMeshRefinement::initialMeshPopulation
const cellSizeAndAlignmentControl& controlFunction =
controlFunctions[fI];
const Switch& forceInsertion =
controlFunction.forceInitialPointInsertion();
Info<< "Inserting points from " << controlFunction.name()
<< " (" << controlFunction.type() << ")" << endl;
Info<< " Force insertion is " << forceInsertion.asText() << endl;
pointField pts;
scalarField sizes;
@ -271,6 +266,8 @@ void Foam::controlMeshRefinement::initialMeshPopulation
controlFunction.initialVertices(pts, sizes, alignments);
Info<< " Got initial vertices list" << endl;
List<Vb> vertices(pts.size());
// Clip the minimum size
@ -279,12 +276,18 @@ void Foam::controlMeshRefinement::initialMeshPopulation
vertices[vI] = Vb(pts[vI], Vb::vtInternalNearBoundary);
vertices[vI].targetCellSize() = max
(
sizes[vI],
min
(
sizes[vI],
shapeController_.cellSize(pts[vI])
),
shapeController_.minimumCellSize()
);
vertices[vI].alignment() = alignments[vI];
}
Info<< " Clipped minimum size" << endl;
pts.clear();
sizes.clear();
alignments.clear();
@ -293,19 +296,25 @@ void Foam::controlMeshRefinement::initialMeshPopulation
PackedBoolList keepVertex(vertices.size(), true);
if (Pstream::parRun())
forAll(vertices, vI)
{
forAll(vertices, vI)
{
const bool onProc = decomposition().positionOnThisProcessor
(
topoint(vertices[vI].point())
);
bool keep = true;
if (!onProc)
{
keepVertex[vI] = false;
}
pointFromPoint pt = topoint(vertices[vI].point());
if (Pstream::parRun())
{
keep = decomposition().positionOnThisProcessor(pt);
}
if (geometryToConformTo_.wellOutside(pt, SMALL))
{
keep = false;
}
if (!keep)
{
keepVertex[vI] = false;
}
}
@ -313,6 +322,8 @@ void Foam::controlMeshRefinement::initialMeshPopulation
const label preInsertedSize = mesh_.number_of_vertices();
Info<< " Check sizes" << endl;
forAll(vertices, vI)
{
bool insertPoint = false;
@ -368,7 +379,7 @@ void Foam::controlMeshRefinement::initialMeshPopulation
insertPoint = true;
}
if (insertPoint)
if (forceInsertion || insertPoint)
{
mesh_.insert
(
@ -492,9 +503,6 @@ Foam::label Foam::controlMeshRefinement::refineMesh
linePointRef l(ptA, ptB);
// Info<< "Edge " << l << endl;
// Info<< vA->info() << vB->info();
const pointHit hitPt = findDiscontinuities(l);
if (hitPt.hit())

View File

@ -48,7 +48,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class controlMeshRefinement Declaration
Class controlMeshRefinement Declaration
\*---------------------------------------------------------------------------*/
class controlMeshRefinement
@ -124,10 +124,6 @@ public:
// Member Functions
// Access
// Check
// Edit
void initialMeshPopulation
@ -139,10 +135,6 @@ public:
(
const autoPtr<backgroundMeshDecomposition>& decomposition
);
// Write
};