mergePoints.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "ListOps.H"
27 #include "point.H"
28 #include "Field.H"
29 
30 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
31 
32 template<class Type>
34 (
35  const UList<Type>& points,
36  const scalar mergeTol,
37  const bool verbose,
38  labelList& pointMap,
39  const Type& origin
40 )
41 {
42  Type compareOrigin = origin;
43 
44  if (origin == Type::max)
45  {
46  if (points.size())
47  {
48  compareOrigin = sum(points)/points.size();
49  }
50  }
51 
52  // Create a old to new point mapping array
53  pointMap.setSize(points.size());
54  pointMap = -1;
55 
56  if (points.empty())
57  {
58  return points.size();
59  }
60 
61  // We're comparing distance squared to origin first.
62  // Say if starting from two close points:
63  // x, y, z
64  // x+mergeTol, y+mergeTol, z+mergeTol
65  // Then the magSqr of both will be
66  // x^2+y^2+z^2
67  // x^2+y^2+z^2 + 2*mergeTol*(x+z+y) + mergeTol^2*...
68  // so the difference will be 2*mergeTol*(x+y+z)
69 
70  const scalar mergeTolSqr = Foam::sqr(scalar(mergeTol));
71 
72  // Sort points by magSqr
73  const Field<Type> d(points - compareOrigin);
74 
75  List<scalar> magSqrD(d.size());
76  forAll(d, pointi)
77  {
78  magSqrD[pointi] = magSqr(d[pointi]);
79  }
80  labelList order;
81  sortedOrder(magSqrD, order);
82 
83 
84  Field<scalar> sortedTol(points.size());
85  forAll(order, sortI)
86  {
87  label pointi = order[sortI];
88 
89  // Convert to scalar precision
90  const point pt
91  (
92  scalar(d[pointi].x()),
93  scalar(d[pointi].y()),
94  scalar(d[pointi].z())
95  );
96  sortedTol[sortI] = 2*mergeTol*(mag(pt.x())+mag(pt.y())+mag(pt.z()));
97  }
98 
99  label newPointi = 0;
100 
101 
102  // Handle 0th point separately (is always unique)
103  label pointi = order[0];
104  pointMap[pointi] = newPointi++;
105 
106 
107  for (label sortI = 1; sortI < order.size(); sortI++)
108  {
109  // Get original point index
110  label pointi = order[sortI];
111  const scalar mag2 = magSqrD[order[sortI]];
112  // Convert to scalar precision
113  const point pt
114  (
115  scalar(points[pointi].x()),
116  scalar(points[pointi].y()),
117  scalar(points[pointi].z())
118  );
119 
120 
121  // Compare to previous points to find equal one.
122  label equalPointi = -1;
123 
124  for
125  (
126  label prevSortI = sortI - 1;
127  prevSortI >= 0
128  && (mag(magSqrD[order[prevSortI]] - mag2) <= sortedTol[sortI]);
129  prevSortI--
130  )
131  {
132  label prevPointi = order[prevSortI];
133  const point prevPt
134  (
135  scalar(points[prevPointi].x()),
136  scalar(points[prevPointi].y()),
137  scalar(points[prevPointi].z())
138  );
139 
140  if (magSqr(pt - prevPt) <= mergeTolSqr)
141  {
142  // Found match.
143  equalPointi = prevPointi;
144 
145  break;
146  }
147  }
148 
149 
150  if (equalPointi != -1)
151  {
152  // Same coordinate as equalPointi. Map to same new point.
153  pointMap[pointi] = pointMap[equalPointi];
154 
155  if (verbose)
156  {
157  Pout<< "Foam::mergePoints : Merging points "
158  << pointi << " and " << equalPointi
159  << " with coordinates:" << points[pointi]
160  << " and " << points[equalPointi]
161  << endl;
162  }
163  }
164  else
165  {
166  // Differs. Store new point.
167  pointMap[pointi] = newPointi++;
168  }
169  }
170 
171  return newPointi;
172 }
173 
174 
175 template<class Type>
177 (
178  const UList<Type>& points,
179  const scalar mergeTol,
180  const bool verbose,
181  labelList& pointMap,
182  List<Type>& newPoints,
183  const Type& origin
184 )
185 {
186  label nUnique = mergePoints
187  (
188  points,
189  mergeTol,
190  verbose,
191  pointMap,
192  origin
193  );
194 
195  newPoints.setSize(nUnique);
196  forAll(pointMap, pointi)
197  {
198  newPoints[pointMap[pointi]] = points[pointi];
199  }
200 
201  return (nUnique != points.size());
202 }
203 
204 
205 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Various functions to operate on Lists.
scalar y
List< label > labelList
A List of labels.
Definition: labelList.H:56
dimensioned< scalar > magSqr(const dimensioned< Type > &)
void setSize(const label)
Reset size of List.
Definition: List.C:281
vector point
Point is a vector.
Definition: point.H:41
label newPointi
Definition: readKivaGrid.H:501
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
label mergePoints(const UList< Type > &points, const scalar mergeTol, const bool verbose, labelList &pointMap, const Type &origin=Type::zero)
Sorts and merges points. All points closer than/equal mergeTol get merged.
dimensioned< scalar > mag(const dimensioned< Type > &)
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311