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 | Copyright (C) 2013-2016 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 "cellVolumeWeightMethod.H"
27 #include "indexedOctree.H"
28 #include "treeDataCell.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  defineTypeNameAndDebug(cellVolumeWeightMethod, 0);
37  (
38  meshToMeshMethod,
39  cellVolumeWeightMethod,
40  components
41  );
42 }
43 
44 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
45 
47 (
48  const labelList& srcCellIDs,
49  const boolList& mapFlag,
50  const label startSeedI,
51  label& srcSeedI,
52  label& tgtSeedI
53 ) const
54 {
55  const cellList& srcCells = src_.cells();
56  const faceList& srcFaces = src_.faces();
57  const pointField& srcPts = src_.points();
58 
59  for (label i = startSeedI; i < srcCellIDs.size(); i++)
60  {
61  label srcI = srcCellIDs[i];
62 
63  if (mapFlag[srcI])
64  {
65  const pointField
66  pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
67 
68  forAll(pts, ptI)
69  {
70  const point& pt = pts[ptI];
71  label tgtI = tgt_.cellTree().findInside(pt);
72 
73  if (tgtI != -1 && intersect(srcI, tgtI))
74  {
75  srcSeedI = srcI;
76  tgtSeedI = tgtI;
77 
78  return true;
79  }
80  }
81  }
82  }
83 
84  if (debug)
85  {
86  Pout<< "could not find starting seed" << endl;
87  }
88 
89  return false;
90 }
91 
92 
94 (
95  labelListList& srcToTgtCellAddr,
96  scalarListList& srcToTgtCellWght,
97  labelListList& tgtToSrcCellAddr,
98  scalarListList& tgtToSrcCellWght,
99  const label srcSeedI,
100  const label tgtSeedI,
101  const labelList& srcCellIDs,
102  boolList& mapFlag,
103  label& startSeedI
104 )
105 {
106  label srcCelli = srcSeedI;
107  label tgtCelli = tgtSeedI;
108 
109  List<DynamicList<label>> srcToTgtAddr(src_.nCells());
110  List<DynamicList<scalar>> srcToTgtWght(src_.nCells());
111 
112  List<DynamicList<label>> tgtToSrcAddr(tgt_.nCells());
113  List<DynamicList<scalar>> tgtToSrcWght(tgt_.nCells());
114 
115  // list of tgt cell neighbour cells
116  DynamicList<label> nbrTgtCells(10);
117 
118  // list of tgt cells currently visited for srcCelli to avoid multiple hits
119  DynamicList<label> visitedTgtCells(10);
120 
121  // list to keep track of tgt cells used to seed src cells
122  labelList seedCells(src_.nCells(), -1);
123  seedCells[srcCelli] = tgtCelli;
124 
125  const scalarField& srcVol = src_.cellVolumes();
126 
127  do
128  {
129  nbrTgtCells.clear();
130  visitedTgtCells.clear();
131 
132  // append initial target cell and neighbours
133  nbrTgtCells.append(tgtCelli);
134  appendNbrCells(tgtCelli, tgt_, visitedTgtCells, nbrTgtCells);
135 
136  do
137  {
138  tgtCelli = nbrTgtCells.remove();
139  visitedTgtCells.append(tgtCelli);
140 
141  scalar vol = interVol(srcCelli, tgtCelli);
142 
143  // accumulate addressing and weights for valid intersection
144  if (vol/srcVol[srcCelli] > tolerance_)
145  {
146  // store src/tgt cell pair
147  srcToTgtAddr[srcCelli].append(tgtCelli);
148  srcToTgtWght[srcCelli].append(vol);
149 
150  tgtToSrcAddr[tgtCelli].append(srcCelli);
151  tgtToSrcWght[tgtCelli].append(vol);
152 
153  appendNbrCells(tgtCelli, tgt_, visitedTgtCells, nbrTgtCells);
154 
155  // accumulate intersection volume
156  V_ += vol;
157  }
158  }
159  while (!nbrTgtCells.empty());
160 
161  mapFlag[srcCelli] = false;
162 
163  // find new source seed cell
164  setNextCells
165  (
166  startSeedI,
167  srcCelli,
168  tgtCelli,
169  srcCellIDs,
170  mapFlag,
171  visitedTgtCells,
172  seedCells
173  );
174  }
175  while (srcCelli != -1);
176 
177  // transfer addressing into persistent storage
178  forAll(srcToTgtCellAddr, i)
179  {
180  srcToTgtCellAddr[i].transfer(srcToTgtAddr[i]);
181  srcToTgtCellWght[i].transfer(srcToTgtWght[i]);
182  }
183 
184  forAll(tgtToSrcCellAddr, i)
185  {
186  tgtToSrcCellAddr[i].transfer(tgtToSrcAddr[i]);
187  tgtToSrcCellWght[i].transfer(tgtToSrcWght[i]);
188  }
189 }
190 
191 
193 (
194  label& startSeedI,
195  label& srcCelli,
196  label& tgtCelli,
197  const labelList& srcCellIDs,
198  const boolList& mapFlag,
199  const DynamicList<label>& visitedCells,
200  labelList& seedCells
201 ) const
202 {
203  const labelList& srcNbrCells = src_.cellCells()[srcCelli];
204 
205  // set possible seeds for later use by querying all src cell neighbours
206  // with all visited target cells
207  bool valuesSet = false;
208  forAll(srcNbrCells, i)
209  {
210  label cellS = srcNbrCells[i];
211 
212  if (mapFlag[cellS] && seedCells[cellS] == -1)
213  {
214  forAll(visitedCells, j)
215  {
216  label cellT = visitedCells[j];
217 
218  if (intersect(cellS, cellT))
219  {
220  seedCells[cellS] = cellT;
221 
222  if (!valuesSet)
223  {
224  srcCelli = cellS;
225  tgtCelli = cellT;
226  valuesSet = true;
227  }
228  }
229  }
230  }
231  }
232 
233  // set next src and tgt cells if not set above
234  if (valuesSet)
235  {
236  return;
237  }
238  else
239  {
240  // try to use existing seed
241  bool foundNextSeed = false;
242  for (label i = startSeedI; i < srcCellIDs.size(); i++)
243  {
244  label cellS = srcCellIDs[i];
245 
246  if (mapFlag[cellS])
247  {
248  if (!foundNextSeed)
249  {
250  startSeedI = i;
251  foundNextSeed = true;
252  }
253 
254  if (seedCells[cellS] != -1)
255  {
256  srcCelli = cellS;
257  tgtCelli = seedCells[cellS];
258 
259  return;
260  }
261  }
262  }
263 
264  // perform new search to find match
265  if (debug)
266  {
267  Pout<< "Advancing front stalled: searching for new "
268  << "target cell" << endl;
269  }
270 
271  bool restart =
272  findInitialSeeds
273  (
274  srcCellIDs,
275  mapFlag,
276  startSeedI,
277  srcCelli,
278  tgtCelli
279  );
280 
281  if (restart)
282  {
283  // successfully found new starting seed-pair
284  return;
285  }
286  }
287 
288  // if we have got to here, there are no more src/tgt cell intersections
289  srcCelli = -1;
290  tgtCelli = -1;
291 }
292 
293 
294 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
295 
296 Foam::cellVolumeWeightMethod::cellVolumeWeightMethod
297 (
298  const polyMesh& src,
299  const polyMesh& tgt
300 )
301 :
302  meshToMeshMethod(src, tgt)
303 {}
304 
305 
306 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
307 
309 {}
310 
311 
312 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
313 
315 (
316  labelListList& srcToTgtAddr,
317  scalarListList& srcToTgtWght,
318  labelListList& tgtToSrcAddr,
319  scalarListList& tgtToSrcWght
320 )
321 {
322  bool ok = initialise
323  (
324  srcToTgtAddr,
325  srcToTgtWght,
326  tgtToSrcAddr,
327  tgtToSrcWght
328  );
329 
330  if (!ok)
331  {
332  return;
333  }
334 
335  // (potentially) participating source mesh cells
336  const labelList srcCellIDs(maskCells());
337 
338  // list to keep track of whether src cell can be mapped
339  boolList mapFlag(src_.nCells(), false);
340  UIndirectList<bool>(mapFlag, srcCellIDs) = true;
341 
342  // find initial point in tgt mesh
343  label srcSeedI = -1;
344  label tgtSeedI = -1;
345  label startSeedI = 0;
346 
347  bool startWalk =
348  findInitialSeeds
349  (
350  srcCellIDs,
351  mapFlag,
352  startSeedI,
353  srcSeedI,
354  tgtSeedI
355  );
356 
357  if (startWalk)
358  {
359  calculateAddressing
360  (
361  srcToTgtAddr,
362  srcToTgtWght,
363  tgtToSrcAddr,
364  tgtToSrcWght,
365  srcSeedI,
366  tgtSeedI,
367  srcCellIDs,
368  mapFlag,
369  startSeedI
370  );
371  }
372  else
373  {
374  // if meshes are collocated, after inflating the source mesh bounding
375  // box tgt mesh cells may be transferred, but may still not overlap
376  // with the source mesh
377  return;
378  }
379 }
380 
381 
382 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:313
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 ~cellVolumeWeightMethod()
Destructor.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
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.
Macros for easy insertion into run-time selection tables.
virtual void calculate(labelListList &srcToTgtAddr, scalarListList &srcToTgtWght, labelListList &tgtToTgtAddr, scalarListList &tgtToTgtWght)
Calculate addressing and weights.
const pointField & points
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:124
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:184
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
Base class for mesh-to-mesh calculation methods.
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.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
A List with indirect addressing.
Definition: fvMatrix.H:106
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:224
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 transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
Namespace for OpenFOAM.