matchPoints.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 "matchPoints.H"
27 #include "SortableList.H"
28 #include "ListOps.H"
29 
30 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
31 
33 (
34  const UList<point>& pts0,
35  const UList<point>& pts1,
36  const UList<scalar>& matchDistances,
37  const bool verbose,
38  labelList& from0To1,
39  const point& origin
40 )
41 {
42  from0To1.setSize(pts0.size());
43  from0To1 = -1;
44 
45  bool fullMatch = true;
46 
47  point compareOrigin = origin;
48 
49  if (origin == point(vGreat, vGreat, vGreat))
50  {
51  if (pts1.size())
52  {
53  compareOrigin = sum(pts1)/pts1.size();
54  }
55  }
56 
57  SortableList<scalar> pts0MagSqr(magSqr(pts0 - compareOrigin));
58 
59  SortableList<scalar> pts1MagSqr(magSqr(pts1 - compareOrigin));
60 
61  forAll(pts0MagSqr, i)
62  {
63  scalar dist0 = pts0MagSqr[i];
64 
65  label face0I = pts0MagSqr.indices()[i];
66 
67  scalar matchDist = matchDistances[face0I];
68 
69  label startI = findLower(pts1MagSqr, 0.99999*dist0 - 2*matchDist);
70 
71  if (startI == -1)
72  {
73  startI = 0;
74  }
75 
76 
77  // Go through range of equal mag and find nearest vector.
78  scalar minDistSqr = vGreat;
79  label minFacei = -1;
80 
81  for
82  (
83  label j = startI;
84  (
85  (j < pts1MagSqr.size())
86  && (pts1MagSqr[j] < 1.00001*dist0 + 2*matchDist)
87  );
88  j++
89  )
90  {
91  label facei = pts1MagSqr.indices()[j];
92  // Compare actual vectors
93  scalar distSqr = magSqr(pts0[face0I] - pts1[facei]);
94 
95  if (distSqr <= sqr(matchDist) && distSqr < minDistSqr)
96  {
97  minDistSqr = distSqr;
98  minFacei = facei;
99  }
100  }
101 
102  if (minFacei == -1)
103  {
104  fullMatch = false;
105 
106  if (verbose)
107  {
108  Pout<< "Cannot find point in pts1 matching point " << face0I
109  << " coord:" << pts0[face0I]
110  << " in pts0 when using tolerance " << matchDist << endl;
111 
112  // Go through range of equal mag and find equal vector.
113  Pout<< "Searching started from:" << startI << " in pts1"
114  << endl;
115  for
116  (
117  label j = startI;
118  (
119  (j < pts1MagSqr.size())
120  && (pts1MagSqr[j] < 1.00001*dist0 + 2*matchDist)
121  );
122  j++
123  )
124  {
125  label facei = pts1MagSqr.indices()[j];
126 
127  Pout<< " Compared coord: " << pts1[facei]
128  << " at index " << j
129  << " with difference to point "
130  << mag(pts1[facei] - pts0[face0I]) << endl;
131  }
132  }
133  }
134 
135  from0To1[face0I] = minFacei;
136  }
137 
138  return fullMatch;
139 }
140 
141 
143 (
144  const UList<point>& pts0,
145  const UList<point>& pts1,
146  const UList<point>& pts0Dir,
147  const UList<point>& pts1Dir,
148  const UList<scalar>& matchDistances,
149  const bool verbose,
150  labelList& from0To1,
151  const point& origin
152 )
153 {
154  from0To1.setSize(pts0.size());
155  from0To1 = -1;
156 
157  bool fullMatch = true;
158 
159  point compareOrigin = origin;
160 
161  if (origin == point(vGreat, vGreat, vGreat))
162  {
163  if (pts1.size())
164  {
165  compareOrigin = sum(pts1)/pts1.size();
166  }
167  }
168 
169  SortableList<scalar> pts0MagSqr(magSqr(pts0 - compareOrigin));
170 
171  SortableList<scalar> pts1MagSqr(magSqr(pts1 - compareOrigin));
172 
173  forAll(pts0MagSqr, i)
174  {
175  scalar dist0 = pts0MagSqr[i];
176 
177  label face0I = pts0MagSqr.indices()[i];
178 
179  scalar matchDist = matchDistances[face0I];
180 
181  label startI = findLower(pts1MagSqr, 0.99999*dist0 - 2*matchDist);
182 
183  if (startI == -1)
184  {
185  startI = 0;
186  }
187 
188  // Go through range of equal mag and find nearest vector.
189  scalar minDistSqr = vGreat;
190  scalar minDistNorm = 0;
191  label minFacei = -1;
192 
193  for
194  (
195  label j = startI;
196  (
197  (j < pts1MagSqr.size())
198  && (pts1MagSqr[j] < 1.00001*dist0 + 2*matchDist)
199  );
200  j++
201  )
202  {
203  label facei = pts1MagSqr.indices()[j];
204  // Compare actual vectors
205  scalar distSqr = magSqr(pts0[face0I] - pts1[facei]);
206 
207  scalar distNorm = (pts0Dir[face0I] & pts1Dir[facei]);
208 
209  if
210  (
211  magSqr(pts0Dir[face0I]) < sqr(small)
212  && magSqr(pts1Dir[facei]) < sqr(small)
213  )
214  {
215  distNorm = -1;
216  }
217 
218  if (distSqr <= sqr(matchDist) && distSqr < minDistSqr)
219  {
220  // Check that the normals point in equal and opposite directions
221  if (distNorm < minDistNorm)
222  {
223  minDistNorm = distNorm;
224  minDistSqr = distSqr;
225  minFacei = facei;
226  }
227  }
228  }
229 
230  if (minFacei == -1)
231  {
232  fullMatch = false;
233 
234  if (verbose)
235  {
236  Pout<< "Cannot find point in pts1 matching point " << face0I
237  << " coord:" << pts0[face0I]
238  << " in pts0 when using tolerance " << matchDist << endl;
239 
240  // Go through range of equal mag and find equal vector.
241  Pout<< "Searching started from:" << startI << " in pts1"
242  << endl;
243  for
244  (
245  label j = startI;
246  (
247  (j < pts1MagSqr.size())
248  && (pts1MagSqr[j] < 1.00001*dist0 + 2*matchDist)
249  );
250  j++
251  )
252  {
253  label facei = pts1MagSqr.indices()[j];
254 
255  Pout<< " Compared coord: " << pts1[facei]
256  << " at index " << j
257  << " with difference to point "
258  << mag(pts1[facei] - pts0[face0I]) << endl;
259  }
260  }
261  }
262 
263  from0To1[face0I] = minFacei;
264  }
265 
266  return fullMatch;
267 }
268 
269 
270 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A list that is sorted upon construction or when explicitly requested with the sort() method...
Definition: List.H:80
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Determine correspondence between points. See below.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Various functions to operate on Lists.
bool matchPoints(const UList< point > &pts0, const UList< point > &pts1, const UList< scalar > &matchDistance, const bool verbose, labelList &from0To1, const point &origin=point::zero)
Determine correspondence between pointFields. Gets passed.
Definition: matchPoints.C:33
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
Definition: SortableList.H:96
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
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
label findLower(const ListType &, typename ListType::const_reference, const label stary, const BinaryOp &bop)
Find last element < given value in sorted list and return index,.
dimensioned< scalar > mag(const dimensioned< Type > &)
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311