pairPatchAgglomeration.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) 2011-2015 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 
189 Foam::pairPatchAgglomeration::pairPatchAgglomeration
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  readLabel(controlDict.lookup("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  (
276  "pairPatchAgglomeration::agglomeratePatch"
277  "("
278  "const bPatch&, "
279  "const labelList&, "
280  "const label"
281  ")"
282  ) << "min(fineToCoarse) == -1" << exit(FatalError);
283  }
284 
285  if (fineToCoarse.size() == 0)
286  {
287  return true;
288  }
289 
290  if (fineToCoarse.size() != patch.size())
291  {
293  (
294  "pairPatchAgglomeration::agglomeratePatch"
295  "("
296  "const bPatch&, "
297  "const labelList&, "
298  "const label"
299  ")"
300  ) << "restrict map does not correspond to fine level. " << endl
301  << " Sizes: restrictMap: " << fineToCoarse.size()
302  << " nEqns: " << patch.size()
303  << abort(FatalError);
304  }
305 
306  const label nCoarseI = max(fineToCoarse)+1;
307  List<face> patchFaces(nCoarseI);
308 
309 
310  // Patch faces per agglomeration
311  labelListList coarseToFine(invertOneToMany(nCoarseI, fineToCoarse));
312 
313  for (label coarseI = 0; coarseI < nCoarseI; coarseI++)
314  {
315  const labelList& fineFaces = coarseToFine[coarseI];
316 
317  // Construct single face
319  (
320  IndirectList<face>(patch, fineFaces),
321  patch.points()
322  );
323 
324  if (upp.edgeLoops().size() != 1)
325  {
326  if (fineFaces.size() == 2)
327  {
328  const edge e(fineFaces[0], fineFaces[1]);
329  facePairWeight_[e] = -1.0;
330  }
331  else if (fineFaces.size() == 3)
332  {
333  const edge e(fineFaces[0], fineFaces[1]);
334  const edge e1(fineFaces[0], fineFaces[2]);
335  const edge e2(fineFaces[2], fineFaces[1]);
336  facePairWeight_[e] = -1.0;
337  facePairWeight_[e1] = -1.0;
338  facePairWeight_[e2] = -1.0;
339  }
340 
341  return false;
342  }
343 
344  patchFaces[coarseI] = face
345  (
346  renumber
347  (
348  upp.meshPoints(),
349  upp.edgeLoops()[0]
350  )
351  );
352  }
353 
354  patchLevels_.set
355  (
356  fineLevelIndex,
357  new bPatch
358  (
359  SubList<face>(patchFaces, nCoarseI, 0),
360  patch.points()
361  )
362  );
363 
364  return true;
365 }
366 
367 
369 {
370  label nPairLevels = 0;
371  label nCreatedLevels = 1; // 0 level is the base patch
372  label nCoarseFaces = 0;
373  label nCoarseFacesOld = 0;
374 
375  while (nCreatedLevels < maxLevels_)
376  {
377  const bPatch& patch = patchLevels_[nCreatedLevels - 1];
378  tmp<labelField> finalAgglomPtr(new labelField(patch.size()));
379  bool agglomOK = false;
380 
381  do
382  {
383  label nCoarseFacesPrev = nCoarseFaces;
384 
385  finalAgglomPtr = agglomerateOneLevel
386  (
387  nCoarseFaces,
388  patch
389  );
390 
391  if (nCoarseFaces > 0 && nCoarseFaces != nCoarseFacesPrev)
392  {
393  if
394  (
395  (
396  agglomOK = agglomeratePatch
397  (
398  patch,
399  finalAgglomPtr,
400  nCreatedLevels
401  )
402  )
403  )
404  {
405  restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
406  mapBaseToTopAgglom(nCreatedLevels);
407  setEdgeWeights(nCreatedLevels);
408 
409  if (nPairLevels % mergeLevels_)
410  {
411  combineLevels(nCreatedLevels);
412  }
413  else
414  {
415  nCreatedLevels++;
416  }
417 
418  nPairLevels++;
419  }
420  }
421  else
422  {
423  agglomOK = true;
424  }
425 
426  reduce(nCoarseFaces, sumOp<label>());
427 
428  } while (!agglomOK);
429 
430  nFaces_[nCreatedLevels] = nCoarseFaces;
431 
432  if
433  (
434  !continueAgglomerating(nCoarseFaces)
435  || (nCoarseFacesOld == nCoarseFaces)
436  )
437  {
438  break;
439  }
440 
441  nCoarseFacesOld = nCoarseFaces;
442  }
443 
444  compactLevels(nCreatedLevels);
445 }
446 
447 
448 Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
449 (
450  label& nCoarseFaces,
451  const bPatch& patch
452 )
453 {
454  const label nFineFaces = patch.size();
455 
456  tmp<labelField> tcoarseCellMap(new labelField(nFineFaces, -1));
457  labelField& coarseCellMap = tcoarseCellMap();
458 
459  const labelListList& faceFaces = patch.faceFaces();
460 
461  nCoarseFaces = 0;
462 
463  forAll(faceFaces, facei)
464  {
465  const labelList& fFaces = faceFaces[facei];
466 
467  if (coarseCellMap[facei] < 0)
468  {
469  label matchFaceNo = -1;
470  label matchFaceNeibNo = -1;
471  scalar maxFaceWeight = -GREAT;
472 
473  // Check faces to find ungrouped neighbour with largest face weight
474  forAll(fFaces, i)
475  {
476  label faceNeig = fFaces[i];
477  const edge edgeCommon = edge(facei, faceNeig);
478  if
479  (
480  facePairWeight_[edgeCommon] > maxFaceWeight
481  && coarseCellMap[faceNeig] < 0
482  && facePairWeight_[edgeCommon] != -1.0
483  )
484  {
485  // Match found. Pick up all the necessary data
486  matchFaceNo = facei;
487  matchFaceNeibNo = faceNeig;
488  maxFaceWeight = facePairWeight_[edgeCommon];
489  }
490  }
491 
492  if (matchFaceNo >= 0)
493  {
494  // Make a new group
495  coarseCellMap[matchFaceNo] = nCoarseFaces;
496  coarseCellMap[matchFaceNeibNo] = nCoarseFaces;
497  nCoarseFaces++;
498  }
499  else
500  {
501  // No match. Find the best neighbouring cluster and
502  // put the cell there
503  label clusterMatchFaceNo = -1;
504  scalar clusterMaxFaceCoeff = -GREAT;
505 
506  forAll(fFaces, i)
507  {
508  label faceNeig = fFaces[i];
509  const edge edgeCommon = edge(facei, faceNeig);
510  if
511  (
512  facePairWeight_[edgeCommon] > clusterMaxFaceCoeff
513  && facePairWeight_[edgeCommon] != -1.0
514  && coarseCellMap[faceNeig] > 0
515  )
516  {
517  clusterMatchFaceNo = faceNeig;
518  clusterMaxFaceCoeff = facePairWeight_[edgeCommon];
519  }
520  }
521 
522  if (clusterMatchFaceNo >= 0)
523  {
524  // Add the cell to the best cluster
525  coarseCellMap[facei] = coarseCellMap[clusterMatchFaceNo];
526  }
527  else
528  {
529  // if not create single-cell "clusters" for each
530  coarseCellMap[facei] = nCoarseFaces;
531  nCoarseFaces++;
532  }
533  }
534  }
535  }
536 
537  // Check that all faces are part of clusters,
538 
539  for (label facei=0; facei<nFineFaces; facei++)
540  {
541  if (coarseCellMap[facei] < 0)
542  {
544  (
545  "pairPatchAgglomeration::agglomerateOneLevel "
546  "(label&, const bPatch&) "
547  ) << " face " << facei
548  << " is not part of a cluster"
549  << exit(FatalError);
550  }
551  }
552 
553  return tcoarseCellMap;
554 }
555 
556 
557 void Foam::pairPatchAgglomeration::combineLevels(const label curLevel)
558 {
559  label prevLevel = curLevel - 1;
560 
561  // Set the previous level nCells to the current
562  nFaces_[prevLevel] = nFaces_[curLevel];
563 
564  // Map the restrictAddressing from the coarser level into the previous
565  // finer level
566 
567  const labelList& curResAddr = restrictAddressing_[curLevel];
568  labelList& prevResAddr = restrictAddressing_[prevLevel];
569 
570  forAll(prevResAddr, i)
571  {
572  prevResAddr[i] = curResAddr[prevResAddr[i]];
573  }
574 
575  // Delete the restrictAddressing for the coarser level
576  restrictAddressing_.set(curLevel, NULL);
577 
578  patchLevels_.set(prevLevel, patchLevels_.set(curLevel, NULL));
579 }
580 
581 
582 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
A List with indirect addressing.
Definition: IndirectList.H:102
label maxLevels_
Max number of levels.
const labelListList & faceFaces() const
Return face-face addressing.
EdgeMap< scalar > facePairWeight_
Edge weights.
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:67
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
PtrList< labelField > restrictAddressing_
Cell restriction addressing array.
labelList nFaces_
The number of faces in each level.
Field< label > labelField
Specialisation of Field<T> for label.
Definition: labelField.H:49
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
void resize(const label newSize)
Resize the hash table for efficiency.
Definition: HashTable.C:436
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
label readLabel(Istream &is)
Definition: label.H:64
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const Field< PointType > & localPoints() const
Return pointField of points in patch.
const double e
Elementary charge.
Definition: doubleFloat.H:78
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
void clear()
Clear all entries from table.
Definition: HashTable.C:473
label nFacesInCoarsestLevel_
Number of faces in coarsest level.
label mergeLevels_
Number of levels to merge, 1 = don&#39;t merge, 2 = merge pairs etc.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
label k
Boltzmann constant.
dimensionedScalar cos(const dimensionedScalar &ds)
ListType renumber(const labelUList &oldToNew, const ListType &)
Renumber the values (not the indices) of a list.
Unit conversion functions.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
errorManip< error > abort(error &err)
Definition: errorManip.H:131
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
A list of faces which address into the list of points.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
error FatalError
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
List< label > labelList
A List of labels.
Definition: labelList.H:56
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
A List obtained as a section of another List.
Definition: SubList.H:53
scalar featureAngle_
Feature angle.
const bPatch & patchLevel(const label leveli) const
Return primitivePatch of given level.
void agglomerate()
Agglomerate patch.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
labelList restrictTopBottomAddressing_
Maps from finest to coarsest.
PrimitivePatch< face, List, const pointField > bPatch
PtrList< bPatch > patchLevels_
Hierarchy of patch addressing.
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
A class for managing temporary objects.
Definition: PtrList.H:118
const Field< PointType > & points() const
Return reference to global points.