mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
ENH: improve handling of area calculations in faMesh (#2233)
- previous handling did not correctly account for off-processor
connections (SEGFAULT).
- revised the point/area normal calculations to use the dual of the
faces. This simplifies the logic and reduces the depth of loops.
Point corrections from the neighbouring patches is handled
centrally by the new halo face information, which properly handles
on-processor or off-processor faces irrespective of any explicit
processor patches.
STYLE: local function and add comments for finiteArea area weighting
- following the original Tukovic code, the area-weighted normal
for a vector pair (defining a triangle) is weighted by
(1) area : larger weight for larger areas
(2) sin : lower weight for narrow angles (eg, shards)
(3) invDist squared : lower weights for distant points
Refactored to eliminate intermediate operations, some additional
divide-by-zero protection - similar to vector normalise()
This commit is contained in:
Mark Olesen
committed by
Andrew Heather
Andrew Heather
parent
ea92cb82c4
commit
a95427c28a
@ -52,6 +52,8 @@ const Foam::word Foam::faMesh::prefix("finite-area");
|
||||
|
||||
Foam::word Foam::faMesh::meshSubDir = "faMesh";
|
||||
|
||||
int Foam::faMesh::origPointAreaMethod_ = 0; // Tuning
|
||||
|
||||
const int Foam::faMesh::quadricsFit_ = 0; // Tuning
|
||||
|
||||
|
||||
@ -227,7 +229,7 @@ void Foam::faMesh::clearGeomNotAreas() const
|
||||
deleteDemandDrivenData(edgeCentresPtr_);
|
||||
deleteDemandDrivenData(faceAreaNormalsPtr_);
|
||||
deleteDemandDrivenData(edgeAreaNormalsPtr_);
|
||||
deleteDemandDrivenData(pointAreaNormalsPtr_);
|
||||
pointAreaNormalsPtr_.reset(nullptr);
|
||||
deleteDemandDrivenData(faceCurvaturesPtr_);
|
||||
deleteDemandDrivenData(edgeTransformTensorsPtr_);
|
||||
}
|
||||
@ -703,11 +705,20 @@ const Foam::vectorField& Foam::faMesh::pointAreaNormals() const
|
||||
{
|
||||
if (!pointAreaNormalsPtr_)
|
||||
{
|
||||
calcPointAreaNormals();
|
||||
pointAreaNormalsPtr_.reset(new vectorField(nPoints()));
|
||||
|
||||
if (origPointAreaMethod_)
|
||||
{
|
||||
calcPointAreaNormals_orig(*pointAreaNormalsPtr_);
|
||||
}
|
||||
else
|
||||
{
|
||||
calcPointAreaNormals(*pointAreaNormalsPtr_);
|
||||
}
|
||||
|
||||
if (quadricsFit_ > 0)
|
||||
{
|
||||
calcPointAreaNormalsByQuadricsFit();
|
||||
calcPointAreaNormalsByQuadricsFit(*pointAreaNormalsPtr_);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -279,8 +279,8 @@ class faMesh
|
||||
//- Edge area normals
|
||||
mutable edgeVectorField* edgeAreaNormalsPtr_;
|
||||
|
||||
//- Edge area normals
|
||||
mutable vectorField* pointAreaNormalsPtr_;
|
||||
//- Point area normals
|
||||
mutable std::unique_ptr<vectorField> pointAreaNormalsPtr_;
|
||||
|
||||
//- Face curvatures
|
||||
mutable areaScalarField* faceCurvaturesPtr_;
|
||||
@ -309,6 +309,9 @@ class faMesh
|
||||
|
||||
// Static Private Data
|
||||
|
||||
//- Use the original method for point normals
|
||||
static int origPointAreaMethod_;
|
||||
|
||||
//- Use quadrics fit
|
||||
static const int quadricsFit_;
|
||||
|
||||
@ -372,10 +375,13 @@ class faMesh
|
||||
void calcEdgeAreaNormals() const;
|
||||
|
||||
//- Calculate point area normals
|
||||
void calcPointAreaNormals() const;
|
||||
void calcPointAreaNormals_orig(vectorField& result) const;
|
||||
|
||||
//- Calculate point area normals
|
||||
void calcPointAreaNormals(vectorField& result) const;
|
||||
|
||||
//- Calculate point area normals by quadrics fit
|
||||
void calcPointAreaNormalsByQuadricsFit() const;
|
||||
void calcPointAreaNormalsByQuadricsFit(vectorField& result) const;
|
||||
|
||||
//- Calculate face curvatures
|
||||
void calcFaceCurvatures() const;
|
||||
@ -448,7 +454,6 @@ class faMesh
|
||||
}
|
||||
|
||||
|
||||
|
||||
public:
|
||||
|
||||
// Public Typedefs
|
||||
|
||||
@ -40,9 +40,73 @@ License
|
||||
#include "processorFaPatchFields.H"
|
||||
#include "emptyFaPatchFields.H"
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
// Define an area-weighted normal for three points (defining a triangle)
|
||||
// (p0, p1, p2) are the base, first, second points respectively
|
||||
//
|
||||
// From the original Tukovic code:
|
||||
//
|
||||
// vector n = (d1^d2)/mag(d1^d2);
|
||||
// scalar sinAlpha = mag(d1^d2)/(mag(d1)*mag(d2));
|
||||
// scalar w = sinAlpha/(mag(d1)*mag(d2));
|
||||
//
|
||||
// ie, normal weighted by area, sine angle and inverse distance squared.
|
||||
// - area : larger weight for larger areas
|
||||
// - sin : lower weight for narrow angles (eg, shards)
|
||||
// - inv distance squared : lower weights for distant points
|
||||
//
|
||||
// The above refactored, with 0.5 for area:
|
||||
//
|
||||
// (d1 ^ d2) / (2 * magSqr(d1) * magSqr(d2))
|
||||
|
||||
static inline vector areaInvDistSqrWeightedNormal
|
||||
(
|
||||
const vector& a,
|
||||
const vector& b
|
||||
)
|
||||
{
|
||||
const scalar s(2*magSqr(a)*magSqr(b));
|
||||
|
||||
return s < VSMALL ? Zero : (a ^ b) / s;
|
||||
}
|
||||
|
||||
|
||||
// The area normal for the face dual (around base-point)
|
||||
// described by the right/left edge points and the centre point
|
||||
//
|
||||
// The adjustment for 1/2 edge point (the dual point) is done internally
|
||||
static inline vector areaInvDistSqrWeightedNormalDualEdge
|
||||
(
|
||||
const point& basePoint,
|
||||
const point& rightPoint,
|
||||
const point& leftPoint,
|
||||
const point& centrePoint
|
||||
)
|
||||
{
|
||||
const vector mid(centrePoint - basePoint);
|
||||
|
||||
return
|
||||
(
|
||||
areaInvDistSqrWeightedNormal
|
||||
(
|
||||
0.5*(rightPoint - basePoint), // vector to right 1/2 edge
|
||||
mid
|
||||
)
|
||||
+ areaInvDistSqrWeightedNormal
|
||||
(
|
||||
mid,
|
||||
0.5*(leftPoint - basePoint) // vector to left 1/2 edge
|
||||
)
|
||||
);
|
||||
}
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
@ -61,6 +125,7 @@ static Foam::bitSet markupBoundaryPoints(const uindirectPrimitivePatch& p)
|
||||
return markPoints;
|
||||
}
|
||||
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
|
||||
@ -496,108 +561,40 @@ void Foam::faMesh::calcEdgeAreaNormals() const
|
||||
// Point area normals
|
||||
const vectorField& pointNormals = pointAreaNormals();
|
||||
|
||||
|
||||
// // Primitive patch edge normals
|
||||
// const labelListList& patchPointEdges = patch().pointEdges();
|
||||
|
||||
// vectorField patchEdgeNormals(nEdges(), Zero);
|
||||
|
||||
// forAll(pointNormals, pointI)
|
||||
// {
|
||||
// const labelList& curPointEdges = patchPointEdges[pointI];
|
||||
|
||||
// forAll(curPointEdges, edgeI)
|
||||
// {
|
||||
// label curEdge = curPointEdges[edgeI];
|
||||
|
||||
// patchEdgeNormals[curEdge] += 0.5*pointNormals[pointI];
|
||||
// }
|
||||
// }
|
||||
|
||||
// patchEdgeNormals /= mag(patchEdgeNormals);
|
||||
|
||||
|
||||
// // Edge area normals
|
||||
// label nIntEdges = patch().nInternalEdges();
|
||||
|
||||
// for (label edgeI = 0; edgeI < nIntEdges; ++edgeI)
|
||||
// {
|
||||
// edgeAreaNormals.ref()[edgeI] =
|
||||
// patchEdgeNormals[edgeI];
|
||||
// }
|
||||
|
||||
// forAll(boundary(), patchI)
|
||||
// {
|
||||
// const labelList& edgeLabels = boundary()[patchI];
|
||||
|
||||
// forAll(edgeAreaNormals.boundaryFieldRef()[patchI], edgeI)
|
||||
// {
|
||||
// edgeAreaNormals.boundaryFieldRef()[patchI][edgeI] =
|
||||
// patchEdgeNormals[edgeLabels[edgeI]];
|
||||
// }
|
||||
// }
|
||||
|
||||
|
||||
forAll(edgeAreaNormals.internalField(), edgeI)
|
||||
forAll(edgeAreaNormals.internalField(), edgei)
|
||||
{
|
||||
const vector e = edges()[edgeI].unitVec(points());
|
||||
const edge& e = edges()[edgei];
|
||||
const vector edgeVec = e.unitVec(points());
|
||||
|
||||
// scalar wStart =
|
||||
// 1.0 - sqr(mag(e^pointNormals[edges()[edgeI].end()]));
|
||||
vector& n = edgeAreaNormals.ref()[edgei];
|
||||
|
||||
// scalar wEnd =
|
||||
// 1.0 - sqr(mag(e^pointNormals[edges()[edgeI].start()]));
|
||||
n = (pointNormals[e.first()] + pointNormals[e.second()]);
|
||||
|
||||
// wStart = 1.0;
|
||||
// wEnd = 1.0;
|
||||
n -= edgeVec*(edgeVec & n);
|
||||
|
||||
// edgeAreaNormals.ref()[edgeI] =
|
||||
// wStart*pointNormals[edges()[edgeI].start()]
|
||||
// + wEnd*pointNormals[edges()[edgeI].end()];
|
||||
|
||||
// vector eC =
|
||||
// 0.5
|
||||
// *(
|
||||
// points()[edges()[edgeI].start()]
|
||||
// + points()[edges()[edgeI].end()]
|
||||
// );
|
||||
|
||||
// vector eCp = 0.5*
|
||||
// (
|
||||
// points()[edges()[edgeI].start()]
|
||||
// + pointNormals[edges()[edgeI].start()]
|
||||
// points()[edges()[edgeI].end()] +
|
||||
// );
|
||||
|
||||
edgeAreaNormals.ref()[edgeI] =
|
||||
pointNormals[edges()[edgeI].start()]
|
||||
+ pointNormals[edges()[edgeI].end()];
|
||||
|
||||
edgeAreaNormals.ref()[edgeI] -=
|
||||
e*(e&edgeAreaNormals.internalField()[edgeI]);
|
||||
n.normalise();
|
||||
}
|
||||
|
||||
edgeAreaNormals.ref() /= mag(edgeAreaNormals.internalField());
|
||||
|
||||
forAll(boundary(), patchI)
|
||||
forAll(boundary(), patchi)
|
||||
{
|
||||
const edgeList::subList patchEdges =
|
||||
boundary()[patchI].patchSlice(edges());
|
||||
boundary()[patchi].patchSlice(edges());
|
||||
|
||||
forAll(patchEdges, edgeI)
|
||||
vectorField& edgeNorms = edgeAreaNormals.boundaryFieldRef()[patchi];
|
||||
|
||||
forAll(patchEdges, edgei)
|
||||
{
|
||||
edgeAreaNormals.boundaryFieldRef()[patchI][edgeI] =
|
||||
pointNormals[patchEdges[edgeI].start()]
|
||||
+ pointNormals[patchEdges[edgeI].end()];
|
||||
const edge& e = patchEdges[edgei];
|
||||
const vector edgeVec = e.unitVec(points());
|
||||
|
||||
const vector e = patchEdges[edgeI].unitVec(points());
|
||||
vector& n = edgeNorms[edgei];
|
||||
|
||||
edgeAreaNormals.boundaryFieldRef()[patchI][edgeI] -=
|
||||
e*(e&edgeAreaNormals.boundaryField()[patchI][edgeI]);
|
||||
n = (pointNormals[e.first()] + pointNormals[e.second()]);
|
||||
|
||||
n -= edgeVec*(edgeVec & n);
|
||||
|
||||
n.normalise();
|
||||
}
|
||||
|
||||
edgeAreaNormals.boundaryFieldRef()[patchI] /=
|
||||
mag(edgeAreaNormals.boundaryField()[patchI]);
|
||||
}
|
||||
}
|
||||
|
||||
@ -880,19 +877,53 @@ Foam::labelList Foam::faMesh::boundaryPoints() const
|
||||
}
|
||||
|
||||
|
||||
void Foam::faMesh::calcPointAreaNormals() const
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
// Point normal calculations
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
// Original method (general)
|
||||
// -------------------------
|
||||
// - For each point, obtain the list of connected patch faces
|
||||
// (from point-to-face addressing).
|
||||
//
|
||||
// - Create a primitive patch for those faces and use that to obtain the
|
||||
// outer edge loop(s). This is effectively an agglomeration of the patch
|
||||
// faces connected to a point.
|
||||
//
|
||||
// - Perform a pair-wise walk of the edge loop to obtain triangles from
|
||||
// the originating point outwards (fan-like triangulation).
|
||||
// Calculate an area-weighted value for each triangle.
|
||||
//
|
||||
// NOTE: not sure why internal and boundary point agglomeration was
|
||||
// handled separately.
|
||||
//
|
||||
// Problems:
|
||||
// - possibly susceptible to edge-loop errors (issue #2233) that cause
|
||||
// the agglomeration logic to include the current point twice?
|
||||
// - use of outer edge loop makes it more sensitive to face warpage.
|
||||
// - relatively expensive with point-face connectivity,
|
||||
// creation/destruction of a primitive-patch around each point.
|
||||
//
|
||||
// Original method (boundary correction)
|
||||
// -------------------------------------
|
||||
//
|
||||
// - correct wedge directly, use processor patch information to exchange
|
||||
// the current summed values
|
||||
//
|
||||
// - explicit correction of other boundaries.
|
||||
// Polls the patch for the ngbPolyPatchPointNormals(), which internally
|
||||
// calls ngbPolyPatchFaces and can return -1 for unmatched edges.
|
||||
// This occurs when the outside perimeter of the faPatch aligns with
|
||||
// a polyMesh processor. The neighbour face is off-processor and cannot
|
||||
// be found. Accessing the mesh face at -1 == SEGFAULT.
|
||||
|
||||
void Foam::faMesh::calcPointAreaNormals_orig(vectorField& result) const
|
||||
{
|
||||
if (pointAreaNormalsPtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "pointAreaNormalsPtr_ already allocated"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
DebugInFunction
|
||||
<< "Calculating pointAreaNormals : original method" << endl;
|
||||
|
||||
|
||||
pointAreaNormalsPtr_ = new vectorField(nPoints(), Zero);
|
||||
|
||||
vectorField& result = *pointAreaNormalsPtr_;
|
||||
result.resize(nPoints());
|
||||
result = Zero;
|
||||
|
||||
labelList intPoints(internalPoints());
|
||||
labelList bndPoints(boundaryPoints());
|
||||
@ -1004,14 +1035,14 @@ void Foam::faMesh::calcPointAreaNormals() const
|
||||
|
||||
const labelList& patchPoints = wedgePatch.pointLabels();
|
||||
|
||||
vector N =
|
||||
const vector N
|
||||
(
|
||||
transform
|
||||
(
|
||||
wedgePatch.edgeT(),
|
||||
wedgePatch.centreNormal()
|
||||
);
|
||||
|
||||
N /= mag(N);
|
||||
).normalise()
|
||||
);
|
||||
|
||||
for (const label pti : patchPoints)
|
||||
{
|
||||
@ -1157,14 +1188,302 @@ void Foam::faMesh::calcPointAreaNormals() const
|
||||
}
|
||||
}
|
||||
|
||||
result /= mag(result);
|
||||
for (vector& n : result)
|
||||
{
|
||||
n.normalise();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::faMesh::calcPointAreaNormalsByQuadricsFit() const
|
||||
{
|
||||
vectorField& result = *pointAreaNormalsPtr_;
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
// Point normal calculations
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
// Revised method (general)
|
||||
// ------------------------
|
||||
//
|
||||
// - For each patch face obtain a centre point (mathematical avg)
|
||||
// and use that to define the face dual as a pair of triangles:
|
||||
// - tri1: base-point / mid-side of right edge / face centre
|
||||
// - tri2: base-point / face centre / mid-side of left edge
|
||||
//
|
||||
// - Walk all face points, inserting directly into the corresponding
|
||||
// locations. No distinction between internal or boundary points (yet).
|
||||
//
|
||||
// Revised method (boundary correction)
|
||||
// ------------------------------------
|
||||
//
|
||||
// - correct wedge directly, use processor patch information to exchange
|
||||
// the current summed values. [No change from original].
|
||||
//
|
||||
// - explicit correction of other boundaries.
|
||||
// Use the new boundary halo information for the face normals.
|
||||
// Calculate the equivalent face-point normals locally and apply
|
||||
// correction as before.
|
||||
|
||||
void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
|
||||
{
|
||||
DebugInFunction
|
||||
<< "Calculating pointAreaNormals : face dual method" << endl;
|
||||
|
||||
result.resize(nPoints());
|
||||
result = Zero;
|
||||
|
||||
const pointField& points = patch().localPoints();
|
||||
const faceList& faces = patch().localFaces();
|
||||
|
||||
|
||||
// Loop over all faces
|
||||
|
||||
for (const face& f : faces)
|
||||
{
|
||||
const label nVerts(f.size());
|
||||
|
||||
point centrePoint(Zero);
|
||||
for (label i = 0; i < nVerts; ++i)
|
||||
{
|
||||
centrePoint += points[f[i]];
|
||||
}
|
||||
centrePoint /= nVerts;
|
||||
|
||||
for (label i = 0; i < nVerts; ++i)
|
||||
{
|
||||
const label pt0 = f.thisLabel(i); // base
|
||||
|
||||
result[pt0] +=
|
||||
areaInvDistSqrWeightedNormalDualEdge
|
||||
(
|
||||
points[pt0], // base
|
||||
points[f.nextLabel(i)], // right
|
||||
points[f.prevLabel(i)], // left
|
||||
centrePoint
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Handle the boundary edges
|
||||
|
||||
bitSet nbrBoundaryAdjust(boundary().size(), true);
|
||||
|
||||
forAll(boundary(), patchi)
|
||||
{
|
||||
const faPatch& fap = boundary()[patchi];
|
||||
|
||||
if (isA<wedgeFaPatch>(fap))
|
||||
{
|
||||
// Correct wedge points
|
||||
|
||||
const auto& wedgePatch = refCast<const wedgeFaPatch>(fap);
|
||||
|
||||
const labelList& patchPoints = wedgePatch.pointLabels();
|
||||
|
||||
const vector N
|
||||
(
|
||||
transform
|
||||
(
|
||||
wedgePatch.edgeT(),
|
||||
wedgePatch.centreNormal()
|
||||
).normalise()
|
||||
);
|
||||
|
||||
for (const label pti : patchPoints)
|
||||
{
|
||||
result[pti] -= N*(N&result[pti]);
|
||||
}
|
||||
|
||||
// Axis point correction
|
||||
if (wedgePatch.axisPoint() > -1)
|
||||
{
|
||||
result[wedgePatch.axisPoint()] =
|
||||
wedgePatch.axis()
|
||||
*(
|
||||
wedgePatch.axis()
|
||||
&result[wedgePatch.axisPoint()]
|
||||
);
|
||||
}
|
||||
|
||||
// Handled
|
||||
nbrBoundaryAdjust.unset(patchi);
|
||||
}
|
||||
else if (Pstream::parRun() && isA<processorFaPatch>(fap))
|
||||
{
|
||||
// Correct processor patch points
|
||||
|
||||
const auto& procPatch = refCast<const processorFaPatch>(fap);
|
||||
|
||||
const labelList& patchPoints = procPatch.pointLabels();
|
||||
const labelList& nbrPatchPoints = procPatch.neighbPoints();
|
||||
|
||||
const labelList& nonGlobalPatchPoints =
|
||||
procPatch.nonGlobalPatchPoints();
|
||||
|
||||
// Send my values
|
||||
|
||||
vectorField patchPointNormals
|
||||
(
|
||||
UIndirectList<vector>(result, patchPoints)
|
||||
);
|
||||
|
||||
{
|
||||
OPstream::write
|
||||
(
|
||||
Pstream::commsTypes::blocking,
|
||||
procPatch.neighbProcNo(),
|
||||
reinterpret_cast<const char*>(patchPointNormals.cdata()),
|
||||
patchPointNormals.size_bytes()
|
||||
);
|
||||
}
|
||||
|
||||
// Receive neighbour values
|
||||
patchPointNormals.resize(nbrPatchPoints.size());
|
||||
|
||||
{
|
||||
IPstream::read
|
||||
(
|
||||
Pstream::commsTypes::blocking,
|
||||
procPatch.neighbProcNo(),
|
||||
reinterpret_cast<char*>(patchPointNormals.data()),
|
||||
patchPointNormals.size_bytes()
|
||||
);
|
||||
}
|
||||
|
||||
for (const label pti : nonGlobalPatchPoints)
|
||||
{
|
||||
result[patchPoints[pti]] +=
|
||||
patchPointNormals[nbrPatchPoints[pti]];
|
||||
}
|
||||
|
||||
// Handled
|
||||
nbrBoundaryAdjust.unset(patchi);
|
||||
}
|
||||
else if (fap.coupled())
|
||||
{
|
||||
// Coupled - no further action for neighbour side
|
||||
nbrBoundaryAdjust.unset(patchi);
|
||||
}
|
||||
// TBD:
|
||||
/// else if (isA<emptyFaPatch>(fap))
|
||||
/// {
|
||||
/// // Ignore this boundary
|
||||
/// nbrBoundaryAdjust.unset(patchi);
|
||||
/// }
|
||||
else if (!correctPatchPointNormals(patchi))
|
||||
{
|
||||
// No corrections
|
||||
nbrBoundaryAdjust.unset(patchi);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Correct global points
|
||||
if (globalData().nGlobalPoints())
|
||||
{
|
||||
const labelList& spLabels = globalData().sharedPointLabels();
|
||||
const labelList& addr = globalData().sharedPointAddr();
|
||||
|
||||
vectorField spNormals
|
||||
(
|
||||
UIndirectList<vector>(result, spLabels)
|
||||
);
|
||||
|
||||
vectorField gpNormals
|
||||
(
|
||||
globalData().nGlobalPoints(),
|
||||
Zero
|
||||
);
|
||||
|
||||
forAll(addr, i)
|
||||
{
|
||||
gpNormals[addr[i]] += spNormals[i];
|
||||
}
|
||||
|
||||
combineReduce(gpNormals, plusEqOp<vectorField>());
|
||||
|
||||
// Extract local data
|
||||
forAll(addr, i)
|
||||
{
|
||||
spNormals[i] = gpNormals[addr[i]];
|
||||
}
|
||||
|
||||
forAll(spNormals, pointI)
|
||||
{
|
||||
result[spLabels[pointI]] = spNormals[pointI];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (returnReduce(nbrBoundaryAdjust.any(), orOp<bool>()))
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
PoutInFunction
|
||||
<< "Apply " << nbrBoundaryAdjust.count()
|
||||
<< " boundary neighbour corrections" << nl;
|
||||
}
|
||||
|
||||
// Apply boundary points correction
|
||||
|
||||
// Collect face normals as point normals
|
||||
|
||||
const auto& haloNormals = this->haloFaceNormals();
|
||||
|
||||
Map<vector> fpNormals(4*nBoundaryEdges());
|
||||
|
||||
for (const label patchi : nbrBoundaryAdjust)
|
||||
{
|
||||
const faPatch& fap = boundary()[patchi];
|
||||
const labelList& edgeLabels = fap.edgeLabels();
|
||||
|
||||
if (fap.ngbPolyPatchIndex() < 0)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Neighbour polyPatch index is not defined "
|
||||
<< "for faPatch " << fap.name()
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
for (const label edgei : edgeLabels)
|
||||
{
|
||||
const edge& e = patch().edges()[edgei];
|
||||
|
||||
// Halo face unitNormal at boundary edge (starts as 0)
|
||||
const vector& fnorm = haloNormals[edgei - nInternalEdges_];
|
||||
|
||||
fpNormals(e.first()) += fnorm;
|
||||
fpNormals(e.second()) += fnorm;
|
||||
}
|
||||
}
|
||||
|
||||
// Apply the correction
|
||||
|
||||
// Note from Zeljko Tukovic:
|
||||
//
|
||||
// This posibility is used for free-surface tracking
|
||||
// calculations to enforce 90 deg contact angle between
|
||||
// finite-area mesh and neighbouring polyPatch. It is very
|
||||
// important for curvature calculation to have correct normal
|
||||
// at contact line points.
|
||||
|
||||
forAllConstIters(fpNormals, iter)
|
||||
{
|
||||
const label pointi = iter.key();
|
||||
vector fpnorm = iter.val();
|
||||
|
||||
fpnorm.normalise();
|
||||
result[pointi] -= fpnorm*(fpnorm & result[pointi]);
|
||||
}
|
||||
}
|
||||
|
||||
for (vector& n : result)
|
||||
{
|
||||
n.normalise();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
|
||||
{
|
||||
const labelList intPoints(internalPoints());
|
||||
const labelList bndPoints(boundaryPoints());
|
||||
|
||||
@ -1724,7 +2043,10 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit() const
|
||||
}
|
||||
}
|
||||
|
||||
result /= mag(result);
|
||||
for (vector& n : result)
|
||||
{
|
||||
n.normalise();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -32,6 +32,7 @@ License
|
||||
#include "faMesh.H"
|
||||
#include "areaFields.H"
|
||||
#include "edgeFields.H"
|
||||
#include "edgeHashes.H"
|
||||
#include "polyMesh.H"
|
||||
#include "demandDrivenData.H"
|
||||
|
||||
@ -213,7 +214,7 @@ void Foam::faPatch::calcPointLabels() const
|
||||
{
|
||||
SLList<label> labels;
|
||||
|
||||
UList<edge> edges = patchSlice(boundaryMesh().mesh().edges());
|
||||
const edgeList::subList edges = patchSlice(boundaryMesh().mesh().edges());
|
||||
|
||||
forAll(edges, edgeI)
|
||||
{
|
||||
@ -302,58 +303,6 @@ const Foam::labelListList& Foam::faPatch::pointEdges() const
|
||||
}
|
||||
|
||||
|
||||
Foam::labelList Foam::faPatch::ngbPolyPatchFaces() const
|
||||
{
|
||||
if (nbrPolyPatchId_ < 0)
|
||||
{
|
||||
return labelList();
|
||||
}
|
||||
|
||||
labelList ngbFaces(faPatch::size());
|
||||
|
||||
const faMesh& aMesh = boundaryMesh().mesh();
|
||||
const polyMesh& pMesh = aMesh.mesh();
|
||||
const auto& patch = aMesh.patch();
|
||||
|
||||
const labelListList& edgeFaces = pMesh.edgeFaces();
|
||||
|
||||
const labelList meshEdges
|
||||
(
|
||||
patch.meshEdges(pMesh.edges(), pMesh.pointEdges())
|
||||
);
|
||||
|
||||
forAll(ngbFaces, edgeI)
|
||||
{
|
||||
ngbFaces[edgeI] = -1;
|
||||
|
||||
label curEdge = (*this)[edgeI];
|
||||
|
||||
label curPMeshEdge = meshEdges[curEdge];
|
||||
|
||||
forAll(edgeFaces[curPMeshEdge], faceI)
|
||||
{
|
||||
label curFace = edgeFaces[curPMeshEdge][faceI];
|
||||
|
||||
label curPatchID = pMesh.boundaryMesh().whichPatch(curFace);
|
||||
|
||||
if (curPatchID == nbrPolyPatchId_)
|
||||
{
|
||||
ngbFaces[edgeI] = curFace;
|
||||
}
|
||||
}
|
||||
|
||||
if (ngbFaces[edgeI] == -1)
|
||||
{
|
||||
WarningInFunction
|
||||
<< "Problem with determination of edge ngb faces!"
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
|
||||
return ngbFaces;
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const
|
||||
{
|
||||
if (nbrPolyPatchId_ < 0)
|
||||
@ -361,24 +310,7 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const
|
||||
return tmp<vectorField>::New();
|
||||
}
|
||||
|
||||
auto tfN = tmp<vectorField>::New();
|
||||
auto& fN = tfN.ref();
|
||||
|
||||
fN.setSize(faPatch::size());
|
||||
|
||||
labelList ngbFaces = ngbPolyPatchFaces();
|
||||
|
||||
const polyMesh& pMesh = boundaryMesh().mesh()();
|
||||
|
||||
const faceList& faces = pMesh.faces();
|
||||
const pointField& points = pMesh.points();
|
||||
|
||||
forAll(fN, faceI)
|
||||
{
|
||||
fN[faceI] = faces[ngbFaces[faceI]].unitNormal(points);
|
||||
}
|
||||
|
||||
return tfN;
|
||||
return boundaryMesh().mesh().haloFaceNormals(this->index());
|
||||
}
|
||||
|
||||
|
||||
@ -389,24 +321,31 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchPointNormals() const
|
||||
return tmp<vectorField>::New();
|
||||
}
|
||||
|
||||
// Unit normals for the neighbour patch faces
|
||||
const vectorField faceNormals
|
||||
(
|
||||
boundaryMesh().mesh().haloFaceNormals(this->index())
|
||||
);
|
||||
|
||||
const labelListList& pntEdges = pointEdges();
|
||||
|
||||
auto tpN = tmp<vectorField>::New(pntEdges.size(), Zero);
|
||||
auto& pN = tpN.ref();
|
||||
auto tpointNorm = tmp<vectorField>::New(pntEdges.size());
|
||||
auto& pointNorm = tpointNorm.ref();
|
||||
|
||||
const vectorField faceNormals(ngbPolyPatchFaceNormals());
|
||||
|
||||
forAll(pN, pointI)
|
||||
forAll(pointNorm, pointi)
|
||||
{
|
||||
forAll(pntEdges[pointI], edgeI)
|
||||
vector& n = pointNorm[pointi];
|
||||
n = Zero;
|
||||
|
||||
for (const label bndEdgei : pntEdges[pointi])
|
||||
{
|
||||
pN[pointI] += faceNormals[pntEdges[pointI][edgeI]];
|
||||
n += faceNormals[bndEdgei];
|
||||
}
|
||||
|
||||
n.normalise();
|
||||
}
|
||||
|
||||
pN /= mag(pN);
|
||||
|
||||
return tpN;
|
||||
return tpointNorm;
|
||||
}
|
||||
|
||||
|
||||
@ -444,11 +383,14 @@ const Foam::scalarField& Foam::faPatch::magEdgeLengths() const
|
||||
|
||||
Foam::tmp<Foam::vectorField> Foam::faPatch::edgeNormals() const
|
||||
{
|
||||
tmp<vectorField> eN(new vectorField(size()));
|
||||
tmp<vectorField> tedgeNorm(edgeLengths());
|
||||
|
||||
eN.ref() = edgeLengths()/magEdgeLengths();
|
||||
for (vector& n : tedgeNorm.ref())
|
||||
{
|
||||
n.normalise();
|
||||
}
|
||||
|
||||
return eN;
|
||||
return tedgeNorm;
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -311,10 +311,8 @@ public:
|
||||
//- Return patch point-edge addressing
|
||||
const labelListList& pointEdges() const;
|
||||
|
||||
//- Return edge neighbour polyPatch faces
|
||||
labelList ngbPolyPatchFaces() const;
|
||||
|
||||
//- Return normals of neighbour polyPatch faces
|
||||
// Same as faMesh::haloFaceNormals()
|
||||
tmp<vectorField> ngbPolyPatchFaceNormals() const;
|
||||
|
||||
//- Return normals of neighbour polyPatch joined points
|
||||
|
||||
Reference in New Issue
Block a user