From 4094e2c58f38c614e3072dad6209e353647cb821 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 12 Jun 2017 14:36:45 +0200
Subject: [PATCH] ENH: reduce memory overhead in mergePoints (issue #495)
- can avoid the intermediate point distance field in exchange for
calculating a point subtraction twice.
---
src/OpenFOAM/meshes/meshTools/mergePoints.C | 61 +++++++++++----------
1 file changed, 32 insertions(+), 29 deletions(-)
diff --git a/src/OpenFOAM/meshes/meshTools/mergePoints.C b/src/OpenFOAM/meshes/meshTools/mergePoints.C
index 2b99235b72b..041a6f81944 100644
--- a/src/OpenFOAM/meshes/meshTools/mergePoints.C
+++ b/src/OpenFOAM/meshes/meshTools/mergePoints.C
@@ -41,22 +41,28 @@ Foam::label Foam::mergePoints
{
typedef typename PointList::value_type point_type;
- // Create a old to new point mapping array
- pointMap.setSize(points.size());
+ const label nPoints = points.size();
+
+ // Create an old to new point mapping array
+ pointMap.setSize(nPoints);
pointMap = -1;
- if (points.empty())
+ if (!nPoints)
{
return 0;
}
- // Explicitly convert to Field to support various list types
- tmp<Field<point_type>> tPoints(new Field<point_type>(points));
-
point_type compareOrigin = origin;
if (origin == point_type::max)
{
- compareOrigin = sum(tPoints())/points.size();
+ // Use average of input points to define a comparison origin.
+ // Same as sum(points)/nPoints, but handles different list types
+ compareOrigin = points[0];
+ for (label pointi=1; pointi < nPoints; ++pointi)
+ {
+ compareOrigin += points[pointi];
+ }
+ compareOrigin /= nPoints;
}
// We're comparing distance squared to origin first.
@@ -68,34 +74,31 @@ Foam::label Foam::mergePoints
// x^2+y^2+z^2 + 2*mergeTol*(x+z+y) + mergeTol^2*...
// so the difference will be 2*mergeTol*(x+y+z)
- const scalar mergeTolSqr = Foam::sqr(scalar(mergeTol));
+ const scalar mergeTolSqr = Foam::sqr(mergeTol);
// Sort points by magSqr
- const Field<point_type> d(tPoints - compareOrigin);
-
- List<scalar> magSqrD(d.size());
- forAll(d, pointi)
+ List<scalar> magSqrDist(nPoints);
+ forAll(points, pointi)
{
- magSqrD[pointi] = magSqr(d[pointi]);
+ magSqrDist[pointi] = magSqr(points[pointi] - compareOrigin);
}
labelList order;
- Foam::sortedOrder(magSqrD, order);
+ Foam::sortedOrder(magSqrDist, order);
- Field<scalar> sortedTol(points.size());
+ Field<scalar> sortedTol(nPoints);
forAll(order, sortI)
{
- const label pointi = order[sortI];
+ const point_type& pt = points[order[sortI]];
- // Convert to scalar precision
- // NOTE: not yet using point_type template parameter
- const point pt
- (
- scalar(d[pointi].x()),
- scalar(d[pointi].y()),
- scalar(d[pointi].z())
- );
- sortedTol[sortI] = 2*mergeTol*(mag(pt.x())+mag(pt.y())+mag(pt.z()));
+ // Use scalar precision
+ sortedTol[sortI] =
+ 2*mergeTol*
+ (
+ mag(scalar(pt.x() - compareOrigin.x())),
+ + mag(scalar(pt.y() - compareOrigin.y())),
+ + mag(scalar(pt.z() - compareOrigin.z()))
+ );
}
label newPointi = 0;
@@ -104,11 +107,11 @@ Foam::label Foam::mergePoints
label pointi = order[0];
pointMap[pointi] = newPointi++;
- for (label sortI = 1; sortI < order.size(); sortI++)
+ for (label sortI = 1; sortI < order.size(); ++sortI)
{
// Get original point index
const label pointi = order[sortI];
- const scalar mag2 = magSqrD[order[sortI]];
+ const scalar mag2 = magSqrDist[order[sortI]];
// Convert to scalar precision
// NOTE: not yet using point_type template parameter
@@ -127,8 +130,8 @@ Foam::label Foam::mergePoints
(
label prevSortI = sortI - 1;
prevSortI >= 0
- && (mag(magSqrD[order[prevSortI]] - mag2) <= sortedTol[sortI]);
- prevSortI--
+ && (mag(magSqrDist[order[prevSortI]] - mag2) <= sortedTol[sortI]);
+ --prevSortI
)
{
const label prevPointi = order[prevSortI];
--
GitLab