Commit 60f07d1f authored by andy's avatar andy
Browse files

Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

parents ef2dcd9a b66b9637
......@@ -295,14 +295,14 @@ void Foam::controlMeshRefinement::initialMeshPopulation
shapeController_.minimumCellSize()
);
}
else if (maxPriority == controlFunction.maxPriority())
{
vertices[vI].targetCellSize() = max
(
min(sizes[vI], size),
shapeController_.minimumCellSize()
);
}
// else if (maxPriority == controlFunction.maxPriority())
// {
// vertices[vI].targetCellSize() = max
// (
// min(sizes[vI], size),
// shapeController_.minimumCellSize()
// );
// }
else
{
vertices[vI].targetCellSize() = max
......
......@@ -1487,8 +1487,13 @@ void Foam::conformalVoronoiMesh::move()
// Save displacements to file.
if (foamyHexMeshControls().objOutput() && time().outputTime())
{
Pout<< "Writing point displacement vectors to file." << endl;
OFstream str("displacements_" + runTime_.timeName() + ".obj");
Info<< "Writing point displacement vectors to file." << endl;
OFstream str
(
time().path()
+ "displacements_" + runTime_.timeName()
+ ".obj"
);
for
(
......
......@@ -550,53 +550,130 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
}
}
// for
// (
// Delaunay::Finite_edges_iterator eit = finite_edges_begin();
// eit != finite_edges_end();
// ++eit
// )
// {
// Cell_handle c = eit->first;
// Vertex_handle vA = c->vertex(eit->second);
// Vertex_handle vB = c->vertex(eit->third);
//
// if
// (
// vA->referred()
// || vB->referred()
// )
// {
// continue;
// }
//
// if
// (
// (vA->internalPoint() && vA->referred())
// || (vB->internalPoint() && vB->referred())
// )
// {
// continue;
// }
//
// if
// (
// (vA->internalPoint() && vB->externalBoundaryPoint())
// || (vB->internalPoint() && vA->externalBoundaryPoint())
// || (vA->internalBoundaryPoint() && vB->internalBoundaryPoint())
// )
// {
// pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
// pointIndexHit surfHit;
// label hitSurface;
//
// geometryToConformTo_.findSurfaceNearest
// (
// (
// vA->internalPoint()
// ? topoint(vA->point())
// : topoint(vB->point())
// ),
// magSqr(topoint(vA->point()) - topoint(vB->point())),
// surfHit,
// hitSurface
// );
//
// if (surfHit.hit())
// {
// surfaceIntersections.append
// (
// pointIndexHitAndFeature(surfHit, hitSurface)
// );
//
// addSurfaceAndEdgeHits
// (
// (
// vA->internalPoint()
// ? topoint(vA->point())
// : topoint(vB->point())
// ),
// surfaceIntersections,
// surfacePtReplaceDistCoeffSqr,
// edgeSearchDistCoeffSqr,
// surfaceHits,
// featureEdgeHits,
// surfaceToTreeShape,
// edgeToTreeShape,
// surfacePtToEdgePtDist,
// false
// );
// }
// }
// }
for
(
Delaunay::Finite_edges_iterator eit = finite_edges_begin();
eit != finite_edges_end();
++eit
Delaunay::Finite_cells_iterator cit = finite_cells_begin();
cit != finite_cells_end();
++cit
)
{
Cell_handle c = eit->first;
Vertex_handle vA = c->vertex(eit->second);
Vertex_handle vB = c->vertex(eit->third);
if
(
vA->referred()
|| vB->referred()
)
label nInternal = 0;
for (label vI = 0; vI < 4; vI++)
{
continue;
if (cit->vertex(vI)->internalPoint())
{
nInternal++;
}
}
if
(
(vA->internalPoint() && vA->referred())
|| (vB->internalPoint() && vB->referred())
)
if (nInternal == 1 && cit->hasBoundaryPoint())
//if (cit->boundaryDualVertex() && !cit->hasReferredPoint())
{
continue;
}
const Foam::point& pt = cit->dual();
const Foam::point cellCentre =
topoint
(
CGAL::centroid
(
CGAL::Tetrahedron_3<baseK>
(
cit->vertex(0)->point(),
cit->vertex(1)->point(),
cit->vertex(2)->point(),
cit->vertex(3)->point()
)
)
);
if
(
(vA->internalPoint() && vB->externalBoundaryPoint())
|| (vB->internalPoint() && vA->externalBoundaryPoint())
)
{
pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
pointIndexHit surfHit;
label hitSurface;
geometryToConformTo_.findSurfaceNearest
geometryToConformTo_.findSurfaceNearestIntersection
(
(
vA->internalPoint()
? topoint(vA->point())
: topoint(vB->point())
),
magSqr(topoint(vA->point()) - topoint(vB->point())),
cellCentre,
pt,
surfHit,
hitSurface
);
......@@ -610,11 +687,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
addSurfaceAndEdgeHits
(
(
vA->internalPoint()
? topoint(vA->point())
: topoint(vB->point())
),
pt,
surfaceIntersections,
surfacePtReplaceDistCoeffSqr,
edgeSearchDistCoeffSqr,
......@@ -1613,6 +1686,7 @@ void Foam::conformalVoronoiMesh::limitDisplacement
// Do not allow infinite recursion
if (callCount > 7)
{
displacement = vector::zero;
return;
}
......@@ -1661,6 +1735,7 @@ void Foam::conformalVoronoiMesh::limitDisplacement
if (magSqr(pt - surfHit.hitPoint()) <= searchDistanceSqr)
{
// Cannot limit displacement, point closer than tolerance
displacement = vector::zero;
return;
}
}
......@@ -1748,17 +1823,19 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
const Foam::point& pt = pHit.first().hitPoint();
pointIndexHit closePoint;
const bool closeToSurfacePt = pointIsNearSurfaceLocation(pt, closePoint);
if
(
closeToSurfacePt
&& mag(pt - closePoint.hitPoint()) > pointPairDistance(pt)
&& (
magSqr(pt - closePoint.hitPoint())
> sqr(pointPairDistance(pt))
)
)
{
const scalar cosAngle
= angleBetweenSurfacePoints(pt, closePoint.hitPoint());
const scalar cosAngle =
angleBetweenSurfacePoints(pt, closePoint.hitPoint());
// @todo make this tolerance run-time selectable?
if (cosAngle < searchAngleOppositeSurface)
......@@ -1768,11 +1845,6 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
const scalar searchDist = targetCellSize(closePoint.hitPoint());
if (searchDist < SMALL)
{
Pout<< "WARNING: SMALL CELL SIZE" << endl;
}
geometryToConformTo_.findSurfaceNearest
(
closePoint.hitPoint(),
......@@ -1789,15 +1861,15 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
norm
);
const vector nA = norm[0];
const vector& nA = norm[0];
pointIndexHit oppositeHit;
label oppositeSurfaceHit = -1;
geometryToConformTo_.findSurfaceNearestIntersection
(
closePoint.hitPoint() + SMALL*nA,
closePoint.hitPoint() + mag(pt - closePoint.hitPoint())*nA,
closePoint.hitPoint() + 0.5*pointPairDistance(pt)*nA,
closePoint.hitPoint() + 5*targetCellSize(pt)*nA,
oppositeHit,
oppositeSurfaceHit
);
......@@ -1808,21 +1880,10 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
pHit.first() = oppositeHit;
pHit.second() = oppositeSurfaceHit;
// appendToSurfacePtTree(pHit.first().hitPoint());
// surfaceToTreeShape.append
// (
// existingSurfacePtLocations_.size() - 1
// );
return !closeToSurfacePt;
}
}
}
else
{
// appendToSurfacePtTree(pt);
// surfaceToTreeShape.append(existingSurfacePtLocations_.size() - 1);
}
return closeToSurfacePt;
}
......
......@@ -230,13 +230,13 @@ inline void Foam::conformalVoronoiMesh::createPointPair
{
vector ppDistn = ppDist*n;
const Foam::point internalPt = surfPt - ppDistn;
const Foam::point externalPt = surfPt + ppDistn;
// const Foam::point internalPt = surfPt - ppDistn;
// const Foam::point externalPt = surfPt + ppDistn;
bool internalInside = geometryToConformTo_.inside(internalPt);
bool externalOutside = geometryToConformTo_.outside(externalPt);
// bool internalInside = geometryToConformTo_.inside(internalPt);
// bool externalOutside = geometryToConformTo_.outside(externalPt);
if (internalInside && externalOutside)
// if (internalInside && externalOutside)
{
pts.append
(
......@@ -259,15 +259,24 @@ inline void Foam::conformalVoronoiMesh::createPointPair
Pstream::myProcNo()
)
);
// if (ptPair)
// {
// ptPairs_.addPointPair
// (
// pts[pts.size() - 2].index(),
// pts[pts.size() - 1].index() // external 0 -> slave
// );
// }
}
else
{
Info<< "Warning: point pair not inside/outside" << nl
<< " surfPt = " << surfPt << nl
<< " internal = " << internalPt << " " << internalInside << nl
<< " external = " << externalPt << " " << externalOutside
<< endl;
}
// else
// {
// Info<< "Warning: point pair not inside/outside" << nl
// << " surfPt = " << surfPt << nl
// << " internal = " << internalPt << " " << internalInside << nl
// << " external = " << externalPt << " " << externalOutside
// << endl;
// }
}
......
......@@ -183,6 +183,7 @@ public:
inline bool hasSeedPoint() const;
inline bool hasInternalPoint() const;
inline bool hasBoundaryPoint() const;
inline bool hasConstrainedPoint() const;
......
......@@ -228,6 +228,19 @@ inline bool CGAL::indexedCell<Gt, Cb>::hasInternalPoint() const
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::hasBoundaryPoint() const
{
return
(
this->vertex(0)->boundaryPoint()
|| this->vertex(1)->boundaryPoint()
|| this->vertex(2)->boundaryPoint()
|| this->vertex(3)->boundaryPoint()
);
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::hasConstrainedPoint() const
{
......
......@@ -683,14 +683,18 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
vector hitDir = info[0].rawPoint() - samplePts[i];
hitDir /= mag(hitDir) + SMALL;
if
pointIndexHit surfHit;
label hitSurface;
findSurfaceNearestIntersection
(
findSurfaceAnyIntersection
(
samplePts[i],
info[0].rawPoint() - 1e-3*mag(hitDir)*(hitDir)
)
)
samplePts[i],
info[0].rawPoint(),
surfHit,
hitSurface
);
if (surfHit.hit() && hitSurface != surfaces_[s])
{
continue;
}
......
......@@ -37,6 +37,90 @@ namespace Foam
defineTypeNameAndDebug(rayShooting, 0);
addToRunTimeSelectionTable(initialPointsMethod, rayShooting, dictionary);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void rayShooting::splitLine
(
const line<point, point>& l,
const scalar& pert,
DynamicList<Vb::Point>& initialPoints
) const
{
Foam::point midPoint(l.centre());
const scalar localCellSize(cellShapeControls().cellSize(midPoint));
const scalar lineLength(l.mag());
const scalar minDistFromSurfaceSqr
(
minimumSurfaceDistanceCoeffSqr_
*sqr(localCellSize)
);
if
(
magSqr(midPoint - l.start()) > minDistFromSurfaceSqr
&& magSqr(midPoint - l.end()) > minDistFromSurfaceSqr
)
{
// Add extra points if line length is much bigger than local cell size
// if (lineLength > 4.0*localCellSize)
// {
// splitLine
// (
// line<point, point>(l.start(), midPoint),
// pert,
// initialPoints
// );
//
// splitLine
// (
// line<point, point>(midPoint, l.end()),
// pert,
// initialPoints
// );
// }
if (randomiseInitialGrid_)
{
Foam::point newPt
(
midPoint.x() + pert*(rndGen().scalar01() - 0.5),
midPoint.y() + pert*(rndGen().scalar01() - 0.5),
midPoint.z() + pert*(rndGen().scalar01() - 0.5)
);
if
(
!geometryToConformTo().findSurfaceAnyIntersection
(
midPoint,
newPt
)
)
{
midPoint = newPt;
}
else
{
WarningIn
(
"rayShooting::splitLine"
"("
" const line<point,point>&,"
" const scalar&,"
" DynamicList<Vb::Point>&"
")"
) << "Point perturbation crosses a surface. Not inserting."
<< endl;
}
}
initialPoints.append(toPoint(midPoint));
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
rayShooting::rayShooting
......@@ -75,7 +159,7 @@ List<Vb::Point> rayShooting::initialPoints() const
const searchableSurfaces& surfaces = geometryToConformTo().geometry();
const labelList& surfacesToConformTo = geometryToConformTo().surfaces();
const scalar maxRayLength = surfaces.bounds().mag();
const scalar maxRayLength(surfaces.bounds().mag());
// Initialise points list
label initialPointsSize = 0;
......@@ -94,8 +178,7 @@ List<Vb::Point> rayShooting::initialPoints() const
const pointField& faceCentres = faceCentresTmp();
Info<< " Shoot rays from " << s.name() << nl
<< " nRays = " << faceCentres.size() << endl;
<< " nRays = " << faceCentres.size() << endl;
forAll(faceCentres, fcI)
{
......@@ -110,9 +193,11 @@ List<Vb::Point> rayShooting::initialPoints() const
continue;
}
const scalar pert =
const scalar pert
(
randomPerturbationCoeff_
*cellShapeControls().cellSize(fC);
*cellShapeControls().cellSize(fC)
);
pointIndexHit surfHitStart;
label hitSurfaceStart;
......@@ -181,27 +266,12 @@ List<Vb::Point> rayShooting::initialPoints() const
}
}
Foam::point midPoint(l.centre());
const scalar minDistFromSurfaceSqr =
minimumSurfaceDistanceCoeffSqr_
*sqr(cellShapeControls().cellSize(midPoint));
if (randomiseInitialGrid_)
{
midPoint.x() += pert*(rndGen().scalar01() - 0.5);
midPoint.y() += pert*(rndGen().scalar01() - 0.5);
midPoint.z() += pert*(rndGen().scalar01() - 0.5);
}