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