nearestCellsToCells.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-2023 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 "nearestCellsToCells.H"
27 #include "pointIndexHit.H"
28 #include "indexedOctree.H"
29 #include "treeDataCell.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace cellsToCellss
37 {
40 }
41 }
42 
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
46 bool Foam::cellsToCellss::nearest::findInitialSeeds
47 (
48  const polyMesh& srcMesh,
49  const polyMesh& tgtMesh,
50  const labelList& srcCellIDs,
51  const boolList& mapFlag,
52  const label startSeedI,
53  label& srcSeedI,
54  label& tgtSeedI
55 ) const
56 {
57  const vectorField& srcCcs = srcMesh.cellCentres();
58 
59  for (label i = startSeedI; i < srcCellIDs.size(); i++)
60  {
61  label srcI = srcCellIDs[i];
62 
63  if (mapFlag[srcI])
64  {
65  const point& srcCc = srcCcs[srcI];
66  pointIndexHit hit =
67  tgtMesh.cellTree().findNearest(srcCc, great);
68 
69  if (hit.hit())
70  {
71  srcSeedI = srcI;
72  tgtSeedI = hit.index();
73 
74  return true;
75  }
76  else
77  {
79  << "Unable to find nearest target cell"
80  << " for source cell " << srcI
81  << " with centre " << srcCc
82  << abort(FatalError);
83  }
84 
85  break;
86  }
87  }
88 
89  if (debug)
90  {
91  Pout<< "could not find starting seed" << endl;
92  }
93 
94  return false;
95 }
96 
97 
98 Foam::scalar Foam::cellsToCellss::nearest::calculateAddressing
99 (
100  const polyMesh& srcMesh,
101  const polyMesh& tgtMesh,
102  labelListList& srcToTgtCellAddr,
103  scalarListList& srcToTgtCellWght,
104  labelListList& tgtToSrcCellAddr,
105  scalarListList& tgtToSrcCellWght,
106  const label srcSeedI,
107  const label tgtSeedI,
108  const labelList& srcCellIDs,
109  boolList& mapFlag,
110  label& startSeedI
111 )
112 {
113  scalar V = 0;
114 
115  List<DynamicList<label>> srcToTgt(srcMesh.nCells());
116  List<DynamicList<label>> tgtToSrc(tgtMesh.nCells());
117 
118  const scalarField& srcVc = srcMesh.cellVolumes();
119  const scalarField& tgtVc = tgtMesh.cellVolumes();
120 
121  {
122  label srcCelli = srcSeedI;
123  label tgtCelli = tgtSeedI;
124 
125  do
126  {
127  // find nearest tgt cell
128  findNearestCell(srcMesh, tgtMesh, 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  srcMesh,
144  tgtMesh,
145  startSeedI,
146  srcCelli,
147  tgtCelli,
148  mapFlag,
149  srcCellIDs
150  );
151  }
152  while (srcCelli >= 0);
153  }
154 
155  // for the case of multiple source cells per target cell, select the
156  // nearest source cell only and discard the others
157  const vectorField& srcCc = srcMesh.cellCentres();
158  const vectorField& tgtCc = tgtMesh.cellCentres();
159 
160  forAll(tgtToSrc, targetCelli)
161  {
162  if (tgtToSrc[targetCelli].size() > 1)
163  {
164  const vector& tgtC = tgtCc[targetCelli];
165 
166  DynamicList<label>& srcCells = tgtToSrc[targetCelli];
167 
168  label srcCelli = srcCells[0];
169  scalar d = magSqr(tgtC - srcCc[srcCelli]);
170 
171  for (label i = 1; i < srcCells.size(); i++)
172  {
173  label srcI = srcCells[i];
174  scalar dNew = magSqr(tgtC - srcCc[srcI]);
175  if (dNew < d)
176  {
177  d = dNew;
178  srcCelli = srcI;
179  }
180  }
181 
182  srcCells.clear();
183  srcCells.append(srcCelli);
184  }
185  }
186 
187  // If there are more target cells than source cells, some target cells
188  // might not yet be mapped
189  forAll(tgtToSrc, tgtCelli)
190  {
191  if (tgtToSrc[tgtCelli].empty())
192  {
193  label srcCelli =
194  findMappedSrcCell(srcMesh, tgtMesh, tgtCelli, tgtToSrc);
195 
196  findNearestCell(tgtMesh, srcMesh, tgtCelli, srcCelli);
197 
198  tgtToSrc[tgtCelli].append(srcCelli);
199  }
200  }
201 
202  // transfer addressing into persistent storage
203  forAll(srcToTgtCellAddr, i)
204  {
205  srcToTgtCellWght[i] = scalarList(srcToTgt[i].size(), srcVc[i]);
206  srcToTgtCellAddr[i].transfer(srcToTgt[i]);
207  }
208 
209  forAll(tgtToSrcCellAddr, i)
210  {
211  tgtToSrcCellWght[i] = scalarList(tgtToSrc[i].size(), tgtVc[i]);
212  tgtToSrcCellAddr[i].transfer(tgtToSrc[i]);
213  }
214 
215  return V;
216 }
217 
218 
219 void Foam::cellsToCellss::nearest::findNearestCell
220 (
221  const polyMesh& mesh1,
222  const polyMesh& mesh2,
223  const label cell1,
224  label& cell2
225 ) const
226 {
227  const vectorField& Cc1 = mesh1.cellCentres();
228  const vectorField& Cc2 = mesh2.cellCentres();
229 
230  const vector& p1 = Cc1[cell1];
231 
232  DynamicList<label> cells2(10);
233  cells2.append(cell2);
234 
235  DynamicList<label> visitedCells(10);
236 
237  scalar d = great;
238 
239  do
240  {
241  label c2 = cells2.remove();
242  visitedCells.append(c2);
243 
244  scalar dTest = magSqr(Cc2[c2] - p1);
245  if (dTest < d)
246  {
247  cell2 = c2;
248  d = dTest;
249  appendNbrCells(cell2, mesh2, visitedCells, cells2);
250  }
251 
252  } while (cells2.size() > 0);
253 }
254 
255 
256 void Foam::cellsToCellss::nearest::setNextNearestCells
257 (
258  const polyMesh& srcMesh,
259  const polyMesh& tgtMesh,
260  label& startSeedI,
261  label& srcCelli,
262  label& tgtCelli,
263  boolList& mapFlag,
264  const labelList& srcCellIDs
265 ) const
266 {
267  const labelList& srcNbr = srcMesh.cellCells()[srcCelli];
268 
269  srcCelli = -1;
270  forAll(srcNbr, i)
271  {
272  label celli = srcNbr[i];
273  if (mapFlag[celli])
274  {
275  srcCelli = celli;
276  return;
277  }
278  }
279 
280  for (label i = startSeedI; i < srcCellIDs.size(); i++)
281  {
282  label celli = srcCellIDs[i];
283  if (mapFlag[celli])
284  {
285  startSeedI = i;
286  break;
287  }
288  }
289 
290  findInitialSeeds
291  (
292  srcMesh,
293  tgtMesh,
294  srcCellIDs,
295  mapFlag,
296  startSeedI,
297  srcCelli,
298  tgtCelli
299  );
300 }
301 
302 
303 Foam::label Foam::cellsToCellss::nearest::findMappedSrcCell
304 (
305  const polyMesh& srcMesh,
306  const polyMesh& tgtMesh,
307  const label tgtCelli,
308  const List<DynamicList<label>>& tgtToSrc
309 ) const
310 {
311  DynamicList<label> testCells(10);
312  DynamicList<label> visitedCells(10);
313 
314  testCells.append(tgtCelli);
315 
316  do
317  {
318  // search target tgtCelli neighbours for match with source cell
319  label tgtI = testCells.remove();
320 
321  if (findIndex(visitedCells, tgtI) == -1)
322  {
323  visitedCells.append(tgtI);
324 
325  if (tgtToSrc[tgtI].size())
326  {
327  return tgtToSrc[tgtI][0];
328  }
329  else
330  {
331  const labelList& nbrCells = tgtMesh.cellCells()[tgtI];
332 
333  forAll(nbrCells, i)
334  {
335  if (findIndex(visitedCells, nbrCells[i]) == -1)
336  {
337  testCells.append(nbrCells[i]);
338  }
339  }
340  }
341  }
342  } while (testCells.size());
343 
344  // did not find any match - should not be possible to get here!
345  return -1;
346 }
347 
348 
349 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
350 
352 (
353  const polyMesh& srcMesh,
354  const polyMesh& tgtMesh
355 )
356 {
357  initialise(srcMesh, tgtMesh);
358 
359  // Determine (potentially) participating source mesh cells
360  const labelList srcCellIDs(maskCells(srcMesh, tgtMesh));
361 
362  // Initialise list to keep track of whether src cell can be mapped
363  boolList mapFlag(srcMesh.nCells(), false);
364  UIndirectList<bool>(mapFlag, srcCellIDs) = true;
365 
366  // Find initial point in tgt mesh
367  label srcSeedI = -1;
368  label tgtSeedI = -1;
369  label startSeedI = 0;
370  bool startWalk =
371  findInitialSeeds
372  (
373  srcMesh,
374  tgtMesh,
375  srcCellIDs,
376  mapFlag,
377  startSeedI,
378  srcSeedI,
379  tgtSeedI
380  );
381 
382  if (startWalk)
383  {
384  return
385  calculateAddressing
386  (
387  srcMesh,
388  tgtMesh,
389  srcLocalTgtCells_,
390  srcWeights_,
391  tgtLocalSrcCells_,
392  tgtWeights_,
393  srcSeedI,
394  tgtSeedI,
395  srcCellIDs,
396  mapFlag,
397  startSeedI
398  );
399  }
400  else
401  {
402  return 0;
403  }
404 }
405 
406 
408 (
409  const polyMesh& srcMesh,
410  labelListList& srcToTgtAddr,
411  scalarListList& srcToTgtWght
412 ) const
413 {
414  forAll(srcToTgtWght, srcCelli)
415  {
416  if (srcToTgtWght[srcCelli].size() > 1)
417  {
418  srcToTgtAddr[srcCelli].resize(1);
419  srcToTgtWght[srcCelli].resize(1);
420  srcToTgtWght[srcCelli][0] = 1;
421  }
422  }
423 }
424 
425 
426 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
427 
429 :
430  cellsToCells()
431 {}
432 
433 
434 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
435 
437 {}
438 
439 
440 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
void resize(const label)
Alias for setSize(const label)
Definition: ListI.H:138
A List with indirect addressing.
Definition: UIndirectList.H:60
Class to calculate interpolative addressing and weights between the cells of two overlapping meshes.
Definition: cellsToCells.H:56
Nearest cells-to-cells interpolation class.
virtual scalar calculate(const polyMesh &srcMesh, const polyMesh &tgtMesh)
Calculate the addressing and weights.
virtual void normalise(const polyMesh &mesh, labelListList &localOtherCells, scalarListList &weights) const
Normalise the weights for a given mesh.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
label nCells() const
A class for handling words, derived from string.
Definition: word.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
defineTypeNameAndDebug(intersection, 0)
addToRunTimeSelectionTable(cellsToCells, intersection, word)
const dimensionedScalar c2
Second radiation constant: default SI units: [m K].
Namespace for OpenFOAM.
List< scalarList > scalarListList
Definition: scalarList.H:51
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
errorManip< error > abort(error &err)
Definition: errorManip.H:131
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
vector point
Point is a vector.
Definition: point.H:41
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
Field< vector > vectorField
Specialisation of Field<T> for vector.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
error FatalError
dimensioned< scalar > magSqr(const dimensioned< Type > &)