intersectionCellsToCells.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 
27 #include "indexedOctree.H"
28 #include "treeDataCell.H"
29 #include "tetOverlapVolume.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace cellsToCellss
37 {
40 }
41 }
42 
43 
44 const Foam::scalar Foam::cellsToCellss::intersection::tolerance_ = 1e-6;
45 
46 
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 
49 bool Foam::cellsToCellss::intersection::intersect
50 (
51  const polyMesh& srcMesh,
52  const polyMesh& tgtMesh,
53  const label srcCelli,
54  const label tgtCelli
55 ) const
56 {
57  return
58  tetOverlapVolume().cellCellOverlapMinDecomp
59  (
60  srcMesh,
61  srcCelli,
62  tgtMesh,
63  tgtCelli,
64  treeBoundBox(tgtMesh.points(), tgtMesh.cellPoints()[tgtCelli]),
65  tolerance_*srcMesh.cellVolumes()[srcCelli]
66  );
67 }
68 
69 
70 Foam::scalar Foam::cellsToCellss::intersection::interVol
71 (
72  const polyMesh& srcMesh,
73  const polyMesh& tgtMesh,
74  const label srcCelli,
75  const label tgtCelli
76 ) const
77 {
78  return
79  tetOverlapVolume().cellCellOverlapVolumeMinDecomp
80  (
81  srcMesh,
82  srcCelli,
83  tgtMesh,
84  tgtCelli,
85  treeBoundBox(tgtMesh.points(), tgtMesh.cellPoints()[tgtCelli])
86  );
87 }
88 
89 
90 bool Foam::cellsToCellss::intersection::findInitialSeeds
91 (
92  const polyMesh& srcMesh,
93  const polyMesh& tgtMesh,
94  const labelList& srcCellIDs,
95  const boolList& mapFlag,
96  const label startSeedI,
97  label& srcSeedI,
98  label& tgtSeedI
99 ) const
100 {
101  const cellList& srcCells = srcMesh.cells();
102  const faceList& srcFaces = srcMesh.faces();
103  const pointField& srcPts = srcMesh.points();
104 
105  for (label i = startSeedI; i < srcCellIDs.size(); i++)
106  {
107  const label srcI = srcCellIDs[i];
108 
109  if (mapFlag[srcI])
110  {
111  const labelList tgtIDs
112  (
113  tgtMesh.cellTree().findBox
114  (
115  treeBoundBox(srcCells[srcI].bb(srcPts, srcFaces))
116  )
117  );
118 
119  forAll(tgtIDs, j)
120  {
121  const label tgtI = tgtIDs[j];
122 
123  if (intersect(srcMesh, tgtMesh, srcI, tgtI))
124  {
125  srcSeedI = srcI;
126  tgtSeedI = tgtI;
127 
128  return true;
129  }
130  }
131  }
132  }
133 
134  if (debug)
135  {
136  Pout<< "could not find starting seed" << endl;
137  }
138 
139  return false;
140 }
141 
142 
143 Foam::scalar Foam::cellsToCellss::intersection::calculateAddressing
144 (
145  const polyMesh& srcMesh,
146  const polyMesh& tgtMesh,
147  labelListList& srcToTgtCellAddr,
148  scalarListList& srcToTgtCellWght,
149  labelListList& tgtToSrcCellAddr,
150  scalarListList& tgtToSrcCellWght,
151  const label srcSeedI,
152  const label tgtSeedI,
153  const labelList& srcCellIDs,
154  boolList& mapFlag,
155  label& startSeedI
156 )
157 {
158  scalar V = 0;
159 
160  label srcCelli = srcSeedI;
161  label tgtCelli = tgtSeedI;
162 
163  List<DynamicList<label>> srcToTgtAddr(srcMesh.nCells());
164  List<DynamicList<scalar>> srcToTgtWght(srcMesh.nCells());
165 
166  List<DynamicList<label>> tgtToSrcAddr(tgtMesh.nCells());
167  List<DynamicList<scalar>> tgtToSrcWght(tgtMesh.nCells());
168 
169  // list of tgt cell neighbour cells
170  DynamicList<label> nbrTgtCells(10);
171 
172  // list of tgt cells currently visited for srcCelli to avoid multiple hits
173  DynamicList<label> visitedTgtCells(10);
174 
175  // list to keep track of tgt cells used to seed src cells
176  labelList seedCells(srcMesh.nCells(), -1);
177  seedCells[srcCelli] = tgtCelli;
178 
179  const scalarField& srcVol = srcMesh.cellVolumes();
180 
181  do
182  {
183  nbrTgtCells.clear();
184  visitedTgtCells.clear();
185 
186  // append initial target cell and neighbours
187  nbrTgtCells.append(tgtCelli);
188  appendNbrCells(tgtCelli, tgtMesh, visitedTgtCells, nbrTgtCells);
189 
190  do
191  {
192  tgtCelli = nbrTgtCells.remove();
193  visitedTgtCells.append(tgtCelli);
194 
195  scalar vol = interVol(srcMesh, tgtMesh, srcCelli, tgtCelli);
196 
197  // accumulate addressing and weights for valid intersection
198  if (vol/srcVol[srcCelli] > tolerance_)
199  {
200  // store src/tgt cell pair
201  srcToTgtAddr[srcCelli].append(tgtCelli);
202  srcToTgtWght[srcCelli].append(vol);
203 
204  tgtToSrcAddr[tgtCelli].append(srcCelli);
205  tgtToSrcWght[tgtCelli].append(vol);
206 
207  appendNbrCells(tgtCelli, tgtMesh, visitedTgtCells, nbrTgtCells);
208 
209  // accumulate intersection volume
210  V += vol;
211  }
212  }
213  while (!nbrTgtCells.empty());
214 
215  mapFlag[srcCelli] = false;
216 
217  // find new source seed cell
218  setNextCells
219  (
220  srcMesh,
221  tgtMesh,
222  startSeedI,
223  srcCelli,
224  tgtCelli,
225  srcCellIDs,
226  mapFlag,
227  visitedTgtCells,
228  seedCells
229  );
230  }
231  while (srcCelli != -1);
232 
233  // transfer addressing into persistent storage
234  forAll(srcToTgtCellAddr, i)
235  {
236  srcToTgtCellAddr[i].transfer(srcToTgtAddr[i]);
237  srcToTgtCellWght[i].transfer(srcToTgtWght[i]);
238  }
239 
240  forAll(tgtToSrcCellAddr, i)
241  {
242  tgtToSrcCellAddr[i].transfer(tgtToSrcAddr[i]);
243  tgtToSrcCellWght[i].transfer(tgtToSrcWght[i]);
244  }
245 
246  return V;
247 }
248 
249 
250 void Foam::cellsToCellss::intersection::setNextCells
251 (
252  const polyMesh& srcMesh,
253  const polyMesh& tgtMesh,
254  label& startSeedI,
255  label& srcCelli,
256  label& tgtCelli,
257  const labelList& srcCellIDs,
258  const boolList& mapFlag,
259  const DynamicList<label>& visitedCells,
260  labelList& seedCells
261 ) const
262 {
263  const labelList& srcNbrCells = srcMesh.cellCells()[srcCelli];
264 
265  // set possible seeds for later use by querying all src cell neighbours
266  // with all visited target cells
267  bool valuesSet = false;
268  forAll(srcNbrCells, i)
269  {
270  label cellS = srcNbrCells[i];
271 
272  if (mapFlag[cellS] && seedCells[cellS] == -1)
273  {
274  forAll(visitedCells, j)
275  {
276  label cellT = visitedCells[j];
277 
278  if (intersect(srcMesh, tgtMesh, cellS, cellT))
279  {
280  seedCells[cellS] = cellT;
281 
282  if (!valuesSet)
283  {
284  srcCelli = cellS;
285  tgtCelli = cellT;
286  valuesSet = true;
287  }
288  }
289  }
290  }
291  }
292 
293  // set next src and tgt cells if not set above
294  if (valuesSet)
295  {
296  return;
297  }
298  else
299  {
300  // try to use existing seed
301  bool foundNextSeed = false;
302  for (label i = startSeedI; i < srcCellIDs.size(); i++)
303  {
304  label cellS = srcCellIDs[i];
305 
306  if (mapFlag[cellS])
307  {
308  if (!foundNextSeed)
309  {
310  startSeedI = i;
311  foundNextSeed = true;
312  }
313 
314  if (seedCells[cellS] != -1)
315  {
316  srcCelli = cellS;
317  tgtCelli = seedCells[cellS];
318 
319  return;
320  }
321  }
322  }
323 
324  // perform new search to find match
325  if (debug)
326  {
327  Pout<< "Advancing front stalled: searching for new "
328  << "target cell" << endl;
329  }
330 
331  bool restart =
332  findInitialSeeds
333  (
334  srcMesh,
335  tgtMesh,
336  srcCellIDs,
337  mapFlag,
338  startSeedI,
339  srcCelli,
340  tgtCelli
341  );
342 
343  if (restart)
344  {
345  // successfully found new starting seed-pair
346  return;
347  }
348  }
349 
350  // if we have got to here, there are no more src/tgt cell intersections
351  srcCelli = -1;
352  tgtCelli = -1;
353 }
354 
355 
356 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
357 
359 (
360  const polyMesh& srcMesh,
361  const polyMesh& tgtMesh
362 )
363 {
364  initialise(srcMesh, tgtMesh);
365 
366  // Determine (potentially) participating source mesh cells
367  const labelList srcCellIDs(maskCells(srcMesh, tgtMesh));
368 
369  // Initialise list to keep track of whether src cell can be mapped
370  boolList mapFlag(srcMesh.nCells(), false);
371  UIndirectList<bool>(mapFlag, srcCellIDs) = true;
372 
373  // Find initial point in tgt mesh
374  label srcSeedI = -1;
375  label tgtSeedI = -1;
376  label startSeedI = 0;
377  bool startWalk =
378  findInitialSeeds
379  (
380  srcMesh,
381  tgtMesh,
382  srcCellIDs,
383  mapFlag,
384  startSeedI,
385  srcSeedI,
386  tgtSeedI
387  );
388 
389  if (startWalk)
390  {
391  return
392  calculateAddressing
393  (
394  srcMesh,
395  tgtMesh,
396  srcLocalTgtCells_,
397  srcWeights_,
398  tgtLocalSrcCells_,
399  tgtWeights_,
400  srcSeedI,
401  tgtSeedI,
402  srcCellIDs,
403  mapFlag,
404  startSeedI
405  );
406  }
407  else
408  {
409  return 0;
410  }
411 }
412 
413 
415 (
416  const polyMesh& srcMesh,
417  labelListList& srcToTgtAddr,
418  scalarListList& srcToTgtWght
419 ) const
420 {
421  tetOverlapVolume overlapEngine;
422 
423  forAll(srcToTgtWght, srcCelli)
424  {
425  const scalar v = overlapEngine.cellVolumeMinDecomp(srcMesh, srcCelli);
426 
427  forAll(srcToTgtWght[srcCelli], i)
428  {
429  srcToTgtWght[srcCelli][i] /= v;
430  }
431  }
432 }
433 
434 
435 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
436 
438 :
439  cellsToCells()
440 {}
441 
442 
443 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
444 
446 {}
447 
448 
449 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:125
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
Intersection-based cells-to-cells interpolation class. Volume conservative.
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
Calculates the overlap volume of two cells using tetrahedral decomposition.
scalar cellVolumeMinDecomp(const primitiveMesh &mesh, const label celli) const
Calculates the cell volume.
A class for handling words, derived from string.
Definition: word.H:62
defineTypeNameAndDebug(intersection, 0)
addToRunTimeSelectionTable(cellsToCells, intersection, word)
const dimensionedScalar e
Elementary charge.
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
List< cell > cellList
list of cells
Definition: cellList.H:42
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
List< face > faceList
Definition: faceListFwd.H:41