mapNearestMethod.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) 2013-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 "mapNearestMethod.H"
27 #include "pointIndexHit.H"
28 #include "indexedOctree.H"
29 #include "treeDataCell.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(mapNearestMethod, 0);
37  addToRunTimeSelectionTable(meshToMeshMethod, mapNearestMethod, components);
38 }
39 
40 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
41 
43 (
44  const labelList& srcCellIDs,
45  const boolList& mapFlag,
46  const label startSeedI,
47  label& srcSeedI,
48  label& tgtSeedI
49 ) const
50 {
51  const vectorField& srcCcs = src_.cellCentres();
52 
53  for (label i = startSeedI; i < srcCellIDs.size(); i++)
54  {
55  label srcI = srcCellIDs[i];
56 
57  if (mapFlag[srcI])
58  {
59  const point& srcCc = srcCcs[srcI];
60  pointIndexHit hit =
61  tgt_.cellTree().findNearest(srcCc, great);
62 
63  if (hit.hit())
64  {
65  srcSeedI = srcI;
66  tgtSeedI = hit.index();
67 
68  return true;
69  }
70  else
71  {
73  << "Unable to find nearest target cell"
74  << " for source cell " << srcI
75  << " with centre " << srcCc
76  << abort(FatalError);
77  }
78 
79  break;
80  }
81  }
82 
83  if (debug)
84  {
85  Pout<< "could not find starting seed" << endl;
86  }
87 
88  return false;
89 }
90 
91 
93 (
94  labelListList& srcToTgtCellAddr,
95  scalarListList& srcToTgtCellWght,
96  labelListList& tgtToSrcCellAddr,
97  scalarListList& tgtToSrcCellWght,
98  const label srcSeedI,
99  const label tgtSeedI,
100  const labelList& srcCellIDs,
101  boolList& mapFlag,
102  label& startSeedI
103 )
104 {
105  List<DynamicList<label>> srcToTgt(src_.nCells());
106  List<DynamicList<label>> tgtToSrc(tgt_.nCells());
107 
108  const scalarField& srcVc = src_.cellVolumes();
109  const scalarField& tgtVc = tgt_.cellVolumes();
110 
111  {
112  label srcCelli = srcSeedI;
113  label tgtCelli = tgtSeedI;
114 
115  do
116  {
117  // find nearest tgt cell
118  findNearestCell(src_, tgt_, srcCelli, tgtCelli);
119 
120  // store src/tgt cell pair
121  srcToTgt[srcCelli].append(tgtCelli);
122  tgtToSrc[tgtCelli].append(srcCelli);
123 
124  // mark source cell srcCelli and tgtCelli as matched
125  mapFlag[srcCelli] = false;
126 
127  // accumulate intersection volume
128  V_ += srcVc[srcCelli];
129 
130  // find new source cell
131  setNextNearestCells
132  (
133  startSeedI,
134  srcCelli,
135  tgtCelli,
136  mapFlag,
137  srcCellIDs
138  );
139  }
140  while (srcCelli >= 0);
141  }
142 
143  // for the case of multiple source cells per target cell, select the
144  // nearest source cell only and discard the others
145  const vectorField& srcCc = src_.cellCentres();
146  const vectorField& tgtCc = tgt_.cellCentres();
147 
148  forAll(tgtToSrc, targetCelli)
149  {
150  if (tgtToSrc[targetCelli].size() > 1)
151  {
152  const vector& tgtC = tgtCc[targetCelli];
153 
154  DynamicList<label>& srcCells = tgtToSrc[targetCelli];
155 
156  label srcCelli = srcCells[0];
157  scalar d = magSqr(tgtC - srcCc[srcCelli]);
158 
159  for (label i = 1; i < srcCells.size(); i++)
160  {
161  label srcI = srcCells[i];
162  scalar dNew = magSqr(tgtC - srcCc[srcI]);
163  if (dNew < d)
164  {
165  d = dNew;
166  srcCelli = srcI;
167  }
168  }
169 
170  srcCells.clear();
171  srcCells.append(srcCelli);
172  }
173  }
174 
175  // If there are more target cells than source cells, some target cells
176  // might not yet be mapped
177  forAll(tgtToSrc, tgtCelli)
178  {
179  if (tgtToSrc[tgtCelli].empty())
180  {
181  label srcCelli = findMappedSrcCell(tgtCelli, tgtToSrc);
182 
183  findNearestCell(tgt_, src_, tgtCelli, srcCelli);
184 
185  tgtToSrc[tgtCelli].append(srcCelli);
186  }
187  }
188 
189  // transfer addressing into persistent storage
190  forAll(srcToTgtCellAddr, i)
191  {
192  srcToTgtCellWght[i] = scalarList(srcToTgt[i].size(), srcVc[i]);
193  srcToTgtCellAddr[i].transfer(srcToTgt[i]);
194  }
195 
196  forAll(tgtToSrcCellAddr, i)
197  {
198  tgtToSrcCellWght[i] = scalarList(tgtToSrc[i].size(), tgtVc[i]);
199  tgtToSrcCellAddr[i].transfer(tgtToSrc[i]);
200  }
201 }
202 
203 
205 (
206  const polyMesh& mesh1,
207  const polyMesh& mesh2,
208  const label cell1,
209  label& cell2
210 ) const
211 {
212  const vectorField& Cc1 = mesh1.cellCentres();
213  const vectorField& Cc2 = mesh2.cellCentres();
214 
215  const vector& p1 = Cc1[cell1];
216 
217  DynamicList<label> cells2(10);
218  cells2.append(cell2);
219 
220  DynamicList<label> visitedCells(10);
221 
222  scalar d = great;
223 
224  do
225  {
226  label c2 = cells2.remove();
227  visitedCells.append(c2);
228 
229  scalar dTest = magSqr(Cc2[c2] - p1);
230  if (dTest < d)
231  {
232  cell2 = c2;
233  d = dTest;
234  appendNbrCells(cell2, mesh2, visitedCells, cells2);
235  }
236 
237  } while (cells2.size() > 0);
238 }
239 
240 
242 (
243  label& startSeedI,
244  label& srcCelli,
245  label& tgtCelli,
246  boolList& mapFlag,
247  const labelList& srcCellIDs
248 ) const
249 {
250  const labelList& srcNbr = src_.cellCells()[srcCelli];
251 
252  srcCelli = -1;
253  forAll(srcNbr, i)
254  {
255  label celli = srcNbr[i];
256  if (mapFlag[celli])
257  {
258  srcCelli = celli;
259  return;
260  }
261  }
262 
263  for (label i = startSeedI; i < srcCellIDs.size(); i++)
264  {
265  label celli = srcCellIDs[i];
266  if (mapFlag[celli])
267  {
268  startSeedI = i;
269  break;
270  }
271  }
272 
273  (void)findInitialSeeds
274  (
275  srcCellIDs,
276  mapFlag,
277  startSeedI,
278  srcCelli,
279  tgtCelli
280  );
281 }
282 
283 
285 (
286  const label tgtCelli,
287  const List<DynamicList<label>>& tgtToSrc
288 ) const
289 {
290  DynamicList<label> testCells(10);
291  DynamicList<label> visitedCells(10);
292 
293  testCells.append(tgtCelli);
294 
295  do
296  {
297  // search target tgtCelli neighbours for match with source cell
298  label tgtI = testCells.remove();
299 
300  if (findIndex(visitedCells, tgtI) == -1)
301  {
302  visitedCells.append(tgtI);
303 
304  if (tgtToSrc[tgtI].size())
305  {
306  return tgtToSrc[tgtI][0];
307  }
308  else
309  {
310  const labelList& nbrCells = tgt_.cellCells()[tgtI];
311 
312  forAll(nbrCells, i)
313  {
314  if (findIndex(visitedCells, nbrCells[i]) == -1)
315  {
316  testCells.append(nbrCells[i]);
317  }
318  }
319  }
320  }
321  } while (testCells.size());
322 
323  // did not find any match - should not be possible to get here!
324  return -1;
325 }
326 
327 
328 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
329 
331 (
332  const polyMesh& src,
333  const polyMesh& tgt
334 )
335 :
336  meshToMeshMethod(src, tgt)
337 {}
338 
339 
340 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
341 
343 {}
344 
345 
346 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
347 
349 (
350  labelListList& srcToTgtAddr,
351  scalarListList& srcToTgtWght,
352  labelListList& tgtToSrcAddr,
353  scalarListList& tgtToSrcWght
354 )
355 {
356  bool ok = initialise
357  (
358  srcToTgtAddr,
359  srcToTgtWght,
360  tgtToSrcAddr,
361  tgtToSrcWght
362  );
363 
364  if (!ok)
365  {
366  return;
367  }
368 
369  // (potentially) participating source mesh cells
370  const labelList srcCellIDs(maskCells());
371 
372  // list to keep track of whether src cell can be mapped
373  boolList mapFlag(src_.nCells(), false);
374  UIndirectList<bool>(mapFlag, srcCellIDs) = true;
375 
376  // find initial point in tgt mesh
377  label srcSeedI = -1;
378  label tgtSeedI = -1;
379  label startSeedI = 0;
380 
381  bool startWalk =
382  findInitialSeeds
383  (
384  srcCellIDs,
385  mapFlag,
386  startSeedI,
387  srcSeedI,
388  tgtSeedI
389  );
390 
391  if (startWalk)
392  {
393  calculateAddressing
394  (
395  srcToTgtAddr,
396  srcToTgtWght,
397  tgtToSrcAddr,
398  tgtToSrcWght,
399  srcSeedI,
400  tgtSeedI,
401  srcCellIDs,
402  mapFlag,
403  startSeedI
404  );
405  }
406  else
407  {
408  // if meshes are collocated, after inflating the source mesh bounding
409  // box tgt mesh cells may be transferred, but may still not overlap
410  // with the source mesh
411  return;
412  }
413 }
414 
415 
416 // ************************************************************************* //
virtual void calculateAddressing(labelListList &srcToTgtCellAddr, scalarListList &srcToTgtCellWght, labelListList &tgtToSrcCellAddr, scalarListList &tgtToSrcCellWght, const label srcSeedI, const label tgtSeedI, const labelList &srcCellIDs, boolList &mapFlag, label &startSeedI)
Calculate the mesh-to-mesh addressing and weights.
#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
virtual void setNextNearestCells(label &startSeedI, label &srcCelli, label &tgtCelli, boolList &mapFlag, const labelList &srcCellIDs) const
Set the next cells for the marching front algorithm.
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
virtual void calculate(labelListList &srcToTgtAddr, scalarListList &srcToTgtWght, labelListList &tgtToTgtAddr, scalarListList &tgtToTgtWght)
Calculate addressing and weights.
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
virtual void findNearestCell(const polyMesh &mesh1, const polyMesh &mesh2, const label cell1, label &cell2) const
Find the nearest cell on mesh2 for cell1 on mesh1.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
Macros for easy insertion into run-time selection tables.
const dimensionedScalar c2
Second radiation constant: default SI units: [m K].
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
bool hit() const
Is there a hit.
const vectorField & cellCentres() const
errorManip< error > abort(error &err)
Definition: errorManip.H:131
dimensioned< scalar > magSqr(const dimensioned< Type > &)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
Base class for mesh-to-mesh calculation methods.
T remove()
Remove and return the top element.
Definition: DynamicListI.H:351
virtual ~mapNearestMethod()
Destructor.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
A List with indirect addressing.
Definition: fvMatrix.H:106
virtual bool findInitialSeeds(const labelList &srcCellIDs, const boolList &mapFlag, const label startSeedI, label &srcSeedI, label &tgtSeedI) const
Find indices of overlapping cells in src and tgt meshes - returns.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:236
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
label index() const
Return index.
virtual label findMappedSrcCell(const label tgtCelli, const List< DynamicList< label >> &tgtToSrc) const
Find a source cell mapped to target cell tgtCelli.
Namespace for OpenFOAM.
mapNearestMethod(const polyMesh &src, const polyMesh &tgt)
Construct from source and target meshes.