cellVolumeWeightMethod.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) 2013-2016 OpenFOAM Foundation
9  Copyright (C) 2015 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 "cellVolumeWeightMethod.H"
30 #include "indexedOctree.H"
31 #include "treeDataCell.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(cellVolumeWeightMethod, 0);
40  (
41  meshToMeshMethod,
42  cellVolumeWeightMethod,
43  components
44  );
45 }
46 
47 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
48 
50 (
51  const labelList& srcCellIDs,
52  const boolList& mapFlag,
53  const label startSeedI,
54  label& srcSeedI,
55  label& tgtSeedI
56 ) const
57 {
58  const cellList& srcCells = src_.cells();
59  const faceList& srcFaces = src_.faces();
60  const pointField& srcPts = src_.points();
61 
62  for (label i = startSeedI; i < srcCellIDs.size(); ++i)
63  {
64  const label srcI = srcCellIDs[i];
65 
66  if (mapFlag[srcI])
67  {
68  const pointField pts(srcCells[srcI].points(srcFaces, srcPts));
69 
70  for (const point& pt : pts)
71  {
72  const label tgtI = tgt_.cellTree().findInside(pt);
73 
74  if (tgtI != -1 && intersect(srcI, tgtI))
75  {
76  srcSeedI = srcI;
77  tgtSeedI = tgtI;
78 
79  return true;
80  }
81  }
82  }
83  }
84 
85  if (debug)
86  {
87  Pout<< "could not find starting seed" << endl;
88  }
89 
90  return false;
91 }
92 
93 
95 (
96  labelListList& srcToTgtCellAddr,
97  scalarListList& srcToTgtCellWght,
98  labelListList& tgtToSrcCellAddr,
99  scalarListList& tgtToSrcCellWght,
100  const label srcSeedI,
101  const label tgtSeedI,
102  const labelList& srcCellIDs,
103  boolList& mapFlag,
104  label& startSeedI
105 )
106 {
107  label srcCelli = srcSeedI;
108  label tgtCelli = tgtSeedI;
109 
110  List<DynamicList<label>> srcToTgtAddr(src_.nCells());
111  List<DynamicList<scalar>> srcToTgtWght(src_.nCells());
112 
113  List<DynamicList<label>> tgtToSrcAddr(tgt_.nCells());
114  List<DynamicList<scalar>> tgtToSrcWght(tgt_.nCells());
115 
116  // list of tgt cell neighbour cells
117  DynamicList<label> nbrTgtCells(10);
118 
119  // list of tgt cells currently visited for srcCelli to avoid multiple hits
120  DynamicList<label> visitedTgtCells(10);
121 
122  // list to keep track of tgt cells used to seed src cells
123  labelList seedCells(src_.nCells(), -1);
124  seedCells[srcCelli] = tgtCelli;
125 
126  const scalarField& srcVol = src_.cellVolumes();
127 
128  do
129  {
130  nbrTgtCells.clear();
131  visitedTgtCells.clear();
132 
133  // append initial target cell and neighbours
134  nbrTgtCells.append(tgtCelli);
135  appendNbrCells(tgtCelli, tgt_, visitedTgtCells, nbrTgtCells);
136 
137  do
138  {
139  tgtCelli = nbrTgtCells.remove();
140  visitedTgtCells.append(tgtCelli);
141 
142  scalar vol = interVol(srcCelli, tgtCelli);
143 
144  // accumulate addressing and weights for valid intersection
145  if (vol/srcVol[srcCelli] > tolerance_)
146  {
147  // store src/tgt cell pair
148  srcToTgtAddr[srcCelli].append(tgtCelli);
149  srcToTgtWght[srcCelli].append(vol);
150 
151  tgtToSrcAddr[tgtCelli].append(srcCelli);
152  tgtToSrcWght[tgtCelli].append(vol);
153 
154  appendNbrCells(tgtCelli, tgt_, visitedTgtCells, nbrTgtCells);
155 
156  // accumulate intersection volume
157  V_ += vol;
158  }
159  }
160  while (!nbrTgtCells.empty());
161 
162  mapFlag[srcCelli] = false;
163 
164  // find new source seed cell
165  setNextCells
166  (
167  startSeedI,
168  srcCelli,
169  tgtCelli,
170  srcCellIDs,
171  mapFlag,
172  visitedTgtCells,
173  seedCells
174  );
175  }
176  while (srcCelli != -1);
177 
178  // transfer addressing into persistent storage
179  forAll(srcToTgtCellAddr, i)
180  {
181  srcToTgtCellAddr[i].transfer(srcToTgtAddr[i]);
182  srcToTgtCellWght[i].transfer(srcToTgtWght[i]);
183  }
184 
185  forAll(tgtToSrcCellAddr, i)
186  {
187  tgtToSrcCellAddr[i].transfer(tgtToSrcAddr[i]);
188  tgtToSrcCellWght[i].transfer(tgtToSrcWght[i]);
189  }
190 
191 
192  if (debug%2)
193  {
194  // At this point the overlaps are still in volume so we could
195  // get out the relative error
196  forAll(srcToTgtCellAddr, cellI)
197  {
198  scalar srcVol = src_.cellVolumes()[cellI];
199  scalar tgtVol = sum(srcToTgtCellWght[cellI]);
200 
201  if (mag(srcVol) > ROOTVSMALL && mag((tgtVol-srcVol)/srcVol) > 1e-6)
202  {
204  << "At cell " << cellI << " cc:"
205  << src_.cellCentres()[cellI]
206  << " vol:" << srcVol
207  << " total overlap volume:" << tgtVol
208  << endl;
209  }
210  }
211 
212  forAll(tgtToSrcCellAddr, cellI)
213  {
214  scalar tgtVol = tgt_.cellVolumes()[cellI];
215  scalar srcVol = sum(tgtToSrcCellWght[cellI]);
216 
217  if (mag(tgtVol) > ROOTVSMALL && mag((srcVol-tgtVol)/tgtVol) > 1e-6)
218  {
220  << "At cell " << cellI << " cc:"
221  << tgt_.cellCentres()[cellI]
222  << " vol:" << tgtVol
223  << " total overlap volume:" << srcVol
224  << endl;
225  }
226  }
227  }
228 }
229 
230 
232 (
233  label& startSeedI,
234  label& srcCelli,
235  label& tgtCelli,
236  const labelList& srcCellIDs,
237  const boolList& mapFlag,
238  const DynamicList<label>& visitedCells,
239  labelList& seedCells
240 ) const
241 {
242  const labelList& srcNbrCells = src_.cellCells()[srcCelli];
243 
244  // set possible seeds for later use by querying all src cell neighbours
245  // with all visited target cells
246  bool valuesSet = false;
247  forAll(srcNbrCells, i)
248  {
249  label cellS = srcNbrCells[i];
250 
251  if (mapFlag[cellS] && seedCells[cellS] == -1)
252  {
253  forAll(visitedCells, j)
254  {
255  label cellT = visitedCells[j];
256 
257  if (intersect(cellS, cellT))
258  {
259  seedCells[cellS] = cellT;
260 
261  if (!valuesSet)
262  {
263  srcCelli = cellS;
264  tgtCelli = cellT;
265  valuesSet = true;
266  }
267  }
268  }
269  }
270  }
271 
272  // set next src and tgt cells if not set above
273  if (valuesSet)
274  {
275  return;
276  }
277  else
278  {
279  // try to use existing seed
280  bool foundNextSeed = false;
281  for (label i = startSeedI; i < srcCellIDs.size(); i++)
282  {
283  label cellS = srcCellIDs[i];
284 
285  if (mapFlag[cellS])
286  {
287  if (!foundNextSeed)
288  {
289  startSeedI = i;
290  foundNextSeed = true;
291  }
292 
293  if (seedCells[cellS] != -1)
294  {
295  srcCelli = cellS;
296  tgtCelli = seedCells[cellS];
297 
298  return;
299  }
300  }
301  }
302 
303  // perform new search to find match
304  if (debug)
305  {
306  Pout<< "Advancing front stalled: searching for new "
307  << "target cell" << endl;
308  }
309 
310  bool restart =
311  findInitialSeeds
312  (
313  srcCellIDs,
314  mapFlag,
315  startSeedI,
316  srcCelli,
317  tgtCelli
318  );
319 
320  if (restart)
321  {
322  // successfully found new starting seed-pair
323  return;
324  }
325  }
326 
327  // if we have got to here, there are no more src/tgt cell intersections
328  srcCelli = -1;
329  tgtCelli = -1;
330 }
331 
332 
333 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
334 
335 Foam::cellVolumeWeightMethod::cellVolumeWeightMethod
336 (
337  const polyMesh& src,
338  const polyMesh& tgt
339 )
340 :
341  meshToMeshMethod(src, tgt)
342 {}
343 
344 
345 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
346 
348 {}
349 
350 
351 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
352 
354 (
355  labelListList& srcToTgtAddr,
356  scalarListList& srcToTgtWght,
357  pointListList& srcToTgtVec,
358  labelListList& tgtToSrcAddr,
359  scalarListList& tgtToSrcWght,
360  pointListList& tgtToSrcVec
361 )
362 {
363  bool ok = initialise
364  (
365  srcToTgtAddr,
366  srcToTgtWght,
367  tgtToSrcAddr,
368  tgtToSrcWght
369  );
370 
371  if (!ok)
372  {
373  return;
374  }
375 
376  // (potentially) participating source mesh cells
377  const labelList srcCellIDs(maskCells());
378 
379  // list to keep track of whether src cell can be mapped
380  boolList mapFlag(src_.nCells(), false);
381  boolUIndList(mapFlag, srcCellIDs) = true;
382 
383  // find initial point in tgt mesh
384  label srcSeedI = -1;
385  label tgtSeedI = -1;
386  label startSeedI = 0;
387 
388  bool startWalk =
389  findInitialSeeds
390  (
391  srcCellIDs,
392  mapFlag,
393  startSeedI,
394  srcSeedI,
395  tgtSeedI
396  );
397 
398  if (startWalk)
399  {
400  calculateAddressing
401  (
402  srcToTgtAddr,
403  srcToTgtWght,
404  tgtToSrcAddr,
405  tgtToSrcWght,
406  srcSeedI,
407  tgtSeedI,
408  srcCellIDs,
409  mapFlag,
410  startSeedI
411  );
412  }
413  else
414  {
415  // if meshes are collocated, after inflating the source mesh bounding
416  // box tgt mesh cells may be transferred, but may still not overlap
417  // with the source mesh
418  return;
419  }
420 }
421 
422 
423 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
Foam::cellVolumeWeightMethod::~cellVolumeWeightMethod
virtual ~cellVolumeWeightMethod()
Destructor.
Definition: cellVolumeWeightMethod.C:347
Foam::DynamicList< label >
indexedOctree.H
Foam::List::append
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
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.
Foam::DynamicList::clear
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:391
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::meshToMeshMethod
Base class for mesh-to-mesh calculation methods.
Definition: meshToMeshMethod.H:51
Foam::boolUIndList
UIndirectList< bool > boolUIndList
UIndirectList of bools.
Definition: UIndirectList.H:54
Foam::Field< vector >
Foam::cellVolumeWeightMethod::findInitialSeeds
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.
Definition: cellVolumeWeightMethod.C:50
Foam::DynamicList::append
DynamicList< T, SizeMin > & append(const T &val)
Append an element to the end of this list.
Definition: DynamicListI.H:511
Foam::List::transfer
void transfer(List< T > &list)
Definition: List.C:456
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
treeDataCell.H
Foam::cellVolumeWeightMethod::setNextCells
void setNextCells(label &startSeedI, label &srcCelli, label &tgtCelli, const labelList &srcCellIDs, const boolList &mapFlag, const DynamicList< label > &visitedCells, labelList &seedCells) const
Set the next cells in the advancing front algorithm.
Definition: cellVolumeWeightMethod.C:232
Foam::cellVolumeWeightMethod::calculate
virtual void calculate(labelListList &srcToTgtAddr, scalarListList &srcToTgtWght, pointListList &srcToTgtVec, labelListList &tgtToSrcAddr, scalarListList &tgtToSrcWght, pointListList &tgtToSrcVec)
Calculate addressing and weights and optionally offset vectors.
Definition: cellVolumeWeightMethod.C:354
Foam::Vector< scalar >
Foam::List< label >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
cellVolumeWeightMethod.H
Foam::cellVolumeWeightMethod::calculateAddressing
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.
Definition: cellVolumeWeightMethod.C:95
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328