faceAreaWeightAMI.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-2018 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 "faceAreaWeightAMI.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  defineTypeNameAndDebug(faceAreaWeightAMI, 0);
34  addToRunTimeSelectionTable(AMIMethod, faceAreaWeightAMI, components);
35 }
36 
37 
38 // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
39 
41 (
42  List<DynamicList<label>>& srcAddr,
43  List<DynamicList<scalar>>& srcWght,
44  List<DynamicList<label>>& tgtAddr,
45  List<DynamicList<scalar>>& tgtWght,
46  label srcFacei,
47  label tgtFacei
48 )
49 {
50  // construct weights and addressing
51  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
52 
53  label nFacesRemaining = srcAddr.size();
54 
55  // list of tgt face neighbour faces
56  DynamicList<label> nbrFaces(10);
57 
58  // list of faces currently visited for srcFacei to avoid multiple hits
59  DynamicList<label> visitedFaces(10);
60 
61  // list to keep track of tgt faces used to seed src faces
62  labelList seedFaces(nFacesRemaining, -1);
63  seedFaces[srcFacei] = tgtFacei;
64 
65  // list to keep track of whether src face can be mapped
66  boolList mapFlag(nFacesRemaining, true);
67 
68  // reset starting seed
69  label startSeedI = 0;
70 
71  DynamicList<label> nonOverlapFaces;
72  do
73  {
74  // Do advancing front starting from srcFacei,tgtFacei
75  bool faceProcessed = processSourceFace
76  (
77  srcFacei,
78  tgtFacei,
79 
80  nbrFaces,
81  visitedFaces,
82 
83  srcAddr,
84  srcWght,
85  tgtAddr,
86  tgtWght
87  );
88 
89  mapFlag[srcFacei] = false;
90 
91  nFacesRemaining--;
92 
93  if (!faceProcessed)
94  {
95  nonOverlapFaces.append(srcFacei);
96  }
97 
98  // choose new src face from current src face neighbour
99  if (nFacesRemaining > 0)
100  {
101  setNextFaces
102  (
103  startSeedI,
104  srcFacei,
105  tgtFacei,
106  mapFlag,
107  seedFaces,
108  visitedFaces
109  );
110  }
111  } while (nFacesRemaining > 0);
112 
113  this->srcNonOverlap_.transfer(nonOverlapFaces);
114 }
115 
116 
118 (
119  const label srcFacei,
120  const label tgtStartFacei,
121 
122  // list of tgt face neighbour faces
123  DynamicList<label>& nbrFaces,
124  // list of faces currently visited for srcFacei to avoid multiple hits
125  DynamicList<label>& visitedFaces,
126 
127  // temporary storage for addressing and weights
128  List<DynamicList<label>>& srcAddr,
129  List<DynamicList<scalar>>& srcWght,
130  List<DynamicList<label>>& tgtAddr,
131  List<DynamicList<scalar>>& tgtWght
132 )
133 {
134  if (tgtStartFacei == -1)
135  {
136  return false;
137  }
138 
139  nbrFaces.clear();
140  visitedFaces.clear();
141 
142  // append initial target face and neighbours
143  nbrFaces.append(tgtStartFacei);
144  this->appendNbrFaces
145  (
146  tgtStartFacei,
147  this->tgtPatch_,
148  visitedFaces,
149  nbrFaces
150  );
151 
152  const scalar srcArea = this->srcMagSf_[srcFacei];
153 
154  bool faceProcessed = false;
155 
156  do
157  {
158  // process new target face
159  const label tgtFacei = nbrFaces.remove();
160  visitedFaces.append(tgtFacei);
161  const scalar tgtArea = this->tgtMagSf_[tgtFacei];
162 
163  // calculate the intersection area
164  const scalar area = interArea(srcFacei, tgtFacei);
165 
166  // store when intersection fractional area > min weight
167  if (area/srcArea > minWeight())
168  {
169  srcAddr[srcFacei].append(tgtFacei);
170  srcWght[srcFacei].append(area/srcArea);
171 
172  tgtAddr[tgtFacei].append(srcFacei);
173  tgtWght[tgtFacei].append(area/tgtArea);
174 
175  this->appendNbrFaces
176  (
177  tgtFacei,
178  this->tgtPatch_,
179  visitedFaces,
180  nbrFaces
181  );
182 
183  faceProcessed = true;
184  }
185 
186  } while (nbrFaces.size() > 0);
187 
188  return faceProcessed;
189 }
190 
191 
193 (
194  label& startSeedI,
195  label& srcFacei,
196  label& tgtFacei,
197  const boolList& mapFlag,
198  labelList& seedFaces,
199  const DynamicList<label>& visitedFaces,
200  bool errorOnNotFound
201 ) const
202 {
203  const labelList& srcNbrFaces = this->srcPatch_.faceFaces()[srcFacei];
204 
205  // initialise tgtFacei
206  tgtFacei = -1;
207 
208  // set possible seeds for later use
209  bool valuesSet = false;
210  forAll(srcNbrFaces, i)
211  {
212  label faceS = srcNbrFaces[i];
213 
214  if (mapFlag[faceS] && seedFaces[faceS] == -1)
215  {
216  forAll(visitedFaces, j)
217  {
218  label faceT = visitedFaces[j];
219  scalar area = interArea(faceS, faceT);
220  scalar areaTotal = this->srcMagSf_[srcFacei];
221 
222  // Check that faces have enough overlap for robust walking
223  if (area/areaTotal > minWeight())
224  {
225  // TODO - throwing area away - re-use in next iteration?
226 
227  seedFaces[faceS] = faceT;
228 
229  if (!valuesSet)
230  {
231  srcFacei = faceS;
232  tgtFacei = faceT;
233  valuesSet = true;
234  }
235  }
236  }
237  }
238  }
239 
240  // set next src and tgt faces if not set above
241  if (valuesSet)
242  {
243  return;
244  }
245  else
246  {
247  // try to use existing seed
248  bool foundNextSeed = false;
249  for (label facei = startSeedI; facei < mapFlag.size(); facei++)
250  {
251  if (mapFlag[facei])
252  {
253  if (!foundNextSeed)
254  {
255  startSeedI = facei;
256  foundNextSeed = true;
257  }
258 
259  if (seedFaces[facei] != -1)
260  {
261  srcFacei = facei;
262  tgtFacei = seedFaces[facei];
263 
264  return;
265  }
266  }
267  }
268 
269  // perform new search to find match
270  if (debug)
271  {
272  Pout<< "Advancing front stalled: searching for new "
273  << "target face" << endl;
274  }
275 
276  foundNextSeed = false;
277  for (label facei = startSeedI; facei < mapFlag.size(); facei++)
278  {
279  if (mapFlag[facei])
280  {
281  if (!foundNextSeed)
282  {
283  startSeedI = facei + 1;
284  foundNextSeed = true;
285  }
286 
287  srcFacei = facei;
288  tgtFacei = this->findTargetFace(srcFacei);
289 
290  if (tgtFacei >= 0)
291  {
292  return;
293  }
294  }
295  }
296 
297  if (errorOnNotFound)
298  {
300  << "Unable to set source and target faces" << abort(FatalError);
301  }
302  }
303 }
304 
305 
307 (
308  const label srcFacei,
309  const label tgtFacei
310 ) const
311 {
312  scalar area = 0;
313 
314  const pointField& srcPoints = this->srcPatch_.points();
315  const pointField& tgtPoints = this->tgtPatch_.points();
316 
317  // references to candidate faces
318  const face& src = this->srcPatch_[srcFacei];
319  const face& tgt = this->tgtPatch_[tgtFacei];
320 
321  // quick reject if either face has zero area
322  // Note: do not use stored face areas for target patch
323  const scalar tgtMag = tgt.mag(tgtPoints);
324  if ((this->srcMagSf_[srcFacei] < rootVSmall) || (tgtMag < rootVSmall))
325  {
326  return area;
327  }
328 
329  // create intersection object
330  faceAreaIntersect inter(srcPoints, tgtPoints, this->reverseTarget_);
331 
332  // crude resultant norm
333  vector n(-this->srcPatch_.faceNormals()[srcFacei]);
334  if (this->reverseTarget_)
335  {
336  n -= this->tgtPatch_.faceNormals()[tgtFacei];
337  }
338  else
339  {
340  n += this->tgtPatch_.faceNormals()[tgtFacei];
341  }
342  scalar magN = mag(n);
343 
344  if (magN > rootVSmall)
345  {
346  area = inter.calc(src, tgt, n/magN, this->triMode_);
347  }
348  else
349  {
351  << "Invalid normal for source face " << srcFacei
352  << " points " << UIndirectList<point>(srcPoints, src)
353  << " target face " << tgtFacei
354  << " points " << UIndirectList<point>(tgtPoints, tgt)
355  << endl;
356  }
357 
358 
359  if ((debug > 1) && (area > 0))
360  {
361  this->writeIntersectionOBJ(area, src, tgt, srcPoints, tgtPoints);
362  }
363 
364  return area;
365 }
366 
367 
369 (
370  List<DynamicList<label>>& srcAddr,
371  List<DynamicList<scalar>>& srcWght,
372  List<DynamicList<label>>& tgtAddr,
373  List<DynamicList<scalar>>& tgtWght
374 )
375 {
376  // Collect all src faces with a low weight
377  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
378 
379  labelHashSet lowWeightFaces(100);
380  forAll(srcWght, srcFacei)
381  {
382  const scalar s = sum(srcWght[srcFacei]);
383 
384  if (s < 0.5)
385  {
386  lowWeightFaces.insert(srcFacei);
387  }
388  }
389 
390  if (debug)
391  {
392  Pout<< "faceAreaWeightAMI: restarting search on "
393  << lowWeightFaces.size() << " faces since sum of weights < 0.5"
394  << endl;
395  }
396 
397  if (lowWeightFaces.size() > 0)
398  {
399  // Erase all the lowWeight source faces from the target
400  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
401 
402  DynamicList<label> okSrcFaces(10);
403  DynamicList<scalar> okSrcWeights(10);
404  forAll(tgtAddr, tgtFacei)
405  {
406  okSrcFaces.clear();
407  okSrcWeights.clear();
408  DynamicList<label>& srcFaces = tgtAddr[tgtFacei];
409  DynamicList<scalar>& srcWeights = tgtWght[tgtFacei];
410  forAll(srcFaces, i)
411  {
412  if (!lowWeightFaces.found(srcFaces[i]))
413  {
414  okSrcFaces.append(srcFaces[i]);
415  okSrcWeights.append(srcWeights[i]);
416  }
417  }
418  if (okSrcFaces.size() < srcFaces.size())
419  {
420  srcFaces.transfer(okSrcFaces);
421  srcWeights.transfer(okSrcWeights);
422  }
423  }
424 
425 
426 
427  // Restart search from best hit
428  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
429 
430  // list of tgt face neighbour faces
431  DynamicList<label> nbrFaces(10);
432 
433  // list of faces currently visited for srcFacei to avoid multiple hits
434  DynamicList<label> visitedFaces(10);
435 
436  forAllConstIter(labelHashSet, lowWeightFaces, iter)
437  {
438  label srcFacei = iter.key();
439  label tgtFacei = this->findTargetFace(srcFacei);
440  if (tgtFacei != -1)
441  {
442  // bool faceProcessed =
443  processSourceFace
444  (
445  srcFacei,
446  tgtFacei,
447 
448  nbrFaces,
449  visitedFaces,
450 
451  srcAddr,
452  srcWght,
453  tgtAddr,
454  tgtWght
455  );
456  // ? Check faceProcessed to see if restarting has worked.
457  }
458  }
459  }
460 }
461 
462 
463 Foam::scalar
465 {
467 }
468 
469 
470 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
471 
473 (
474  const primitivePatch& srcPatch,
475  const primitivePatch& tgtPatch,
476  const scalarField& srcMagSf,
477  const scalarField& tgtMagSf,
479  const bool reverseTarget,
480  const bool requireMatch,
481  const bool restartUncoveredSourceFace
482 )
483 :
484  AMIMethod
485  (
486  srcPatch,
487  tgtPatch,
488  srcMagSf,
489  tgtMagSf,
490  triMode,
491  reverseTarget,
492  requireMatch
493  ),
494  restartUncoveredSourceFace_(restartUncoveredSourceFace)
495 {}
496 
497 
498 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
499 
501 {}
502 
503 
504 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
505 
507 (
508  labelListList& srcAddress,
509  scalarListList& srcWeights,
510  labelListList& tgtAddress,
511  scalarListList& tgtWeights,
512  label srcFacei,
513  label tgtFacei
514 )
515 {
516  bool ok =
517  this->initialise
518  (
519  srcAddress,
520  srcWeights,
521  tgtAddress,
522  tgtWeights,
523  srcFacei,
524  tgtFacei
525  );
526 
527  if (!ok)
528  {
529  return;
530  }
531 
532  // temporary storage for addressing and weights
533  List<DynamicList<label>> srcAddr(this->srcPatch_.size());
534  List<DynamicList<scalar>> srcWght(srcAddr.size());
535  List<DynamicList<label>> tgtAddr(this->tgtPatch_.size());
536  List<DynamicList<scalar>> tgtWght(tgtAddr.size());
537 
538  calcAddressing
539  (
540  srcAddr,
541  srcWght,
542  tgtAddr,
543  tgtWght,
544  srcFacei,
545  tgtFacei
546  );
547 
548  if (debug && !this->srcNonOverlap_.empty())
549  {
550  Pout<< " AMI: " << this->srcNonOverlap_.size()
551  << " non-overlap faces identified"
552  << endl;
553  }
554 
555 
556  // Check for badly covered faces
557  if (restartUncoveredSourceFace_)
558  {
559  restartUncoveredSourceFace
560  (
561  srcAddr,
562  srcWght,
563  tgtAddr,
564  tgtWght
565  );
566  }
567 
568 
569  // transfer data to persistent storage
570  forAll(srcAddr, i)
571  {
572  srcAddress[i].transfer(srcAddr[i]);
573  srcWeights[i].transfer(srcWght[i]);
574  }
575  forAll(tgtAddr, i)
576  {
577  tgtAddress[i].transfer(tgtAddr[i]);
578  tgtWeights[i].transfer(tgtWght[i]);
579  }
580 }
581 
582 
583 // ************************************************************************* //
virtual void calcAddressing(List< DynamicList< label >> &srcAddress, List< DynamicList< scalar >> &srcWeights, List< DynamicList< label >> &tgtAddress, List< DynamicList< scalar >> &tgtWeights, label srcFacei, label tgtFacei)
Calculate addressing and weights using temporary storage.
virtual scalar interArea(const label srcFacei, const label tgtFacei) const
Area of intersection between source and target faces.
virtual bool processSourceFace(const label srcFacei, const label tgtStartFacei, DynamicList< label > &nbrFaces, DynamicList< label > &visitedFaces, List< DynamicList< label >> &srcAddr, List< DynamicList< scalar >> &srcWght, List< DynamicList< label >> &tgtAddr, List< DynamicList< scalar >> &tgtWght)
Determine overlap contributions for source face srcFacei.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
static scalar & tolerance()
Fraction of local length scale to use as intersection tolerance.
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
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
virtual ~faceAreaWeightAMI()
Destructor.
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
virtual void calculate(labelListList &srcAddress, scalarListList &srcWeights, labelListList &tgtAddress, scalarListList &tgtWeights, label srcFacei=-1, label tgtFacei=-1)
Update addressing and weights.
virtual scalar minWeight() const
The minimum weight below which connections are discarded.
A list of faces which address into the list of points.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: faceI.H:81
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
scalar calc(const face &faceA, const face &faceB, const vector &n, const triangulationMode &triMode)
Return area of intersection of faceA with faceB.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
virtual void setNextFaces(label &startSeedI, label &srcFacei, label &tgtFacei, const boolList &mapFlag, labelList &seedFaces, const DynamicList< label > &visitedFaces, bool errorOnNotFound=true) const
Set the source and target seed faces.
Base class for Arbitrary Mesh Interface (AMI) methods.
Definition: AMIMethod.H:55
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
Face intersection class.
defineTypeNameAndDebug(combustionModel, 0)
T remove()
Remove and return the top element.
Definition: DynamicListI.H:351
virtual void restartUncoveredSourceFace(List< DynamicList< label >> &srcAddr, List< DynamicList< scalar >> &srcWght, List< DynamicList< label >> &tgtAddr, List< DynamicList< scalar >> &tgtWght)
Attempt to re-evaluate source faces that have not been included.
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
A List with indirect addressing.
Definition: fvMatrix.H:106
faceAreaWeightAMI(const primitivePatch &srcPatch, const primitivePatch &tgtPatch, const scalarField &srcMagSf, const scalarField &tgtMagSf, const faceAreaIntersect::triangulationMode &triMode, const bool reverseTarget=false, const bool requireMatch=true, const bool restartUncoveredSourceFace=true)
Construct from components.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:236
void transfer(List< T > &)
Transfer contents of the argument List into this.
Definition: DynamicListI.H:285
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
Namespace for OpenFOAM.