mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
fc541cee40
- fixed indexing error (crashes under windows). No crash but still incorrect on other systems
1823 lines
55 KiB
C
1823 lines
55 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2015 OpenFOAM Foundation
|
|
Copyright (C) 2015-2020 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
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 "meshRefinement.H"
|
|
#include "Time.H"
|
|
#include "refinementSurfaces.H"
|
|
#include "refinementFeatures.H"
|
|
#include "shellSurfaces.H"
|
|
#include "triSurfaceMesh.H"
|
|
#include "treeDataCell.H"
|
|
#include "searchableSurfaces.H"
|
|
#include "DynamicField.H"
|
|
#include "transportData.H"
|
|
#include "FaceCellWave.H"
|
|
#include "volFields.H"
|
|
#include "zeroGradientFvPatchFields.H"
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
Foam::label Foam::meshRefinement::markSurfaceGapRefinement
|
|
(
|
|
const scalar planarCos,
|
|
|
|
const label nAllowRefine,
|
|
const labelList& neiLevel,
|
|
const pointField& neiCc,
|
|
|
|
labelList& refineCell,
|
|
label& nRefine
|
|
) const
|
|
{
|
|
const labelList& cellLevel = meshCutter_.cellLevel();
|
|
const pointField& cellCentres = mesh_.cellCentres();
|
|
|
|
// Get the gap level for the shells
|
|
const labelList maxLevel(shells_.maxGapLevel());
|
|
|
|
label oldNRefine = nRefine;
|
|
|
|
if (max(maxLevel) > 0)
|
|
{
|
|
// Use cached surfaceIndex_ to detect if any intersection. If so
|
|
// re-intersect to determine level wanted.
|
|
|
|
// Collect candidate faces
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
labelList testFaces(getRefineCandidateFaces(refineCell));
|
|
|
|
// Collect segments
|
|
// ~~~~~~~~~~~~~~~~
|
|
|
|
pointField start(testFaces.size());
|
|
pointField end(testFaces.size());
|
|
{
|
|
labelList minLevel(testFaces.size());
|
|
calcCellCellRays
|
|
(
|
|
neiCc,
|
|
neiLevel,
|
|
testFaces,
|
|
start,
|
|
end,
|
|
minLevel
|
|
);
|
|
}
|
|
|
|
|
|
// Collect cells to test for inside/outside in shell
|
|
labelList cellToCompact(mesh_.nCells(), -1);
|
|
labelList bFaceToCompact(mesh_.nBoundaryFaces(), -1);
|
|
labelList gapShell;
|
|
List<FixedList<label, 3>> shellGapInfo;
|
|
List<volumeType> shellGapMode;
|
|
{
|
|
DynamicField<point> compactToCc(mesh_.nCells()/10);
|
|
DynamicList<label> compactToLevel(compactToCc.capacity());
|
|
forAll(testFaces, i)
|
|
{
|
|
label faceI = testFaces[i];
|
|
label own = mesh_.faceOwner()[faceI];
|
|
if (cellToCompact[own] == -1)
|
|
{
|
|
cellToCompact[own] = compactToCc.size();
|
|
compactToCc.append(cellCentres[own]);
|
|
compactToLevel.append(cellLevel[own]);
|
|
}
|
|
if (mesh_.isInternalFace(faceI))
|
|
{
|
|
label nei = mesh_.faceNeighbour()[faceI];
|
|
if (cellToCompact[nei] == -1)
|
|
{
|
|
cellToCompact[nei] = compactToCc.size();
|
|
compactToCc.append(cellCentres[nei]);
|
|
compactToLevel.append(cellLevel[nei]);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
label bFaceI = faceI - mesh_.nInternalFaces();
|
|
if (bFaceToCompact[bFaceI] == -1)
|
|
{
|
|
bFaceToCompact[bFaceI] = compactToCc.size();
|
|
compactToCc.append(neiCc[bFaceI]);
|
|
compactToLevel.append(neiLevel[bFaceI]);
|
|
}
|
|
}
|
|
}
|
|
|
|
shells_.findHigherGapLevel
|
|
(
|
|
compactToCc,
|
|
compactToLevel,
|
|
|
|
gapShell,
|
|
shellGapInfo,
|
|
shellGapMode
|
|
);
|
|
}
|
|
|
|
|
|
//const fileName dir(mesh_.time().path()/timeName());
|
|
//if (debug)
|
|
//{
|
|
// mkDir(dir);
|
|
// OBJstream insideStr(dir/"insideShell.obj");
|
|
// OBJstream outsideStr(dir/"outsideShell.obj");
|
|
// Pout<< "Writing points to:" << nl
|
|
// << " inside : " << insideStr.name() << nl
|
|
// << " outside: " << outsideStr.name() << nl
|
|
// << endl;
|
|
//
|
|
// forAll(cellToCompact, celli)
|
|
// {
|
|
// const label compacti = cellToCompact[celli];
|
|
//
|
|
// if (compacti != -1)
|
|
// {
|
|
// if (gapShell[compacti] != -1)
|
|
// {
|
|
// insideStr.write(mesh_.cellCentres()[celli]);
|
|
// }
|
|
// else
|
|
// {
|
|
// outsideStr.write(mesh_.cellCentres()[celli]);
|
|
// }
|
|
// }
|
|
// }
|
|
// forAll(bFaceToCompact, bFacei)
|
|
// {
|
|
// const label compacti = bFaceToCompact[bFacei];
|
|
// if (compacti != -1)
|
|
// {
|
|
// if (gapShell[compacti] != -1)
|
|
// {
|
|
// insideStr.write(neiCc[bFacei]);
|
|
// }
|
|
// else
|
|
// {
|
|
// outsideStr.write(neiCc[bFacei]);
|
|
// }
|
|
// }
|
|
// }
|
|
//}
|
|
|
|
|
|
const List<FixedList<label, 3>>& extendedGapLevel =
|
|
surfaces_.extendedGapLevel();
|
|
const List<volumeType>& extendedGapMode =
|
|
surfaces_.extendedGapMode();
|
|
const boolList& extendedGapSelf = surfaces_.gapSelf();
|
|
|
|
labelList ccSurface1;
|
|
List<pointIndexHit> ccHit1;
|
|
labelList ccRegion1;
|
|
vectorField ccNormal1;
|
|
{
|
|
labelList ccSurface2;
|
|
List<pointIndexHit> ccHit2;
|
|
labelList ccRegion2;
|
|
vectorField ccNormal2;
|
|
|
|
surfaces_.findNearestIntersection
|
|
(
|
|
identity(surfaces_.surfaces().size()),
|
|
start,
|
|
end,
|
|
|
|
ccSurface1,
|
|
ccHit1,
|
|
ccRegion1,
|
|
ccNormal1,
|
|
|
|
ccSurface2,
|
|
ccHit2,
|
|
ccRegion2,
|
|
ccNormal2
|
|
);
|
|
}
|
|
|
|
start.clear();
|
|
end.clear();
|
|
|
|
DynamicField<point> rayStart(2*ccSurface1.size());
|
|
DynamicField<point> rayEnd(2*ccSurface1.size());
|
|
DynamicField<scalar> gapSize(2*ccSurface1.size());
|
|
|
|
DynamicField<point> rayStart2(2*ccSurface1.size());
|
|
DynamicField<point> rayEnd2(2*ccSurface1.size());
|
|
DynamicField<scalar> gapSize2(2*ccSurface1.size());
|
|
|
|
DynamicList<label> cellMap(2*ccSurface1.size());
|
|
DynamicList<label> compactMap(2*ccSurface1.size());
|
|
|
|
forAll(ccSurface1, i)
|
|
{
|
|
label surfI = ccSurface1[i];
|
|
|
|
if (surfI != -1)
|
|
{
|
|
label globalRegionI =
|
|
surfaces_.globalRegion(surfI, ccRegion1[i]);
|
|
|
|
label faceI = testFaces[i];
|
|
const point& surfPt = ccHit1[i].hitPoint();
|
|
|
|
label own = mesh_.faceOwner()[faceI];
|
|
if
|
|
(
|
|
cellToCompact[own] != -1
|
|
&& shellGapInfo[cellToCompact[own]][2] > 0
|
|
)
|
|
{
|
|
// Combine info from shell and surface
|
|
label compactI = cellToCompact[own];
|
|
FixedList<label, 3> gapInfo;
|
|
volumeType gapMode;
|
|
mergeGapInfo
|
|
(
|
|
shellGapInfo[compactI],
|
|
shellGapMode[compactI],
|
|
extendedGapLevel[globalRegionI],
|
|
extendedGapMode[globalRegionI],
|
|
|
|
gapInfo,
|
|
gapMode
|
|
);
|
|
|
|
const point& cc = cellCentres[own];
|
|
label nRays = generateRays
|
|
(
|
|
false,
|
|
surfPt,
|
|
ccNormal1[i],
|
|
gapInfo,
|
|
gapMode,
|
|
surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
|
|
cellLevel[own],
|
|
|
|
rayStart,
|
|
rayEnd,
|
|
gapSize,
|
|
|
|
rayStart2,
|
|
rayEnd2,
|
|
gapSize2
|
|
);
|
|
for (label j = 0; j < nRays; j++)
|
|
{
|
|
cellMap.append(own);
|
|
compactMap.append(i);
|
|
}
|
|
}
|
|
if (mesh_.isInternalFace(faceI))
|
|
{
|
|
label nei = mesh_.faceNeighbour()[faceI];
|
|
if
|
|
(
|
|
cellToCompact[nei] != -1
|
|
&& shellGapInfo[cellToCompact[nei]][2] > 0
|
|
)
|
|
{
|
|
// Combine info from shell and surface
|
|
label compactI = cellToCompact[nei];
|
|
FixedList<label, 3> gapInfo;
|
|
volumeType gapMode;
|
|
mergeGapInfo
|
|
(
|
|
shellGapInfo[compactI],
|
|
shellGapMode[compactI],
|
|
extendedGapLevel[globalRegionI],
|
|
extendedGapMode[globalRegionI],
|
|
|
|
gapInfo,
|
|
gapMode
|
|
);
|
|
|
|
const point& cc = cellCentres[nei];
|
|
label nRays = generateRays
|
|
(
|
|
false,
|
|
surfPt,
|
|
ccNormal1[i],
|
|
gapInfo,
|
|
gapMode,
|
|
surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
|
|
cellLevel[nei],
|
|
|
|
rayStart,
|
|
rayEnd,
|
|
gapSize,
|
|
|
|
rayStart2,
|
|
rayEnd2,
|
|
gapSize2
|
|
);
|
|
for (label j = 0; j < nRays; j++)
|
|
{
|
|
cellMap.append(nei);
|
|
compactMap.append(i);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Note: on coupled face. What cell are we going to
|
|
// refine? We've got the neighbouring cell centre
|
|
// and level but we cannot mark it for refinement on
|
|
// this side...
|
|
label bFaceI = faceI - mesh_.nInternalFaces();
|
|
|
|
if
|
|
(
|
|
bFaceToCompact[bFaceI] != -1
|
|
&& shellGapInfo[bFaceToCompact[bFaceI]][2] > 0
|
|
)
|
|
{
|
|
// Combine info from shell and surface
|
|
label compactI = bFaceToCompact[bFaceI];
|
|
FixedList<label, 3> gapInfo;
|
|
volumeType gapMode;
|
|
mergeGapInfo
|
|
(
|
|
shellGapInfo[compactI],
|
|
shellGapMode[compactI],
|
|
extendedGapLevel[globalRegionI],
|
|
extendedGapMode[globalRegionI],
|
|
|
|
gapInfo,
|
|
gapMode
|
|
);
|
|
|
|
const point& cc = neiCc[bFaceI];
|
|
label nRays = generateRays
|
|
(
|
|
false,
|
|
surfPt,
|
|
ccNormal1[i],
|
|
gapInfo,
|
|
gapMode,
|
|
surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
|
|
neiLevel[bFaceI],
|
|
|
|
rayStart,
|
|
rayEnd,
|
|
gapSize,
|
|
|
|
rayStart2,
|
|
rayEnd2,
|
|
gapSize2
|
|
);
|
|
for (label j = 0; j < nRays; j++)
|
|
{
|
|
cellMap.append(-1); // See above.
|
|
compactMap.append(i);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Shooting " << returnReduce(rayStart.size(), sumOp<label>())
|
|
<< " rays from " << returnReduce(testFaces.size(), sumOp<label>())
|
|
<< " intersected faces" << endl;
|
|
|
|
rayStart.shrink();
|
|
rayEnd.shrink();
|
|
gapSize.shrink();
|
|
|
|
rayStart2.shrink();
|
|
rayEnd2.shrink();
|
|
gapSize2.shrink();
|
|
|
|
cellMap.shrink();
|
|
compactMap.shrink();
|
|
|
|
testFaces.clear();
|
|
ccSurface1.clear();
|
|
ccHit1.clear();
|
|
ccRegion1.clear();
|
|
ccNormal1 = UIndirectList<vector>(ccNormal1, compactMap)();
|
|
|
|
|
|
// Do intersections in pairs
|
|
labelList surf1;
|
|
List<pointIndexHit> hit1;
|
|
vectorField normal1;
|
|
surfaces_.findNearestIntersection
|
|
(
|
|
rayStart,
|
|
rayEnd,
|
|
surf1,
|
|
hit1,
|
|
normal1
|
|
);
|
|
|
|
labelList surf2;
|
|
List<pointIndexHit> hit2;
|
|
vectorField normal2;
|
|
surfaces_.findNearestIntersection
|
|
(
|
|
rayStart2,
|
|
rayEnd2,
|
|
surf2,
|
|
hit2,
|
|
normal2
|
|
);
|
|
|
|
forAll(surf1, i)
|
|
{
|
|
// Combine selfProx of shell and surfaces.
|
|
// Ignore regions for now
|
|
const label cellI = cellMap[i];
|
|
|
|
const label shelli =
|
|
(
|
|
(cellI != -1 && cellToCompact[cellI] != -1)
|
|
? gapShell[cellToCompact[cellI]]
|
|
: -1
|
|
);
|
|
|
|
bool selfProx = true;
|
|
if (shelli != -1)
|
|
{
|
|
selfProx = shells_.gapSelf()[shelli][0];
|
|
}
|
|
if (surf1[i] != -1 && selfProx)
|
|
{
|
|
const label globalRegioni = surfaces_.globalRegion(surf1[i], 0);
|
|
selfProx = extendedGapSelf[globalRegioni];
|
|
}
|
|
|
|
if
|
|
(
|
|
surf1[i] != -1
|
|
&& surf2[i] != -1
|
|
&& (surf2[i] != surf1[i] || selfProx)
|
|
)
|
|
{
|
|
// Found intersection with surface. Check opposite normal.
|
|
if
|
|
(
|
|
cellI != -1
|
|
&& (mag(normal1[i]&normal2[i]) > planarCos)
|
|
&& (
|
|
magSqr(hit1[i].hitPoint()-hit2[i].hitPoint())
|
|
< Foam::sqr(gapSize[i])
|
|
)
|
|
)
|
|
{
|
|
if
|
|
(
|
|
!markForRefine
|
|
(
|
|
surf1[i],
|
|
nAllowRefine,
|
|
refineCell[cellI],
|
|
nRefine
|
|
)
|
|
)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if
|
|
(
|
|
returnReduce(nRefine, sumOp<label>())
|
|
> returnReduce(nAllowRefine, sumOp<label>())
|
|
)
|
|
{
|
|
Info<< "Reached refinement limit." << endl;
|
|
}
|
|
}
|
|
|
|
return returnReduce(nRefine-oldNRefine, sumOp<label>());
|
|
}
|
|
|
|
|
|
//Foam::meshRefinement::findNearestOppositeOp::findNearestOppositeOp
|
|
//(
|
|
// const indexedOctree<treeDataTriSurface>& tree,
|
|
// const point& oppositePoint,
|
|
// const vector& oppositeNormal,
|
|
// const scalar minCos
|
|
//)
|
|
//:
|
|
// tree_(tree),
|
|
// oppositePoint_(oppositePoint),
|
|
// oppositeNormal_(oppositeNormal),
|
|
// minCos_(minCos)
|
|
//{}
|
|
//
|
|
//
|
|
//void Foam::meshRefinement::findNearestOppositeOp::operator()
|
|
//(
|
|
// const labelUList& indices,
|
|
// const point& sample,
|
|
// scalar& nearestDistSqr,
|
|
// label& minIndex,
|
|
// point& nearestPoint
|
|
//) const
|
|
//{
|
|
// const treeDataTriSurface& shape = tree_.shapes();
|
|
// const triSurface& patch = shape.patch();
|
|
// const pointField& points = patch.points();
|
|
//
|
|
// forAll(indices, i)
|
|
// {
|
|
// const label index = indices[i];
|
|
// const labelledTri& f = patch[index];
|
|
//
|
|
// pointHit nearHit = f.nearestPoint(sample, points);
|
|
// scalar distSqr = sqr(nearHit.distance());
|
|
//
|
|
// if (distSqr < nearestDistSqr)
|
|
// {
|
|
// // Nearer. Check if
|
|
// // - a bit way from other hit
|
|
// // - in correct search cone
|
|
// vector d(nearHit.rawPoint()-oppositePoint_);
|
|
// scalar normalDist(d&oppositeNormal_);
|
|
//
|
|
// if (normalDist > Foam::sqr(SMALL) && normalDist/mag(d) > minCos_)
|
|
// {
|
|
// nearestDistSqr = distSqr;
|
|
// minIndex = index;
|
|
// nearestPoint = nearHit.rawPoint();
|
|
// }
|
|
// }
|
|
// }
|
|
//}
|
|
//
|
|
//
|
|
//void Foam::meshRefinement::searchCone
|
|
//(
|
|
// const label surfI,
|
|
// labelList& nearMap, // cells
|
|
// scalarField& nearGap, // gap size
|
|
// List<pointIndexHit>& nearInfo, // nearest point on surface
|
|
// List<pointIndexHit>& oppositeInfo // detected point on gap (or miss)
|
|
//) const
|
|
//{
|
|
// const labelList& cellLevel = meshCutter_.cellLevel();
|
|
// const pointField& cellCentres = mesh_.cellCentres();
|
|
// const scalar edge0Len = meshCutter_.level0EdgeLength();
|
|
//
|
|
// const labelList& surfaceIndices = surfaces_.surfaces();
|
|
// const List<FixedList<label, 3>>& extendedGapLevel =
|
|
// surfaces_.extendedGapLevel();
|
|
// const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
|
|
//
|
|
//
|
|
// label geomI = surfaceIndices[surfI];
|
|
// const searchableSurface& geom = surfaces_.geometry()[geomI];
|
|
//
|
|
// const triSurfaceMesh& s = refCast<const triSurfaceMesh>(geom);
|
|
// const indexedOctree<treeDataTriSurface>& tree = s.tree();
|
|
//
|
|
//
|
|
// const scalar searchCos = Foam::cos(degToRad(30.0));
|
|
//
|
|
// // Normals for ray shooting and inside/outside detection
|
|
// vectorField nearNormal;
|
|
// geom.getNormal(nearInfo, nearNormal);
|
|
// // Regions
|
|
// labelList nearRegion;
|
|
// geom.getRegion(nearInfo, nearRegion);
|
|
//
|
|
//
|
|
// // Now loop over all near points and search in the half cone
|
|
// labelList map(nearInfo.size());
|
|
// label compactI = 0;
|
|
//
|
|
// oppositeInfo.setSize(nearInfo.size());
|
|
//
|
|
// forAll(nearInfo, i)
|
|
// {
|
|
// label globalRegionI =
|
|
// surfaces_.globalRegion(surfI, nearRegion[i]);
|
|
//
|
|
// // Get updated gap information now we have the region
|
|
// label nGapCells = extendedGapLevel[globalRegionI][0];
|
|
// label minLevel = extendedGapLevel[globalRegionI][1];
|
|
// label maxLevel = extendedGapLevel[globalRegionI][2];
|
|
// volumeType mode = extendedGapMode[globalRegionI];
|
|
//
|
|
// label cellI = nearMap[i];
|
|
// label cLevel = cellLevel[cellI];
|
|
//
|
|
// if (cLevel >= minLevel && cLevel < maxLevel)
|
|
// {
|
|
// scalar cellSize = edge0Len/pow(2.0, cLevel);
|
|
//
|
|
// // Update gap size
|
|
// nearGap[i] = nGapCells*cellSize;
|
|
//
|
|
// const point& nearPt = nearInfo[i].hitPoint();
|
|
// vector v(cellCentres[cellI]-nearPt);
|
|
// scalar magV = mag(v);
|
|
//
|
|
// // Like with ray shooting we want to
|
|
// // - find triangles up to nearGap away on the wanted side of the
|
|
// // surface
|
|
// // - find triangles up to 0.5*cellSize away on the unwanted side
|
|
// // of the surface. This is for cells straddling the surface
|
|
// // where
|
|
// // the cell centre might be on the wrong side of the surface
|
|
//
|
|
// // Tbd: check that cell centre is inbetween the gap hits
|
|
// // (only if the cell is far enough away)
|
|
//
|
|
// scalar posNormalSize = 0.0;
|
|
// scalar negNormalSize = 0.0;
|
|
//
|
|
// if (mode == volumeType::OUTSIDE)
|
|
// {
|
|
// posNormalSize = nearGap[i];
|
|
// if (magV < 0.5*cellSize)
|
|
// {
|
|
// negNormalSize = 0.5*cellSize;
|
|
// }
|
|
// }
|
|
// else if (mode == volumeType::INSIDE)
|
|
// {
|
|
// if (magV < 0.5*cellSize)
|
|
// {
|
|
// posNormalSize = 0.5*cellSize;
|
|
// }
|
|
// negNormalSize = nearGap[i];
|
|
// }
|
|
// else
|
|
// {
|
|
// posNormalSize = nearGap[i];
|
|
// negNormalSize = nearGap[i];
|
|
// }
|
|
//
|
|
// // Test with positive normal
|
|
// oppositeInfo[compactI] = tree.findNearest
|
|
// (
|
|
// nearPt,
|
|
// sqr(posNormalSize),
|
|
// findNearestOppositeOp
|
|
// (
|
|
// tree,
|
|
// nearPt,
|
|
// nearNormal[i],
|
|
// searchCos
|
|
// )
|
|
// );
|
|
//
|
|
// if (oppositeInfo[compactI].hit())
|
|
// {
|
|
// map[compactI++] = i;
|
|
// }
|
|
// else
|
|
// {
|
|
// // Test with negative normal
|
|
// oppositeInfo[compactI] = tree.findNearest
|
|
// (
|
|
// nearPt,
|
|
// sqr(negNormalSize),
|
|
// findNearestOppositeOp
|
|
// (
|
|
// tree,
|
|
// nearPt,
|
|
// -nearNormal[i],
|
|
// searchCos
|
|
// )
|
|
// );
|
|
//
|
|
// if (oppositeInfo[compactI].hit())
|
|
// {
|
|
// map[compactI++] = i;
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// Info<< "Selected " << returnReduce(compactI, sumOp<label>())
|
|
// << " hits on the correct side out of "
|
|
// << returnReduce(map.size(), sumOp<label>()) << endl;
|
|
// map.setSize(compactI);
|
|
// oppositeInfo.setSize(compactI);
|
|
//
|
|
// nearMap = labelUIndList(nearMap, map)();
|
|
// nearGap = UIndirectList<scalar>(nearGap, map)();
|
|
// nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
|
|
// nearNormal = UIndirectList<vector>(nearNormal, map)();
|
|
//
|
|
// // Exclude hits which aren't opposite enough. E.g. you might find
|
|
// // a point on a perpendicular wall - but this does not constitute a gap.
|
|
// vectorField oppositeNormal;
|
|
// geom.getNormal(oppositeInfo, oppositeNormal);
|
|
//
|
|
// compactI = 0;
|
|
// forAll(oppositeInfo, i)
|
|
// {
|
|
// if ((nearNormal[i] & oppositeNormal[i]) < -0.707)
|
|
// {
|
|
// map[compactI++] = i;
|
|
// }
|
|
// }
|
|
//
|
|
// Info<< "Selected " << returnReduce(compactI, sumOp<label>())
|
|
// << " hits opposite the nearest out of "
|
|
// << returnReduce(map.size(), sumOp<label>()) << endl;
|
|
// map.setSize(compactI);
|
|
//
|
|
// nearMap = labelUIndList(nearMap, map)();
|
|
// nearGap = UIndirectList<scalar>(nearGap, map)();
|
|
// nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
|
|
// oppositeInfo = UIndirectList<pointIndexHit>(oppositeInfo, map)();
|
|
//}
|
|
|
|
|
|
Foam::label Foam::meshRefinement::generateRays
|
|
(
|
|
const point& nearPoint,
|
|
const vector& nearNormal,
|
|
const FixedList<label, 3>& gapInfo,
|
|
const volumeType& mode,
|
|
|
|
const label cLevel,
|
|
|
|
DynamicField<point>& start,
|
|
DynamicField<point>& end
|
|
) const
|
|
{
|
|
label nOldRays = start.size();
|
|
|
|
if (cLevel >= gapInfo[1] && cLevel < gapInfo[2] && gapInfo[0] > 0)
|
|
{
|
|
scalar cellSize = meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
|
|
|
|
// Calculate gap size
|
|
scalar nearGap = gapInfo[0]*cellSize;
|
|
|
|
const vector& n = nearNormal;
|
|
|
|
// Situation 'C' above: cell too close. Use surface
|
|
// -normal and -point to shoot rays
|
|
|
|
if (mode == volumeType::OUTSIDE)
|
|
{
|
|
start.append(nearPoint+1e-6*n);
|
|
end.append(nearPoint+nearGap*n);
|
|
}
|
|
else if (mode == volumeType::INSIDE)
|
|
{
|
|
start.append(nearPoint-1e-6*n);
|
|
end.append(nearPoint-nearGap*n);
|
|
}
|
|
else if (mode == volumeType::MIXED)
|
|
{
|
|
start.append(nearPoint+1e-6*n);
|
|
end.append(nearPoint+nearGap*n);
|
|
|
|
start.append(nearPoint-1e-6*n);
|
|
end.append(nearPoint-nearGap*n);
|
|
}
|
|
}
|
|
|
|
return start.size()-nOldRays;
|
|
}
|
|
|
|
|
|
Foam::label Foam::meshRefinement::generateRays
|
|
(
|
|
const bool useSurfaceNormal,
|
|
|
|
const point& nearPoint,
|
|
const vector& nearNormal,
|
|
const FixedList<label, 3>& gapInfo,
|
|
const volumeType& mode,
|
|
|
|
const point& cc,
|
|
const label cLevel,
|
|
|
|
DynamicField<point>& start,
|
|
DynamicField<point>& end,
|
|
DynamicField<scalar>& gapSize,
|
|
|
|
DynamicField<point>& start2,
|
|
DynamicField<point>& end2,
|
|
DynamicField<scalar>& gapSize2
|
|
) const
|
|
{
|
|
// We want to handle the following cases:
|
|
// - surface: small gap (marked with 'surface'). gap might be
|
|
// on inside or outside of surface.
|
|
// - A: cell well inside the gap.
|
|
// - B: cell well outside the gap.
|
|
// - C: cell straddling the gap. cell centre might be inside
|
|
// or outside
|
|
//
|
|
// +---+
|
|
// | B |
|
|
// +---+
|
|
//
|
|
// +------+
|
|
// | |
|
|
// | C |
|
|
// --------|------|----surface
|
|
// +------+
|
|
//
|
|
// +---+
|
|
// | A |
|
|
// +---+
|
|
//
|
|
//
|
|
// --------------------surface
|
|
//
|
|
// So:
|
|
// - find nearest point on surface
|
|
// - in situation A,B decide if on wanted side of surface
|
|
// - detect if locally a gap (and the cell inside the gap) by
|
|
// shooting a ray from the point on the surface in the direction
|
|
// of
|
|
// - A,B: the cell centre
|
|
// - C: the surface normal and/or negative surface normal
|
|
// and see we hit anything
|
|
//
|
|
// Variations of this scheme:
|
|
// - always shoot in the direction of the surface normal. This needs
|
|
// then an additional check to make sure the cell centre is
|
|
// somewhere inside the gap
|
|
// - instead of ray shooting use a 'constrained' nearest search
|
|
// by e.g. looking inside a search cone (implemented in searchCone).
|
|
// The problem with this constrained nearest is that it still uses
|
|
// the absolute nearest point on each triangle and only afterwards
|
|
// checks if it is inside the search cone.
|
|
|
|
|
|
// Decide which near points are good:
|
|
// - with updated minLevel and maxLevel and nearGap make sure
|
|
// the cell is still a candidate
|
|
// NOTE: inside the gap the nearest point on the surface will
|
|
// be HALF the gap size - otherwise we would have found
|
|
// a point on the opposite side
|
|
// - if the mode is both sides
|
|
// - or if the hit is inside the current cell (situation 'C',
|
|
// magV < 0.5cellSize)
|
|
// - or otherwise if on the correct side
|
|
|
|
label nOldRays = start.size();
|
|
|
|
if (cLevel >= gapInfo[1] && cLevel < gapInfo[2] && gapInfo[0] > 0)
|
|
{
|
|
scalar cellSize = meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
|
|
|
|
// Calculate gap size
|
|
scalar nearGap = gapInfo[0]*cellSize;
|
|
|
|
// Distance to nearest
|
|
vector v(cc-nearPoint);
|
|
scalar magV = mag(v);
|
|
|
|
if (useSurfaceNormal || magV < 0.5*cellSize)
|
|
{
|
|
const vector& n = nearNormal;
|
|
|
|
// Situation 'C' above: cell too close. Use surface
|
|
// -normal and -point to shoot rays
|
|
|
|
if (mode == volumeType::OUTSIDE)
|
|
{
|
|
start.append(nearPoint+1e-6*n);
|
|
end.append(nearPoint+nearGap*n);
|
|
gapSize.append(nearGap);
|
|
// Second vector so we get pairs of intersections
|
|
start2.append(nearPoint+1e-6*n);
|
|
end2.append(nearPoint-1e-6*n);
|
|
gapSize2.append(gapSize.last());
|
|
}
|
|
else if (mode == volumeType::INSIDE)
|
|
{
|
|
start.append(nearPoint-1e-6*n);
|
|
end.append(nearPoint-nearGap*n);
|
|
gapSize.append(nearGap);
|
|
// Second vector so we get pairs of intersections
|
|
start2.append(nearPoint-1e-6*n);
|
|
end2.append(nearPoint+1e-6*n);
|
|
gapSize2.append(gapSize.last());
|
|
}
|
|
else if (mode == volumeType::MIXED)
|
|
{
|
|
// Do both rays:
|
|
// Outside
|
|
{
|
|
start.append(nearPoint+1e-6*n);
|
|
end.append(nearPoint+nearGap*n);
|
|
gapSize.append(nearGap);
|
|
// Second vector so we get pairs of intersections
|
|
start2.append(nearPoint+1e-6*n);
|
|
end2.append(nearPoint-1e-6*n);
|
|
gapSize2.append(gapSize.last());
|
|
}
|
|
// Inside
|
|
{
|
|
start.append(nearPoint-1e-6*n);
|
|
end.append(nearPoint-nearGap*n);
|
|
gapSize.append(nearGap);
|
|
// Second vector so we get pairs of intersections
|
|
start2.append(nearPoint-1e-6*n);
|
|
end2.append(nearPoint+1e-6*n);
|
|
gapSize2.append(gapSize.last());
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Situation 'A' or 'B' above: cell well away. Test if
|
|
// cell on correct side of surface and shoot ray through
|
|
// cell centre. Note: no need to shoot ray in other
|
|
// direction since we're trying to detect cell inside
|
|
// the gap.
|
|
|
|
scalar s = (v&nearNormal);
|
|
|
|
if
|
|
(
|
|
(mode == volumeType::MIXED)
|
|
|| (mode == volumeType::OUTSIDE && s > SMALL)
|
|
|| (mode == volumeType::INSIDE && s < -SMALL)
|
|
)
|
|
{
|
|
//// Use single vector through cell centre
|
|
//vector n(v/(magV+ROOTVSMALL));
|
|
//
|
|
//start.append(cc);
|
|
//end.append(cc+nearGap*n);
|
|
//gapSize.append(nearGap);
|
|
//
|
|
//start2.append(cc);
|
|
//end2.append(cc-nearGap*n);
|
|
//gapSize2.append(nearGap);
|
|
|
|
|
|
//// Shoot some rays through the cell centre
|
|
//// X-direction:
|
|
//start.append(cc);
|
|
//end.append(cc+nearGap*vector(1, 0, 0));
|
|
//gapSize.append(nearGap);
|
|
//
|
|
//start2.append(cc);
|
|
//end2.append(cc-nearGap*vector(1, 0, 0));
|
|
//gapSize2.append(nearGap);
|
|
//
|
|
//// Y-direction:
|
|
//start.append(cc);
|
|
//end.append(cc+nearGap*vector(0, 1, 0));
|
|
//gapSize.append(nearGap);
|
|
//
|
|
//start2.append(cc);
|
|
//end2.append(cc-nearGap*vector(0, 1, 0));
|
|
//gapSize2.append(nearGap);
|
|
//
|
|
//// Z-direction:
|
|
//start.append(cc);
|
|
//end.append(cc+nearGap*vector(0, 0, 1));
|
|
//gapSize.append(nearGap);
|
|
//
|
|
//start2.append(cc);
|
|
//end2.append(cc-nearGap*vector(0, 0, 1));
|
|
//gapSize2.append(nearGap);
|
|
|
|
|
|
// 3 axes aligned with normal
|
|
|
|
// Use vector through cell centre
|
|
vector n(v/(magV+ROOTVSMALL));
|
|
|
|
// Get second vector. Make sure it is sufficiently perpendicular
|
|
vector e2(1, 0, 0);
|
|
scalar s = (e2 & n);
|
|
if (mag(s) < 0.9)
|
|
{
|
|
e2 -= s*n;
|
|
}
|
|
else
|
|
{
|
|
e2 = vector(0, 1, 0);
|
|
e2 -= (e2 & n)*n;
|
|
}
|
|
e2 /= mag(e2);
|
|
|
|
// Third vector
|
|
vector e3 = n ^ e2;
|
|
|
|
|
|
// Rays in first direction
|
|
start.append(cc);
|
|
end.append(cc+nearGap*n);
|
|
gapSize.append(nearGap);
|
|
|
|
start2.append(cc);
|
|
end2.append(cc-nearGap*n);
|
|
gapSize2.append(nearGap);
|
|
|
|
// Rays in second direction
|
|
start.append(cc);
|
|
end.append(cc+nearGap*e2);
|
|
gapSize.append(nearGap);
|
|
|
|
start2.append(cc);
|
|
end2.append(cc-nearGap*e2);
|
|
gapSize2.append(nearGap);
|
|
|
|
// Rays in third direction
|
|
start.append(cc);
|
|
end.append(cc+nearGap*e3);
|
|
gapSize.append(nearGap);
|
|
|
|
start2.append(cc);
|
|
end2.append(cc-nearGap*e3);
|
|
gapSize2.append(nearGap);
|
|
}
|
|
}
|
|
}
|
|
|
|
return start.size()-nOldRays;
|
|
}
|
|
|
|
|
|
void Foam::meshRefinement::selectGapCandidates
|
|
(
|
|
const labelList& refineCell,
|
|
const label nRefine,
|
|
|
|
labelList& cellMap,
|
|
labelList& gapShell,
|
|
List<FixedList<label, 3>>& shellGapInfo,
|
|
List<volumeType>& shellGapMode
|
|
) const
|
|
{
|
|
const labelList& cellLevel = meshCutter_.cellLevel();
|
|
const pointField& cellCentres = mesh_.cellCentres();
|
|
|
|
// Collect cells to test
|
|
cellMap.setSize(cellLevel.size()-nRefine);
|
|
label compactI = 0;
|
|
|
|
forAll(cellLevel, cellI)
|
|
{
|
|
if (refineCell[cellI] == -1)
|
|
{
|
|
cellMap[compactI++] = cellI;
|
|
}
|
|
}
|
|
Info<< "Selected " << returnReduce(compactI, sumOp<label>())
|
|
<< " unmarked cells out of "
|
|
<< mesh_.globalData().nTotalCells() << endl;
|
|
cellMap.setSize(compactI);
|
|
|
|
// Do test to see whether cells are inside/outside shell with
|
|
// applicable specification (minLevel <= celllevel < maxLevel)
|
|
shells_.findHigherGapLevel
|
|
(
|
|
pointField(cellCentres, cellMap),
|
|
labelUIndList(cellLevel, cellMap)(),
|
|
|
|
gapShell,
|
|
shellGapInfo,
|
|
shellGapMode
|
|
);
|
|
|
|
// Compact out hits
|
|
|
|
labelList map(shellGapInfo.size());
|
|
compactI = 0;
|
|
forAll(shellGapInfo, i)
|
|
{
|
|
if (shellGapInfo[i][2] > 0)
|
|
{
|
|
map[compactI++] = i;
|
|
}
|
|
}
|
|
|
|
Info<< "Selected " << returnReduce(compactI, sumOp<label>())
|
|
<< " cells inside gap shells out of "
|
|
<< mesh_.globalData().nTotalCells() << endl;
|
|
|
|
map.setSize(compactI);
|
|
cellMap = labelUIndList(cellMap, map)();
|
|
gapShell = labelUIndList(gapShell, map)();
|
|
shellGapInfo = UIndirectList<FixedList<label, 3>>(shellGapInfo, map)();
|
|
shellGapMode = UIndirectList<volumeType>(shellGapMode, map)();
|
|
}
|
|
|
|
|
|
void Foam::meshRefinement::mergeGapInfo
|
|
(
|
|
const FixedList<label, 3>& shellGapInfo,
|
|
const volumeType shellGapMode,
|
|
const FixedList<label, 3>& surfGapInfo,
|
|
const volumeType surfGapMode,
|
|
|
|
FixedList<label, 3>& gapInfo,
|
|
volumeType& gapMode
|
|
) const
|
|
{
|
|
if (surfGapInfo[0] == 0)
|
|
{
|
|
gapInfo = shellGapInfo;
|
|
gapMode = shellGapMode;
|
|
}
|
|
else if (shellGapInfo[0] == 0)
|
|
{
|
|
gapInfo = surfGapInfo;
|
|
gapMode = surfGapMode;
|
|
}
|
|
else
|
|
{
|
|
// Both specify a level. Does surface level win? Or does information
|
|
// need to be merged?
|
|
|
|
//gapInfo[0] = max(surfGapInfo[0], shellGapInfo[0]);
|
|
//gapInfo[1] = min(surfGapInfo[1], shellGapInfo[1]);
|
|
//gapInfo[2] = max(surfGapInfo[2], shellGapInfo[2]);
|
|
gapInfo = surfGapInfo;
|
|
gapMode = surfGapMode;
|
|
}
|
|
}
|
|
|
|
|
|
Foam::label Foam::meshRefinement::markInternalGapRefinement
|
|
(
|
|
const scalar planarCos,
|
|
const bool spreadGapSize,
|
|
const label nAllowRefine,
|
|
|
|
labelList& refineCell,
|
|
label& nRefine,
|
|
labelList& numGapCells,
|
|
scalarField& detectedGapSize
|
|
) const
|
|
{
|
|
detectedGapSize.setSize(mesh_.nCells());
|
|
detectedGapSize = GREAT;
|
|
numGapCells.setSize(mesh_.nCells());
|
|
numGapCells = -1;
|
|
|
|
const labelList& cellLevel = meshCutter_.cellLevel();
|
|
const pointField& cellCentres = mesh_.cellCentres();
|
|
const scalar edge0Len = meshCutter_.level0EdgeLength();
|
|
|
|
const List<FixedList<label, 3>>& extendedGapLevel =
|
|
surfaces_.extendedGapLevel();
|
|
const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
|
|
const boolList& extendedGapSelf = surfaces_.gapSelf();
|
|
|
|
// Get the gap level for the shells
|
|
const labelList maxLevel(shells_.maxGapLevel());
|
|
|
|
label oldNRefine = nRefine;
|
|
|
|
if (max(maxLevel) > 0)
|
|
{
|
|
// Collect cells to test
|
|
labelList cellMap;
|
|
labelList gapShell;
|
|
List<FixedList<label, 3>> shellGapInfo;
|
|
List<volumeType> shellGapMode;
|
|
selectGapCandidates
|
|
(
|
|
refineCell,
|
|
nRefine,
|
|
|
|
cellMap,
|
|
gapShell,
|
|
shellGapInfo,
|
|
shellGapMode
|
|
);
|
|
|
|
// Find nearest point and normal on the surfaces
|
|
List<pointIndexHit> nearInfo;
|
|
vectorField nearNormal;
|
|
labelList nearSurface;
|
|
labelList nearRegion;
|
|
{
|
|
// Now we have both the cell-level and the gap size information. Use
|
|
// this to calculate the gap size
|
|
scalarField gapSize(cellMap.size());
|
|
forAll(cellMap, i)
|
|
{
|
|
label cellI = cellMap[i];
|
|
scalar cellSize = edge0Len/pow(2.0, cellLevel[cellI]);
|
|
gapSize[i] = shellGapInfo[i][0]*cellSize;
|
|
}
|
|
|
|
surfaces_.findNearestRegion
|
|
(
|
|
identity(surfaces_.surfaces().size()),
|
|
pointField(cellCentres, cellMap),
|
|
sqr(gapSize),
|
|
nearSurface,
|
|
nearInfo,
|
|
nearRegion,
|
|
nearNormal
|
|
);
|
|
}
|
|
|
|
|
|
|
|
DynamicList<label> map(nearInfo.size());
|
|
DynamicField<point> rayStart(nearInfo.size());
|
|
DynamicField<point> rayEnd(nearInfo.size());
|
|
DynamicField<scalar> gapSize(nearInfo.size());
|
|
|
|
DynamicField<point> rayStart2(nearInfo.size());
|
|
DynamicField<point> rayEnd2(nearInfo.size());
|
|
DynamicField<scalar> gapSize2(nearInfo.size());
|
|
|
|
label nTestCells = 0;
|
|
|
|
forAll(nearInfo, i)
|
|
{
|
|
if (nearInfo[i].hit())
|
|
{
|
|
label globalRegionI = surfaces_.globalRegion
|
|
(
|
|
nearSurface[i],
|
|
nearRegion[i]
|
|
);
|
|
|
|
// Combine info from shell and surface
|
|
FixedList<label, 3> gapInfo;
|
|
volumeType gapMode;
|
|
mergeGapInfo
|
|
(
|
|
shellGapInfo[i],
|
|
shellGapMode[i],
|
|
|
|
extendedGapLevel[globalRegionI],
|
|
extendedGapMode[globalRegionI],
|
|
|
|
gapInfo,
|
|
gapMode
|
|
);
|
|
|
|
// Store wanted number of cells in gap
|
|
label cellI = cellMap[i];
|
|
label cLevel = cellLevel[cellI];
|
|
if (cLevel >= gapInfo[1] && cLevel < gapInfo[2])
|
|
{
|
|
numGapCells[cellI] = max(numGapCells[cellI], gapInfo[0]);
|
|
}
|
|
|
|
// Construct one or more rays to test for oppositeness
|
|
label nRays = generateRays
|
|
(
|
|
false,
|
|
nearInfo[i].hitPoint(),
|
|
nearNormal[i],
|
|
gapInfo,
|
|
gapMode,
|
|
|
|
cellCentres[cellI],
|
|
cLevel,
|
|
|
|
rayStart,
|
|
rayEnd,
|
|
gapSize,
|
|
|
|
rayStart2,
|
|
rayEnd2,
|
|
gapSize2
|
|
);
|
|
if (nRays > 0)
|
|
{
|
|
nTestCells++;
|
|
for (label j = 0; j < nRays; j++)
|
|
{
|
|
map.append(i);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
|
|
<< " cells for testing out of "
|
|
<< mesh_.globalData().nTotalCells() << endl;
|
|
map.shrink();
|
|
rayStart.shrink();
|
|
rayEnd.shrink();
|
|
gapSize.shrink();
|
|
|
|
rayStart2.shrink();
|
|
rayEnd2.shrink();
|
|
gapSize2.shrink();
|
|
|
|
cellMap = labelUIndList(cellMap, map)();
|
|
nearNormal = UIndirectList<vector>(nearNormal, map)();
|
|
shellGapInfo.clear();
|
|
shellGapMode.clear();
|
|
nearInfo.clear();
|
|
nearSurface.clear();
|
|
nearRegion.clear();
|
|
|
|
|
|
// Do intersections in pairs
|
|
labelList surf1;
|
|
List<pointIndexHit> hit1;
|
|
vectorField normal1;
|
|
surfaces_.findNearestIntersection
|
|
(
|
|
rayStart,
|
|
rayEnd,
|
|
surf1,
|
|
hit1,
|
|
normal1
|
|
);
|
|
|
|
labelList surf2;
|
|
List<pointIndexHit> hit2;
|
|
vectorField normal2;
|
|
surfaces_.findNearestIntersection
|
|
(
|
|
rayStart2,
|
|
rayEnd2,
|
|
surf2,
|
|
hit2,
|
|
normal2
|
|
);
|
|
|
|
// Extract cell based gap size
|
|
forAll(surf1, i)
|
|
{
|
|
// Combine selfProx of shell and surfaces. Ignore regions for
|
|
// now
|
|
const label shelli = gapShell[map[i]];
|
|
|
|
bool selfProx = true;
|
|
if (shelli != -1)
|
|
{
|
|
selfProx = shells_.gapSelf()[shelli][0];
|
|
}
|
|
if (surf1[i] != -1 && selfProx)
|
|
{
|
|
const label globalRegioni = surfaces_.globalRegion(surf1[i], 0);
|
|
selfProx = extendedGapSelf[globalRegioni];
|
|
}
|
|
|
|
if
|
|
(
|
|
surf1[i] != -1
|
|
&& surf2[i] != -1
|
|
&& (surf2[i] != surf1[i] || selfProx)
|
|
)
|
|
{
|
|
// Found intersections with surface. Check for
|
|
// - small gap
|
|
// - coplanar normals
|
|
|
|
const label cellI = cellMap[i];
|
|
|
|
const scalar d2 = magSqr(hit1[i].hitPoint()-hit2[i].hitPoint());
|
|
|
|
if
|
|
(
|
|
cellI != -1
|
|
&& (mag(normal1[i]&normal2[i]) > planarCos)
|
|
&& (d2 < Foam::sqr(gapSize[i]))
|
|
)
|
|
{
|
|
detectedGapSize[cellI] = min
|
|
(
|
|
detectedGapSize[cellI],
|
|
Foam::sqrt(d2)
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Spread it
|
|
if (spreadGapSize)
|
|
{
|
|
// Field on cells and faces
|
|
List<transportData> cellData(mesh_.nCells());
|
|
List<transportData> faceData(mesh_.nFaces());
|
|
|
|
// Start of walk
|
|
const pointField& faceCentres = mesh_.faceCentres();
|
|
|
|
DynamicList<label> frontFaces(mesh_.nFaces());
|
|
DynamicList<transportData> frontData(mesh_.nFaces());
|
|
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
|
|
{
|
|
label own = mesh_.faceOwner()[faceI];
|
|
label nei = mesh_.faceNeighbour()[faceI];
|
|
|
|
scalar minSize = min
|
|
(
|
|
detectedGapSize[own],
|
|
detectedGapSize[nei]
|
|
);
|
|
|
|
if (minSize < GREAT)
|
|
{
|
|
frontFaces.append(faceI);
|
|
frontData.append
|
|
(
|
|
transportData
|
|
(
|
|
faceCentres[faceI],
|
|
minSize,
|
|
0.0
|
|
)
|
|
);
|
|
}
|
|
}
|
|
for
|
|
(
|
|
label faceI = mesh_.nInternalFaces();
|
|
faceI < mesh_.nFaces();
|
|
faceI++
|
|
)
|
|
{
|
|
label own = mesh_.faceOwner()[faceI];
|
|
|
|
if (detectedGapSize[own] < GREAT)
|
|
{
|
|
frontFaces.append(faceI);
|
|
frontData.append
|
|
(
|
|
transportData
|
|
(
|
|
faceCentres[faceI],
|
|
detectedGapSize[own],
|
|
0.0
|
|
)
|
|
);
|
|
}
|
|
}
|
|
|
|
Info<< "Selected "
|
|
<< returnReduce(frontFaces.size(), sumOp<label>())
|
|
<< " faces for spreading gap size out of "
|
|
<< mesh_.globalData().nTotalFaces() << endl;
|
|
|
|
|
|
transportData::trackData td(surfaceIndex());
|
|
|
|
FaceCellWave<transportData, transportData::trackData> deltaCalc
|
|
(
|
|
mesh_,
|
|
frontFaces,
|
|
frontData,
|
|
faceData,
|
|
cellData,
|
|
mesh_.globalData().nTotalCells()+1,
|
|
td
|
|
);
|
|
|
|
|
|
forAll(cellMap, i)
|
|
{
|
|
label cellI = cellMap[i];
|
|
if
|
|
(
|
|
cellI != -1
|
|
&& cellData[cellI].valid(deltaCalc.data())
|
|
&& numGapCells[cellI] != -1
|
|
)
|
|
{
|
|
// Update transported gap size
|
|
detectedGapSize[cellI] = min
|
|
(
|
|
detectedGapSize[cellI],
|
|
cellData[cellI].data()
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Use it
|
|
forAll(cellMap, i)
|
|
{
|
|
label cellI = cellMap[i];
|
|
|
|
if (cellI != -1 && numGapCells[cellI] != -1)
|
|
{
|
|
// Needed gap size
|
|
label cLevel = cellLevel[cellI];
|
|
scalar cellSize =
|
|
meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
|
|
scalar neededGapSize = numGapCells[cellI]*cellSize;
|
|
|
|
if (neededGapSize > detectedGapSize[cellI])
|
|
{
|
|
if
|
|
(
|
|
!markForRefine
|
|
(
|
|
123,
|
|
nAllowRefine,
|
|
refineCell[cellI],
|
|
nRefine
|
|
)
|
|
)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
if
|
|
(
|
|
returnReduce(nRefine, sumOp<label>())
|
|
> returnReduce(nAllowRefine, sumOp<label>())
|
|
)
|
|
{
|
|
Info<< "Reached refinement limit." << endl;
|
|
}
|
|
}
|
|
|
|
return returnReduce(nRefine-oldNRefine, sumOp<label>());
|
|
}
|
|
|
|
|
|
Foam::label Foam::meshRefinement::markSmallFeatureRefinement
|
|
(
|
|
const scalar planarCos,
|
|
const label nAllowRefine,
|
|
const labelList& neiLevel,
|
|
const pointField& neiCc,
|
|
|
|
labelList& refineCell,
|
|
label& nRefine
|
|
) const
|
|
{
|
|
const labelList& cellLevel = meshCutter_.cellLevel();
|
|
const labelList& surfaceIndices = surfaces_.surfaces();
|
|
const List<FixedList<label, 3>>& extendedGapLevel =
|
|
surfaces_.extendedGapLevel();
|
|
const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
|
|
const boolList& extendedGapSelf = surfaces_.gapSelf();
|
|
|
|
label oldNRefine = nRefine;
|
|
|
|
// Check that we're using any gap refinement
|
|
labelList shellMaxLevel(shells_.maxGapLevel());
|
|
|
|
if (max(shellMaxLevel) == 0)
|
|
{
|
|
return 0;
|
|
}
|
|
|
|
//- Force calculation of tetBasePt
|
|
(void)mesh_.tetBasePtIs();
|
|
(void)mesh_.cellTree();
|
|
|
|
|
|
forAll(surfaceIndices, surfI)
|
|
{
|
|
label geomI = surfaceIndices[surfI];
|
|
const searchableSurface& geom = surfaces_.geometry()[geomI];
|
|
|
|
|
|
// Get the element index in a roundabout way. Problem is e.g.
|
|
// distributed surface where local indices differ from global
|
|
// ones (needed for getRegion call)
|
|
|
|
pointField ctrs;
|
|
labelList region;
|
|
vectorField normal;
|
|
{
|
|
// Representative local coordinates and bounding sphere
|
|
scalarField radiusSqr;
|
|
geom.boundingSpheres(ctrs, radiusSqr);
|
|
|
|
List<pointIndexHit> info;
|
|
geom.findNearest(ctrs, radiusSqr, info);
|
|
|
|
forAll(info, i)
|
|
{
|
|
if (!info[i].hit())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "fc:" << ctrs[i]
|
|
<< " radius:" << radiusSqr[i]
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
|
|
geom.getRegion(info, region);
|
|
geom.getNormal(info, normal);
|
|
}
|
|
|
|
// Do test to see whether triangles are inside/outside shell with
|
|
// applicable specification (minLevel <= celllevel < maxLevel)
|
|
List<FixedList<label, 3>> shellGapInfo;
|
|
List<volumeType> shellGapMode;
|
|
labelList gapShell;
|
|
shells_.findHigherGapLevel
|
|
(
|
|
ctrs,
|
|
labelList(ctrs.size(), Zero),
|
|
|
|
gapShell,
|
|
shellGapInfo,
|
|
shellGapMode
|
|
);
|
|
|
|
|
|
DynamicList<label> map(ctrs.size());
|
|
DynamicList<label> cellMap(ctrs.size());
|
|
|
|
DynamicField<point> rayStart(ctrs.size());
|
|
DynamicField<point> rayEnd(ctrs.size());
|
|
DynamicField<scalar> gapSize(ctrs.size());
|
|
|
|
label nTestCells = 0;
|
|
|
|
forAll(ctrs, i)
|
|
{
|
|
if (shellGapInfo[i][2] > 0)
|
|
{
|
|
label globalRegionI = surfaces_.globalRegion(surfI, region[i]);
|
|
|
|
// Combine info from shell and surface
|
|
FixedList<label, 3> gapInfo;
|
|
volumeType gapMode;
|
|
mergeGapInfo
|
|
(
|
|
shellGapInfo[i],
|
|
shellGapMode[i],
|
|
|
|
extendedGapLevel[globalRegionI],
|
|
extendedGapMode[globalRegionI],
|
|
|
|
gapInfo,
|
|
gapMode
|
|
);
|
|
|
|
//- Option 1: use octree nearest searching inside polyMesh
|
|
//label cellI = mesh_.findCell(pt, polyMesh::CELL_TETS);
|
|
|
|
//- Option 2: use octree 'inside' searching inside polyMesh. Is
|
|
// much faster.
|
|
label cellI = -1;
|
|
const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
|
|
if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
|
|
{
|
|
cellI = tree.findInside(ctrs[i]);
|
|
}
|
|
|
|
if (cellI != -1 && refineCell[cellI] == -1)
|
|
{
|
|
// Construct one or two rays to test for oppositeness
|
|
// Note that we always want to use the surface normal
|
|
// and not the vector from cell centre to surface point
|
|
|
|
label nRays = generateRays
|
|
(
|
|
ctrs[i],
|
|
normal[i],
|
|
gapInfo,
|
|
gapMode,
|
|
|
|
cellLevel[cellI],
|
|
|
|
rayStart,
|
|
rayEnd
|
|
);
|
|
|
|
if (nRays > 0)
|
|
{
|
|
nTestCells++;
|
|
for (label j = 0; j < nRays; j++)
|
|
{
|
|
cellMap.append(cellI);
|
|
map.append(i);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
|
|
<< " cells containing triangle centres out of "
|
|
<< mesh_.globalData().nTotalCells() << endl;
|
|
map.shrink();
|
|
cellMap.shrink();
|
|
rayStart.shrink();
|
|
rayEnd.shrink();
|
|
|
|
ctrs.clear();
|
|
region.clear();
|
|
shellGapInfo.clear();
|
|
shellGapMode.clear();
|
|
normal = UIndirectList<vector>(normal, map)();
|
|
|
|
// Do intersections.
|
|
labelList surfaceHit;
|
|
vectorField surfaceNormal;
|
|
surfaces_.findNearestIntersection
|
|
(
|
|
rayStart,
|
|
rayEnd,
|
|
surfaceHit,
|
|
surfaceNormal
|
|
);
|
|
|
|
|
|
label nOldRefine = 0;
|
|
|
|
|
|
forAll(surfaceHit, i)
|
|
{
|
|
// Combine selfProx of shell and surfaces. Ignore regions for
|
|
// now
|
|
const label shelli = gapShell[map[i]];
|
|
bool selfProx = true;
|
|
if (shelli != -1)
|
|
{
|
|
selfProx = shells_.gapSelf()[shelli][0];
|
|
}
|
|
if (surfI != -1 && selfProx)
|
|
{
|
|
const label globalRegioni = surfaces_.globalRegion(surfI, 0);
|
|
selfProx = extendedGapSelf[globalRegioni];
|
|
}
|
|
|
|
if
|
|
(
|
|
surfaceHit[i] != -1
|
|
&& (surfaceHit[i] != surfI || selfProx)
|
|
)
|
|
{
|
|
// Found intersection with surface. Check coplanar normals.
|
|
label cellI = cellMap[i];
|
|
|
|
if (mag(normal[i]&surfaceNormal[i]) > planarCos)
|
|
{
|
|
if
|
|
(
|
|
!markForRefine
|
|
(
|
|
surfaceHit[i],
|
|
nAllowRefine,
|
|
refineCell[cellI],
|
|
nRefine
|
|
)
|
|
)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "For surface " << geom.name() << " found "
|
|
<< returnReduce(nRefine-nOldRefine, sumOp<label>())
|
|
<< " cells in small gaps" << endl;
|
|
|
|
if
|
|
(
|
|
returnReduce(nRefine, sumOp<label>())
|
|
> returnReduce(nAllowRefine, sumOp<label>())
|
|
)
|
|
{
|
|
Info<< "Reached refinement limit." << endl;
|
|
}
|
|
}
|
|
|
|
return returnReduce(nRefine-oldNRefine, sumOp<label>());
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|