pairPatchAgglomeration.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) 2011-2020 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 "pairPatchAgglomeration.H"
27 #include "meshTools.H"
28 #include "unitConversion.H"
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 void Foam::pairPatchAgglomeration::compactLevels(const label nCreatedLevels)
33 {
34  nFaces_.setSize(nCreatedLevels);
35  restrictAddressing_.setSize(nCreatedLevels);
36  patchLevels_.setSize(nCreatedLevels);
37 }
38 
39 
40 bool Foam::pairPatchAgglomeration::continueAgglomerating
41 (
42  const label nCoarseFaces
43 )
44 {
45  // Check the need for further agglomeration on all processors
46  label localnCoarseFaces = nCoarseFaces;
47  bool contAgg = localnCoarseFaces >= nFacesInCoarsestLevel_;
48  return contAgg;
49 }
50 
51 
52 void Foam::pairPatchAgglomeration::setBasedEdgeWeights()
53 {
54  const bPatch& coarsePatch = patchLevels_[0];
55  forAll(coarsePatch.edges(), i)
56  {
57  if (coarsePatch.isInternalEdge(i))
58  {
59  scalar edgeLength =
60  coarsePatch.edges()[i].mag(coarsePatch.localPoints());
61 
62  const labelList& eFaces = coarsePatch.edgeFaces()[i];
63 
64  if (eFaces.size() == 2)
65  {
66  scalar cosI =
67  coarsePatch.faceNormals()[eFaces[0]]
68  & coarsePatch.faceNormals()[eFaces[1]];
69 
70  const edge edgeCommon = edge(eFaces[0], eFaces[1]);
71 
72  if (facePairWeight_.found(edgeCommon))
73  {
74  facePairWeight_[edgeCommon] += edgeLength;
75  }
76  else
77  {
78  facePairWeight_.insert(edgeCommon, edgeLength);
79  }
80 
81  if (cosI < Foam::cos(degToRad(featureAngle_)))
82  {
83  facePairWeight_[edgeCommon] = -1.0;
84  }
85  }
86  else
87  {
88  forAll(eFaces, j)
89  {
90  for (label k = j+1; k<eFaces.size(); k++)
91  {
93  (
94  edge(eFaces[j], eFaces[k]),
95  -1.0
96  );
97  }
98  }
99  }
100  }
101  }
102 }
103 
104 
105 void Foam::pairPatchAgglomeration::setEdgeWeights
106 (
107  const label fineLevelIndex
108 )
109 {
110 
111  const bPatch& coarsePatch = patchLevels_[fineLevelIndex];
112 
113  const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
114  const label nCoarseI = max(fineToCoarse) + 1;
115  labelListList coarseToFine(invertOneToMany(nCoarseI, fineToCoarse));
116 
117  HashSet<edge, Hash<edge>> fineFeaturedFaces(coarsePatch.nEdges()/10);
118 
119  // Map fine faces with featured edge into coarse faces
120  forAllConstIter(EdgeMap<scalar>, facePairWeight_, iter)
121  {
122  if (iter() == -1.0)
123  {
124  const edge e = iter.key();
125  const edge edgeFeatured
126  (
127  fineToCoarse[e[0]],
128  fineToCoarse[e[1]]
129  );
130  fineFeaturedFaces.insert(edgeFeatured);
131  }
132  }
133 
134  // Clean old weitghs
136  facePairWeight_.resize(coarsePatch.nEdges());
137 
138  forAll(coarsePatch.edges(), i)
139  {
140  if (coarsePatch.isInternalEdge(i))
141  {
142  scalar edgeLength =
143  coarsePatch.edges()[i].mag(coarsePatch.localPoints());
144 
145  const labelList& eFaces = coarsePatch.edgeFaces()[i];
146 
147  if (eFaces.size() == 2)
148  {
149  const edge edgeCommon = edge(eFaces[0], eFaces[1]);
150 
151  if (facePairWeight_.found(edgeCommon))
152  {
153  facePairWeight_[edgeCommon] += edgeLength;
154  }
155  else
156  {
157  facePairWeight_.insert(edgeCommon, edgeLength);
158  }
159 
160  // If the fine 'pair' faces was featured edge so it is
161  // the coarse 'pair'
162  if (fineFeaturedFaces.found(edgeCommon))
163  {
164  facePairWeight_[edgeCommon] = -1.0;
165  }
166  }
167  else
168  {
169  // Set edge as barrier by setting weight to -1
170  forAll(eFaces, j)
171  {
172  for (label k = j+1; k<eFaces.size(); k++)
173  {
175  (
176  edge(eFaces[j], eFaces[k]),
177  -1.0
178  );
179  }
180  }
181  }
182  }
183  }
184 }
185 
186 
187 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
188 
190 (
191  const polyPatch& patch,
192  const dictionary& controlDict,
193  const bool additionalWeights
194 )
195 :
197  (
198  controlDict.lookupOrDefault<label>("mergeLevels", 2)
199  ),
200  maxLevels_(50),
202  (
203  controlDict.lookup<label>("nFacesInCoarsestLevel")
204  ),
206  (
207  controlDict.lookupOrDefault<scalar>("featureAngle", 0)
208  ),
211  restrictTopBottomAddressing_(identity(patch.size())),
213  facePairWeight_(patch.size())
214 {
215  // Set base fine patch
216  patchLevels_.set
217  (
218  0,
219  new bPatch
220  (
221  patch.localFaces(),
222  patch.localPoints()
223  )
224  );
225 
226  // Set number of faces for the base patch
227  nFaces_[0] = patch.size();
228 
229  // Set edge weights for level 0
230  setBasedEdgeWeights();
231 }
232 
233 
234 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
235 
237 {}
238 
239 
240 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
241 
243 (
244  const label i
245 ) const
246 {
247  return patchLevels_[i];
248 }
249 
250 
251 void Foam::pairPatchAgglomeration::mapBaseToTopAgglom
252 (
253  const label fineLevelIndex
254 )
255 {
256  const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
258  {
260  fineToCoarse[restrictTopBottomAddressing_[i]];
261  }
262 }
263 
264 
265 bool Foam::pairPatchAgglomeration::agglomeratePatch
266 (
267  const bPatch& patch,
268  const labelList& fineToCoarse,
269  const label fineLevelIndex
270 )
271 {
272  if (min(fineToCoarse) == -1)
273  {
275  << "min(fineToCoarse) == -1" << exit(FatalError);
276  }
277 
278  if (fineToCoarse.size() == 0)
279  {
280  return true;
281  }
282 
283  if (fineToCoarse.size() != patch.size())
284  {
286  << "restrict map does not correspond to fine level. " << endl
287  << " Sizes: restrictMap: " << fineToCoarse.size()
288  << " nEqns: " << patch.size()
289  << abort(FatalError);
290  }
291 
292  const label nCoarseI = max(fineToCoarse)+1;
293  List<face> patchFaces(nCoarseI);
294 
295 
296  // Patch faces per agglomeration
297  labelListList coarseToFine(invertOneToMany(nCoarseI, fineToCoarse));
298 
299  for (label coarseI = 0; coarseI < nCoarseI; coarseI++)
300  {
301  const labelList& fineFaces = coarseToFine[coarseI];
302 
303  // Construct single face
305  (
306  IndirectList<face>(patch, fineFaces),
307  patch.points()
308  );
309 
310  if (upp.edgeLoops().size() != 1)
311  {
312  if (fineFaces.size() == 2)
313  {
314  const edge e(fineFaces[0], fineFaces[1]);
315  facePairWeight_[e] = -1.0;
316  }
317  else if (fineFaces.size() == 3)
318  {
319  const edge e(fineFaces[0], fineFaces[1]);
320  const edge e1(fineFaces[0], fineFaces[2]);
321  const edge e2(fineFaces[2], fineFaces[1]);
322  facePairWeight_[e] = -1.0;
323  facePairWeight_[e1] = -1.0;
324  facePairWeight_[e2] = -1.0;
325  }
326 
327  return false;
328  }
329 
330  patchFaces[coarseI] = face
331  (
332  renumber
333  (
334  upp.meshPoints(),
335  upp.edgeLoops()[0]
336  )
337  );
338  }
339 
340  patchLevels_.set
341  (
342  fineLevelIndex,
343  new bPatch
344  (
345  SubList<face>(patchFaces, nCoarseI, 0),
346  patch.points()
347  )
348  );
349 
350  return true;
351 }
352 
353 
355 {
356  label nPairLevels = 0;
357  label nCreatedLevels = 1; // 0 level is the base patch
358  label nCoarseFaces = 0;
359  label nCoarseFacesOld = 0;
360 
361  while (nCreatedLevels < maxLevels_)
362  {
363  const bPatch& patch = patchLevels_[nCreatedLevels - 1];
364  tmp<labelField> finalAgglomPtr(new labelField(patch.size()));
365  bool agglomOK = false;
366 
367  do
368  {
369  label nCoarseFacesPrev = nCoarseFaces;
370 
371  finalAgglomPtr = agglomerateOneLevel
372  (
373  nCoarseFaces,
374  patch
375  );
376 
377  if (nCoarseFaces > 0 && nCoarseFaces != nCoarseFacesPrev)
378  {
379  if
380  (
381  (
382  agglomOK = agglomeratePatch
383  (
384  patch,
385  finalAgglomPtr,
386  nCreatedLevels
387  )
388  )
389  )
390  {
391  restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
392  mapBaseToTopAgglom(nCreatedLevels);
393  setEdgeWeights(nCreatedLevels);
394 
395  if (nPairLevels % mergeLevels_)
396  {
397  combineLevels(nCreatedLevels);
398  }
399  else
400  {
401  nCreatedLevels++;
402  }
403 
404  nPairLevels++;
405  }
406  }
407  else
408  {
409  agglomOK = true;
410  }
411 
412  reduce(nCoarseFaces, sumOp<label>());
413  reduce(agglomOK, orOp<bool>());
414 
415  } while (!agglomOK);
416 
417  nFaces_[nCreatedLevels] = nCoarseFaces;
418 
419  if
420  (
421  !continueAgglomerating(nCoarseFaces)
422  || (nCoarseFacesOld == nCoarseFaces)
423  )
424  {
425  break;
426  }
427 
428  nCoarseFacesOld = nCoarseFaces;
429  }
430 
431  compactLevels(nCreatedLevels);
432 }
433 
434 
435 Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
436 (
437  label& nCoarseFaces,
438  const bPatch& patch
439 )
440 {
441  const label nFineFaces = patch.size();
442 
443  tmp<labelField> tcoarseCellMap(new labelField(nFineFaces, -1));
444  labelField& coarseCellMap = tcoarseCellMap.ref();
445 
446  const labelListList& faceFaces = patch.faceFaces();
447 
448  nCoarseFaces = 0;
449 
450  forAll(faceFaces, facei)
451  {
452  const labelList& fFaces = faceFaces[facei];
453 
454  if (coarseCellMap[facei] < 0)
455  {
456  label matchFaceNo = -1;
457  label matchFaceNeibNo = -1;
458  scalar maxFaceWeight = -great;
459 
460  // Check faces to find ungrouped neighbour with largest face weight
461  forAll(fFaces, i)
462  {
463  label faceNeig = fFaces[i];
464  const edge edgeCommon = edge(facei, faceNeig);
465  if
466  (
467  facePairWeight_[edgeCommon] > maxFaceWeight
468  && coarseCellMap[faceNeig] < 0
469  && facePairWeight_[edgeCommon] != -1.0
470  )
471  {
472  // Match found. Pick up all the necessary data
473  matchFaceNo = facei;
474  matchFaceNeibNo = faceNeig;
475  maxFaceWeight = facePairWeight_[edgeCommon];
476  }
477  }
478 
479  if (matchFaceNo >= 0)
480  {
481  // Make a new group
482  coarseCellMap[matchFaceNo] = nCoarseFaces;
483  coarseCellMap[matchFaceNeibNo] = nCoarseFaces;
484  nCoarseFaces++;
485  }
486  else
487  {
488  // No match. Find the best neighbouring cluster and
489  // put the cell there
490  label clusterMatchFaceNo = -1;
491  scalar clusterMaxFaceCoeff = -great;
492 
493  forAll(fFaces, i)
494  {
495  label faceNeig = fFaces[i];
496  const edge edgeCommon = edge(facei, faceNeig);
497  if
498  (
499  facePairWeight_[edgeCommon] > clusterMaxFaceCoeff
500  && facePairWeight_[edgeCommon] != -1.0
501  && coarseCellMap[faceNeig] > 0
502  )
503  {
504  clusterMatchFaceNo = faceNeig;
505  clusterMaxFaceCoeff = facePairWeight_[edgeCommon];
506  }
507  }
508 
509  if (clusterMatchFaceNo >= 0)
510  {
511  // Add the cell to the best cluster
512  coarseCellMap[facei] = coarseCellMap[clusterMatchFaceNo];
513  }
514  else
515  {
516  // if not create single-cell "clusters" for each
517  coarseCellMap[facei] = nCoarseFaces;
518  nCoarseFaces++;
519  }
520  }
521  }
522  }
523 
524  // Check that all faces are part of clusters,
525 
526  for (label facei=0; facei<nFineFaces; facei++)
527  {
528  if (coarseCellMap[facei] < 0)
529  {
531  << " face " << facei
532  << " is not part of a cluster"
533  << exit(FatalError);
534  }
535  }
536 
537  return tcoarseCellMap;
538 }
539 
540 
541 void Foam::pairPatchAgglomeration::combineLevels(const label curLevel)
542 {
543  label prevLevel = curLevel - 1;
544 
545  // Set the previous level nCells to the current
546  nFaces_[prevLevel] = nFaces_[curLevel];
547 
548  // Map the restrictAddressing from the coarser level into the previous
549  // finer level
550 
551  const labelList& curResAddr = restrictAddressing_[curLevel];
552  labelList& prevResAddr = restrictAddressing_[prevLevel];
553 
554  forAll(prevResAddr, i)
555  {
556  prevResAddr[i] = curResAddr[prevResAddr[i]];
557  }
558 
559  // Delete the restrictAddressing for the coarser level
560  restrictAddressing_.set(curLevel, nullptr);
561 
562  patchLevels_.set(prevLevel, patchLevels_.set(curLevel, nullptr));
563 }
564 
565 
566 // ************************************************************************* //
ListType renumber(const labelUList &oldToNew, const ListType &)
Renumber the values (not the indices) of a list.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
Field< label > labelField
Specialisation of Field<T> for label.
Definition: labelField.H:49
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
pairPatchAgglomeration(const polyPatch &patch, const dictionary &controlDict, const bool additionalWeights=false)
Construct given mesh and controls.
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
PtrList< labelField > restrictAddressing_
Cell restriction addressing array.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
Unit conversion functions.
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:67
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
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
label k
Boltzmann constant.
const Field< PointType > & localPoints() const
Return pointField of points in patch.
PtrList< bPatch > patchLevels_
Hierarchy of patch addressing.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
A list of faces which address into the list of points.
A List obtained as a section of another List.
Definition: SubList.H:53
const labelListList & faceFaces() const
Return face-face addressing.
dimensionedScalar cos(const dimensionedScalar &ds)
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
void agglomerate()
Agglomerate patch.
void clear()
Clear all entries from table.
Definition: HashTable.C:468
const Field< PointType > & points() const
Return reference to global points.
label nFacesInCoarsestLevel_
Number of faces in coarsest level.
List< label > labelList
A List of labels.
Definition: labelList.H:56
scalar featureAngle_
Feature angle.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
const bPatch & patchLevel(const label leveli) const
Return primitivePatch of given level.
void setSize(const label)
Reset size of List.
Definition: List.C:281
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
labelList restrictTopBottomAddressing_
Maps from finest to coarsest.
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
EdgeMap< scalar > facePairWeight_
Edge weights.
labelList nFaces_
The number of faces in each level.
PrimitivePatch< faceList, const pointField > bPatch
void resize(const label newSize)
Resize the hash table for efficiency.
Definition: HashTable.C:432
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
A class for managing temporary objects.
Definition: PtrList.H:53
label maxLevels_
Max number of levels.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
A List with indirect addressing.
Definition: IndirectList.H:101
label mergeLevels_
Number of levels to merge, 1 = don&#39;t merge, 2 = merge pairs etc.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844