matchPoints.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "matchPoints.H"
30 #include "SortableList.H"
31 #include "ListOps.H"
32 
33 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
34 
36 (
37  const UList<point>& pts0,
38  const UList<point>& pts1,
39  const UList<scalar>& matchDistances,
40  const bool verbose,
41  labelList& from0To1,
42  const point& origin
43 )
44 {
45  from0To1.setSize(pts0.size());
46  from0To1 = -1;
47 
48  bool fullMatch = true;
49 
50  point compareOrigin = origin;
51 
52  if (origin == point::max)
53  {
54  if (pts1.size())
55  {
56  compareOrigin = sum(pts1)/pts1.size();
57  }
58  else
59  {
60  // Unusable, but avoid comparison with VGREAT!
61  compareOrigin = point::zero;
62  }
63  }
64 
65  SortableList<scalar> pts0MagSqr(magSqr(pts0 - compareOrigin));
66 
67  SortableList<scalar> pts1MagSqr(magSqr(pts1 - compareOrigin));
68 
69  forAll(pts0MagSqr, i)
70  {
71  scalar dist0 = pts0MagSqr[i];
72 
73  label face0I = pts0MagSqr.indices()[i];
74 
75  scalar matchDist = matchDistances[face0I];
76 
77  label startI = findLower(pts1MagSqr, 0.99999*dist0 - 2*matchDist);
78 
79  if (startI == -1)
80  {
81  startI = 0;
82  }
83 
84 
85  // Go through range of equal mag and find nearest vector.
86  scalar minDistSqr = VGREAT;
87  label minFacei = -1;
88 
89  for
90  (
91  label j = startI;
92  (
93  (j < pts1MagSqr.size())
94  && (pts1MagSqr[j] < 1.00001*dist0 + 2*matchDist)
95  );
96  j++
97  )
98  {
99  label facei = pts1MagSqr.indices()[j];
100  // Compare actual vectors
101  scalar distSqr = magSqr(pts0[face0I] - pts1[facei]);
102 
103  if (distSqr <= sqr(matchDist) && distSqr < minDistSqr)
104  {
105  minDistSqr = distSqr;
106  minFacei = facei;
107  }
108  }
109 
110  if (minFacei == -1)
111  {
112  fullMatch = false;
113 
114  if (verbose)
115  {
116  Pout<< "Cannot find point in pts1 matching point " << face0I
117  << " coord:" << pts0[face0I]
118  << " in pts0 when using tolerance " << matchDist << endl;
119 
120  // Go through range of equal mag and find equal vector.
121  Pout<< "Searching started from:" << startI << " in pts1"
122  << endl;
123  for
124  (
125  label j = startI;
126  (
127  (j < pts1MagSqr.size())
128  && (pts1MagSqr[j] < 1.00001*dist0 + 2*matchDist)
129  );
130  j++
131  )
132  {
133  label facei = pts1MagSqr.indices()[j];
134 
135  Pout<< " Compared coord: " << pts1[facei]
136  << " at index " << j
137  << " with difference to point "
138  << mag(pts1[facei] - pts0[face0I]) << endl;
139  }
140  }
141  }
142 
143  from0To1[face0I] = minFacei;
144  }
145 
146  return fullMatch;
147 }
148 
149 
151 (
152  const UList<point>& pts0,
153  const UList<point>& pts1,
154  const UList<point>& pts0Dir,
155  const UList<point>& pts1Dir,
156  const UList<scalar>& matchDistances,
157  const bool verbose,
158  labelList& from0To1,
159  const point& origin
160 )
161 {
162  from0To1.setSize(pts0.size());
163  from0To1 = -1;
164 
165  bool fullMatch = true;
166 
167  point compareOrigin = origin;
168 
169  if (origin == point::max)
170  {
171  if (pts1.size())
172  {
173  compareOrigin = sum(pts1)/pts1.size();
174  }
175  else
176  {
177  // Unusable, but avoid comparison with VGREAT!
178  compareOrigin = point::zero;
179  }
180  }
181 
182  SortableList<scalar> pts0MagSqr(magSqr(pts0 - compareOrigin));
183 
184  SortableList<scalar> pts1MagSqr(magSqr(pts1 - compareOrigin));
185 
186  forAll(pts0MagSqr, i)
187  {
188  scalar dist0 = pts0MagSqr[i];
189 
190  label face0I = pts0MagSqr.indices()[i];
191 
192  scalar matchDist = matchDistances[face0I];
193 
194  label startI = findLower(pts1MagSqr, 0.99999*dist0 - 2*matchDist);
195 
196  if (startI == -1)
197  {
198  startI = 0;
199  }
200 
201  // Go through range of equal mag and find nearest vector.
202  scalar minDistSqr = VGREAT;
203  label minFacei = -1;
204 
205  // Valid candidate points should have opposite normal
206  const scalar minDistNorm = 0;
207 
208  for
209  (
210  label j = startI;
211  (
212  (j < pts1MagSqr.size())
213  && (pts1MagSqr[j] < 1.00001*dist0 + 2*matchDist)
214  );
215  j++
216  )
217  {
218  label facei = pts1MagSqr.indices()[j];
219  // Compare actual vectors
220  scalar distSqr = magSqr(pts0[face0I] - pts1[facei]);
221 
222  scalar distNorm = (pts0Dir[face0I] & pts1Dir[facei]);
223 
224  if
225  (
226  magSqr(pts0Dir[face0I]) < sqr(SMALL)
227  && magSqr(pts1Dir[facei]) < sqr(SMALL)
228  )
229  {
230  distNorm = -1;
231  }
232 
233  if (distSqr <= sqr(matchDist) && distSqr < minDistSqr)
234  {
235  // Check that the normals point in equal and opposite directions
236  if (distNorm < minDistNorm)
237  {
238  minDistSqr = distSqr;
239  minFacei = facei;
240  }
241  }
242  }
243 
244  if (minFacei == -1)
245  {
246  fullMatch = false;
247 
248  if (verbose)
249  {
250  Pout<< "Cannot find point in pts1 matching point " << face0I
251  << " coord:" << pts0[face0I]
252  << " in pts0 when using tolerance " << matchDist << endl;
253 
254  // Go through range of equal mag and find equal vector.
255  Pout<< "Searching started from:" << startI << " in pts1"
256  << endl;
257  for
258  (
259  label j = startI;
260  (
261  (j < pts1MagSqr.size())
262  && (pts1MagSqr[j] < 1.00001*dist0 + 2*matchDist)
263  );
264  j++
265  )
266  {
267  label facei = pts1MagSqr.indices()[j];
268 
269  Pout<< " Compared coord: " << pts1[facei]
270  << " at index " << j
271  << " with difference to point "
272  << mag(pts1[facei] - pts0[face0I]) << endl;
273  }
274  }
275  }
276 
277  from0To1[face0I] = minFacei;
278  }
279 
280  return fullMatch;
281 }
282 
283 
284 // ************************************************************************* //
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
SortableList.H
matchPoints.H
Determine correspondence between points. See below.
Foam::findLower
label findLower(const ListType &input, const T &val, const label start, const ComparePredicate &comp)
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::SortableList
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: List.H:60
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::Vector< scalar >
Foam::List< label >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::UList
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:103
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
Foam::SortableList::indices
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
Definition: SortableList.H:108
ListOps.H
Various functions to operate on Lists.
Foam::UList::size
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Foam::matchPoints
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:36