mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
3bc8ab2839
- return internalEdges() and boundaryEdges() sub lists directly - calculate and return boundaryFaces() to identify faces attached to boundary edges. - minor code cleanup, and add PrimitivePatchBase class for non-templated code. STYLE: mark unused parameter in globalMeshData mergePoints
501 lines
15 KiB
C
501 lines
15 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2011-2016 OpenFOAM Foundation
|
|
Copyright (C) 2020-2021 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/>.
|
|
|
|
Description
|
|
For every point on the patch find the closest face on the target side.
|
|
Return a target face label for each patch point
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "boolList.H"
|
|
#include "PointHit.H"
|
|
#include "objectHit.H"
|
|
#include "bandCompression.H"
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
template<class FaceList, class PointField>
|
|
template<class ToPatch>
|
|
Foam::List<Foam::objectHit>
|
|
Foam::PrimitivePatch<FaceList, PointField>::projectPoints
|
|
(
|
|
const ToPatch& targetPatch,
|
|
const Field
|
|
<
|
|
typename Foam::PrimitivePatch<FaceList, PointField>::point_type
|
|
>& projectionDirection,
|
|
const intersection::algorithm alg,
|
|
const intersection::direction dir
|
|
) const
|
|
{
|
|
// The current patch is slave, i.e. it is being projected onto the target
|
|
|
|
if (projectionDirection.size() != nPoints())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Projection direction field does not correspond to "
|
|
<< "patch points." << endl
|
|
<< "Size: " << projectionDirection.size()
|
|
<< " Number of points: " << nPoints()
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
const labelList& slavePointOrder = localPointOrder();
|
|
|
|
const labelList& slaveMeshPoints = meshPoints();
|
|
|
|
// Result
|
|
List<objectHit> result(nPoints());
|
|
|
|
const labelListList& masterFaceFaces = targetPatch.faceFaces();
|
|
|
|
const ToPatch& masterFaces = targetPatch;
|
|
|
|
const Field<point_type>& masterPoints = targetPatch.points();
|
|
|
|
// Estimate face centre of target side
|
|
Field<point_type> masterFaceCentres(targetPatch.size());
|
|
|
|
forAll(masterFaceCentres, facei)
|
|
{
|
|
masterFaceCentres[facei] =
|
|
average(masterFaces[facei].points(masterPoints));
|
|
}
|
|
|
|
// Algorithm:
|
|
// Loop through all points of the slave side. For every point find the
|
|
// radius for the current contact face. If the contact point falls inside
|
|
// the face and the radius is smaller than for all neighbouring faces,
|
|
// the contact is found. If not, visit the neighbour closest to the
|
|
// calculated contact point. If a single master face is visited more than
|
|
// twice, initiate n-squared search.
|
|
|
|
label curFace = 0;
|
|
label nNSquaredSearches = 0;
|
|
|
|
forAll(slavePointOrder, pointi)
|
|
{
|
|
// Pick up slave point and direction
|
|
const label curLocalPointLabel = slavePointOrder[pointi];
|
|
|
|
const point_type& curPoint =
|
|
points_[slaveMeshPoints[curLocalPointLabel]];
|
|
|
|
const point_type& curProjectionDir =
|
|
projectionDirection[curLocalPointLabel];
|
|
|
|
bool closer;
|
|
|
|
boolList visitedTargetFace(targetPatch.size(), false);
|
|
bool doNSquaredSearch = false;
|
|
|
|
bool foundEligible = false;
|
|
|
|
scalar sqrDistance = GREAT;
|
|
|
|
// Force the full search for the first point to ensure good
|
|
// starting face
|
|
if (pointi == 0)
|
|
{
|
|
doNSquaredSearch = true;
|
|
}
|
|
else
|
|
{
|
|
do
|
|
{
|
|
closer = false;
|
|
doNSquaredSearch = false;
|
|
|
|
// Calculate intersection with curFace
|
|
PointHit<point_type> curHit =
|
|
masterFaces[curFace].ray
|
|
(
|
|
curPoint,
|
|
curProjectionDir,
|
|
masterPoints,
|
|
alg,
|
|
dir
|
|
);
|
|
|
|
visitedTargetFace[curFace] = true;
|
|
|
|
if (curHit.hit())
|
|
{
|
|
result[curLocalPointLabel] = objectHit(true, curFace);
|
|
|
|
break;
|
|
}
|
|
else
|
|
{
|
|
// If a new miss is eligible, it is closer than
|
|
// any previous eligible miss (due to surface walk)
|
|
|
|
// Only grab the miss if it is eligible
|
|
if (curHit.eligibleMiss())
|
|
{
|
|
foundEligible = true;
|
|
result[curLocalPointLabel] = objectHit(false, curFace);
|
|
}
|
|
|
|
// Find the next likely face for intersection
|
|
|
|
// Calculate the miss point on the plane of the
|
|
// face. This is cooked (illogical!) for fastest
|
|
// surface walk.
|
|
//
|
|
point_type missPlanePoint =
|
|
curPoint + curProjectionDir*curHit.distance();
|
|
|
|
const labelList& masterNbrs = masterFaceFaces[curFace];
|
|
|
|
sqrDistance =
|
|
magSqr(missPlanePoint - masterFaceCentres[curFace]);
|
|
|
|
forAll(masterNbrs, nbrI)
|
|
{
|
|
if
|
|
(
|
|
magSqr
|
|
(
|
|
missPlanePoint
|
|
- masterFaceCentres[masterNbrs[nbrI]]
|
|
)
|
|
<= sqrDistance
|
|
)
|
|
{
|
|
closer = true;
|
|
curFace = masterNbrs[nbrI];
|
|
}
|
|
}
|
|
|
|
if (visitedTargetFace[curFace])
|
|
{
|
|
// This face has already been visited.
|
|
// Execute n-squared search
|
|
doNSquaredSearch = true;
|
|
break;
|
|
}
|
|
}
|
|
|
|
DebugInfo << '.';
|
|
} while (closer);
|
|
}
|
|
|
|
if
|
|
(
|
|
doNSquaredSearch || !foundEligible
|
|
)
|
|
{
|
|
nNSquaredSearches++;
|
|
|
|
DebugInfo << "p " << curLocalPointLabel << ": ";
|
|
|
|
result[curLocalPointLabel] = objectHit(false, -1);
|
|
scalar minDistance = GREAT;
|
|
|
|
forAll(masterFaces, facei)
|
|
{
|
|
PointHit<point_type> curHit =
|
|
masterFaces[facei].ray
|
|
(
|
|
curPoint,
|
|
curProjectionDir,
|
|
masterPoints,
|
|
alg,
|
|
dir
|
|
);
|
|
|
|
if (curHit.hit())
|
|
{
|
|
result[curLocalPointLabel] = objectHit(true, facei);
|
|
curFace = facei;
|
|
|
|
break;
|
|
}
|
|
else if (curHit.eligibleMiss())
|
|
{
|
|
// Calculate min distance
|
|
scalar missDist =
|
|
Foam::mag(curHit.missPoint() - curPoint);
|
|
|
|
if (missDist < minDistance)
|
|
{
|
|
minDistance = missDist;
|
|
|
|
result[curLocalPointLabel] = objectHit(false, facei);
|
|
curFace = facei;
|
|
}
|
|
}
|
|
}
|
|
|
|
DebugInfo << result[curLocalPointLabel] << nl;
|
|
}
|
|
else
|
|
{
|
|
DebugInfo << 'x';
|
|
}
|
|
}
|
|
|
|
DebugInfo
|
|
<< nl << "Executed " << nNSquaredSearches
|
|
<< " n-squared searches out of total of "
|
|
<< nPoints() << endl;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
template<class FaceList, class PointField>
|
|
template<class ToPatch>
|
|
Foam::List<Foam::objectHit>
|
|
Foam::PrimitivePatch<FaceList, PointField>::projectFaceCentres
|
|
(
|
|
const ToPatch& targetPatch,
|
|
const Field
|
|
<
|
|
typename Foam::PrimitivePatch<FaceList, PointField>::point_type
|
|
>& projectionDirection,
|
|
const intersection::algorithm alg,
|
|
const intersection::direction dir
|
|
) const
|
|
{
|
|
// The current patch is slave, i.e. it is being projected onto the target
|
|
|
|
if (projectionDirection.size() != this->size())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Projection direction field does not correspond to patch faces."
|
|
<< endl << "Size: " << projectionDirection.size()
|
|
<< " Number of points: " << this->size()
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
labelList slaveFaceOrder = bandCompression(faceFaces());
|
|
|
|
// calculate master face centres
|
|
Field<point_type> masterFaceCentres(targetPatch.size());
|
|
|
|
const labelListList& masterFaceFaces = targetPatch.faceFaces();
|
|
|
|
const ToPatch& masterFaces = targetPatch;
|
|
|
|
const typename ToPatch::PointFieldType& masterPoints = targetPatch.points();
|
|
|
|
forAll(masterFaceCentres, facei)
|
|
{
|
|
masterFaceCentres[facei] =
|
|
masterFaces[facei].centre(masterPoints);
|
|
}
|
|
|
|
// Result
|
|
List<objectHit> result(this->size());
|
|
|
|
const PrimitivePatch<FaceList, PointField>& slaveFaces = *this;
|
|
|
|
const PointField& slaveGlobalPoints = points();
|
|
|
|
// Algorithm:
|
|
// Loop through all points of the slave side. For every point find the
|
|
// radius for the current contact face. If the contact point falls inside
|
|
// the face and the radius is smaller than for all neighbouring faces,
|
|
// the contact is found. If not, visit the neighbour closest to the
|
|
// calculated contact point. If a single master face is visited more than
|
|
// twice, initiate n-squared search.
|
|
|
|
label curFace = 0;
|
|
label nNSquaredSearches = 0;
|
|
|
|
forAll(slaveFaceOrder, facei)
|
|
{
|
|
// pick up slave point and direction
|
|
const label curLocalFaceLabel = slaveFaceOrder[facei];
|
|
|
|
const point& curFaceCentre =
|
|
slaveFaces[curLocalFaceLabel].centre(slaveGlobalPoints);
|
|
|
|
const vector& curProjectionDir =
|
|
projectionDirection[curLocalFaceLabel];
|
|
|
|
bool closer;
|
|
|
|
boolList visitedTargetFace(targetPatch.size(), false);
|
|
bool doNSquaredSearch = false;
|
|
|
|
bool foundEligible = false;
|
|
|
|
scalar sqrDistance = GREAT;
|
|
|
|
// Force the full search for the first point to ensure good
|
|
// starting face
|
|
if (facei == 0)
|
|
{
|
|
doNSquaredSearch = true;
|
|
}
|
|
else
|
|
{
|
|
do
|
|
{
|
|
closer = false;
|
|
doNSquaredSearch = false;
|
|
|
|
// Calculate intersection with curFace
|
|
PointHit<point_type> curHit =
|
|
masterFaces[curFace].ray
|
|
(
|
|
curFaceCentre,
|
|
curProjectionDir,
|
|
masterPoints,
|
|
alg,
|
|
dir
|
|
);
|
|
|
|
visitedTargetFace[curFace] = true;
|
|
|
|
if (curHit.hit())
|
|
{
|
|
result[curLocalFaceLabel] = objectHit(true, curFace);
|
|
|
|
break;
|
|
}
|
|
else
|
|
{
|
|
// If a new miss is eligible, it is closer than
|
|
// any previous eligible miss (due to surface walk)
|
|
|
|
// Only grab the miss if it is eligible
|
|
if (curHit.eligibleMiss())
|
|
{
|
|
foundEligible = true;
|
|
result[curLocalFaceLabel] = objectHit(false, curFace);
|
|
}
|
|
|
|
// Find the next likely face for intersection
|
|
|
|
// Calculate the miss point. This is
|
|
// cooked (illogical!) for fastest surface walk.
|
|
//
|
|
point_type missPlanePoint =
|
|
curFaceCentre + curProjectionDir*curHit.distance();
|
|
|
|
sqrDistance =
|
|
magSqr(missPlanePoint - masterFaceCentres[curFace]);
|
|
|
|
const labelList& masterNbrs = masterFaceFaces[curFace];
|
|
|
|
forAll(masterNbrs, nbrI)
|
|
{
|
|
if
|
|
(
|
|
magSqr
|
|
(
|
|
missPlanePoint
|
|
- masterFaceCentres[masterNbrs[nbrI]]
|
|
)
|
|
<= sqrDistance
|
|
)
|
|
{
|
|
closer = true;
|
|
curFace = masterNbrs[nbrI];
|
|
}
|
|
}
|
|
|
|
if (visitedTargetFace[curFace])
|
|
{
|
|
// This face has already been visited.
|
|
// Execute n-squared search
|
|
doNSquaredSearch = true;
|
|
break;
|
|
}
|
|
}
|
|
|
|
DebugInfo << '.';
|
|
} while (closer);
|
|
}
|
|
|
|
if (doNSquaredSearch || !foundEligible)
|
|
{
|
|
nNSquaredSearches++;
|
|
|
|
DebugInfo << "p " << curLocalFaceLabel << ": ";
|
|
|
|
result[curLocalFaceLabel] = objectHit(false, -1);
|
|
scalar minDistance = GREAT;
|
|
|
|
forAll(masterFaces, facei)
|
|
{
|
|
PointHit<point_type> curHit =
|
|
masterFaces[facei].ray
|
|
(
|
|
curFaceCentre,
|
|
curProjectionDir,
|
|
masterPoints,
|
|
alg,
|
|
dir
|
|
);
|
|
|
|
if (curHit.hit())
|
|
{
|
|
result[curLocalFaceLabel] = objectHit(true, facei);
|
|
curFace = facei;
|
|
|
|
break;
|
|
}
|
|
else if (curHit.eligibleMiss())
|
|
{
|
|
// Calculate min distance
|
|
scalar missDist =
|
|
Foam::mag(curHit.missPoint() - curFaceCentre);
|
|
|
|
if (missDist < minDistance)
|
|
{
|
|
minDistance = missDist;
|
|
|
|
result[curLocalFaceLabel] = objectHit(false, facei);
|
|
curFace = facei;
|
|
}
|
|
}
|
|
}
|
|
|
|
DebugInfo << result[curLocalFaceLabel] << nl;
|
|
}
|
|
else
|
|
{
|
|
DebugInfo << 'x';
|
|
}
|
|
}
|
|
|
|
DebugInfo
|
|
<< nl
|
|
<< "Executed " << nNSquaredSearches
|
|
<< " n-squared searches out of total of "
|
|
<< this->size() << endl;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|