singleCellFvMesh.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-2019 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 "singleCellFvMesh.H"
27 #include "syncTools.H"
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 void Foam::singleCellFvMesh::agglomerateMesh
33 (
34  const fvMesh& mesh,
35  const labelListList& agglom
36 )
37 {
38  // Conversion is a two step process:
39  // - from original (fine) patch faces to agglomerations (aggloms might not
40  // be in correct patch order)
41  // - from agglomerations to coarse patch faces
42 
43  const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
44 
45  // Check agglomeration within patch face range and continuous
46  labelList nAgglom(oldPatches.size(), 0);
47 
48  forAll(oldPatches, patchi)
49  {
50  const polyPatch& pp = oldPatches[patchi];
51  if (pp.size() > 0)
52  {
53  nAgglom[patchi] = max(agglom[patchi])+1;
54 
55  forAll(pp, i)
56  {
57  if (agglom[patchi][i] < 0 || agglom[patchi][i] >= pp.size())
58  {
60  << "agglomeration on patch " << patchi
61  << " is out of range 0.." << pp.size()-1
62  << exit(FatalError);
63  }
64  }
65  }
66  }
67 
68  // Check agglomeration is sync
69  {
70  // Get neighbouring agglomeration
71  labelList nbrAgglom(mesh.nFaces()-mesh.nInternalFaces());
72  forAll(oldPatches, patchi)
73  {
74  const polyPatch& pp = oldPatches[patchi];
75 
76  if (pp.coupled())
77  {
78  label offset = pp.start()-mesh.nInternalFaces();
79  forAll(pp, i)
80  {
81  nbrAgglom[offset+i] = agglom[patchi][i];
82  }
83  }
84  }
85  syncTools::swapBoundaryFaceList(mesh, nbrAgglom);
86 
87 
88  // Get correspondence between this agglomeration and remote one
89  Map<label> localToNbr(nbrAgglom.size()/10);
90 
91  forAll(oldPatches, patchi)
92  {
93  const polyPatch& pp = oldPatches[patchi];
94 
95  if (pp.coupled())
96  {
97  label offset = pp.start()-mesh.nInternalFaces();
98 
99  forAll(pp, i)
100  {
101  label bFacei = offset+i;
102  label myZone = agglom[patchi][i];
103  label nbrZone = nbrAgglom[bFacei];
104 
105  Map<label>::const_iterator iter = localToNbr.find(myZone);
106 
107  if (iter == localToNbr.end())
108  {
109  // First occurrence of this zone. Store correspondence
110  // to remote zone number.
111  localToNbr.insert(myZone, nbrZone);
112  }
113  else
114  {
115  // Check that zone numbers are still the same.
116  if (iter() != nbrZone)
117  {
119  << "agglomeration is not synchronised across"
120  << " coupled patch " << pp.name()
121  << endl
122  << "Local agglomeration " << myZone
123  << ". Remote agglomeration " << nbrZone
124  << exit(FatalError);
125  }
126  }
127  }
128  }
129  }
130  }
131 
132 
133  label coarseI = 0;
134  forAll(nAgglom, patchi)
135  {
136  coarseI += nAgglom[patchi];
137  }
138  // New faces
139  faceList patchFaces(coarseI);
140  // New patch start and size
141  labelList patchStarts(oldPatches.size());
142  labelList patchSizes(oldPatches.size());
143 
144  // From new patch face back to agglomeration
145  patchFaceMap_.setSize(oldPatches.size());
146 
147  // From fine face to coarse face (or -1)
148  reverseFaceMap_.setSize(mesh.nFaces());
149  reverseFaceMap_.labelList::operator=(-1);
150 
151  // Face counter
152  coarseI = 0;
153 
154 
155  forAll(oldPatches, patchi)
156  {
157  patchStarts[patchi] = coarseI;
158 
159  const polyPatch& pp = oldPatches[patchi];
160 
161  if (pp.size() > 0)
162  {
163  patchFaceMap_[patchi].setSize(nAgglom[patchi]);
164 
165  // Patchfaces per agglomeration
166  labelListList agglomToPatch
167  (
168  invertOneToMany(nAgglom[patchi], agglom[patchi])
169  );
170 
171  // From agglomeration to compact patch face
172  labelList agglomToFace(nAgglom[patchi], -1);
173 
174  forAll(pp, i)
175  {
176  label myAgglom = agglom[patchi][i];
177 
178  if (agglomToFace[myAgglom] == -1)
179  {
180  // Agglomeration not yet done. We now have:
181  // - coarseI : current coarse mesh face
182  // - patchStarts[patchi] : coarse mesh patch start
183  // - myAgglom : agglomeration
184  // - agglomToPatch[myAgglom] : fine mesh faces for zone
185  label coarsePatchFacei = coarseI - patchStarts[patchi];
186  patchFaceMap_[patchi][coarsePatchFacei] = myAgglom;
187  agglomToFace[myAgglom] = coarsePatchFacei;
188 
189  const labelList& fineFaces = agglomToPatch[myAgglom];
190 
191  // Create overall map from fine mesh faces to coarseI.
192  forAll(fineFaces, fineI)
193  {
194  reverseFaceMap_[pp.start()+fineFaces[fineI]] = coarseI;
195  }
196 
197  // Construct single face
199  (
200  UIndirectList<face>(pp, fineFaces),
201  pp.points()
202  );
203 
204  if (upp.edgeLoops().size() != 1)
205  {
207  << "agglomeration does not create a"
208  << " single, non-manifold"
209  << " face for agglomeration " << myAgglom
210  << " on patch " << patchi
211  << exit(FatalError);
212  }
213 
214  patchFaces[coarseI++] = face
215  (
216  renumber
217  (
218  upp.meshPoints(),
219  upp.edgeLoops()[0]
220  )
221  );
222  }
223  }
224  }
225 
226  patchSizes[patchi] = coarseI-patchStarts[patchi];
227  }
228 
229  // Pout<< "patchStarts:" << patchStarts << endl;
230  // Pout<< "patchSizes:" << patchSizes << endl;
231 
232  // Compact numbering for points
233  reversePointMap_.setSize(mesh.nPoints());
234  reversePointMap_.labelList::operator=(-1);
235  label newI = 0;
236 
237  forAll(patchFaces, coarseI)
238  {
239  face& f = patchFaces[coarseI];
240 
241  forAll(f, fp)
242  {
243  if (reversePointMap_[f[fp]] == -1)
244  {
245  reversePointMap_[f[fp]] = newI++;
246  }
247 
248  f[fp] = reversePointMap_[f[fp]];
249  }
250  }
251 
252  pointMap_ = invert(newI, reversePointMap_);
253 
254  // Subset used points
255  pointField boundaryPoints(mesh.points(), pointMap_);
256 
257  // Add patches (on still zero sized mesh)
258  List<polyPatch*> newPatches(oldPatches.size());
259  forAll(oldPatches, patchi)
260  {
261  newPatches[patchi] = oldPatches[patchi].clone
262  (
263  boundaryMesh(),
264  patchi,
265  0,
266  0
267  ).ptr();
268  }
269  addFvPatches(newPatches);
270 
271  // Owner, neighbour is trivial
272  labelList owner(patchFaces.size(), 0);
273  labelList neighbour(0);
274 
275 
276  // actually change the mesh
278  (
279  move(boundaryPoints),
280  move(patchFaces),
281  move(owner),
282  move(neighbour),
283  patchSizes,
284  patchStarts,
285  true // syncPar
286  );
287 
288 
289  // Adapt the zones
290  cellZones().clear();
291  cellZones().setSize(mesh.cellZones().size());
292  {
293  forAll(mesh.cellZones(), zoneI)
294  {
295  const cellZone& oldFz = mesh.cellZones()[zoneI];
296 
297  DynamicList<label> newAddressing;
298 
299  // Note: uncomment if you think it makes sense. Note that value
300  // of cell0 is the average.
302  // if (oldFz.localID(0) != -1)
303  //{
304  // newAddressing.append(0);
305  //}
306 
307  cellZones().set
308  (
309  zoneI,
310  oldFz.clone
311  (
312  newAddressing,
313  zoneI,
314  cellZones()
315  )
316  );
317  }
318  }
319 
320  faceZones().clear();
321  faceZones().setSize(mesh.faceZones().size());
322  {
323  forAll(mesh.faceZones(), zoneI)
324  {
325  const faceZone& oldFz = mesh.faceZones()[zoneI];
326 
327  DynamicList<label> newAddressing(oldFz.size());
328  DynamicList<bool> newFlipMap(oldFz.size());
329 
330  forAll(oldFz, i)
331  {
332  label newFacei = reverseFaceMap_[oldFz[i]];
333 
334  if (newFacei != -1)
335  {
336  newAddressing.append(newFacei);
337  newFlipMap.append(oldFz.flipMap()[i]);
338  }
339  }
340 
341  faceZones().set
342  (
343  zoneI,
344  oldFz.clone
345  (
346  newAddressing,
347  newFlipMap,
348  zoneI,
349  faceZones()
350  )
351  );
352  }
353  }
354 
355 
356  pointZones().clear();
357  pointZones().setSize(mesh.pointZones().size());
358  {
359  forAll(mesh.pointZones(), zoneI)
360  {
361  const pointZone& oldFz = mesh.pointZones()[zoneI];
362 
363  DynamicList<label> newAddressing(oldFz.size());
364 
365  forAll(oldFz, i)
366  {
367  label newPointi = reversePointMap_[oldFz[i]];
368  if (newPointi != -1)
369  {
370  newAddressing.append(newPointi);
371  }
372  }
373 
374  pointZones().set
375  (
376  zoneI,
377  oldFz.clone
378  (
379  pointZones(),
380  zoneI,
381  newAddressing
382  )
383  );
384  }
385  }
386 }
387 
388 
389 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
390 
392 (
393  const IOobject& io,
394  const fvMesh& mesh
395 )
396 :
397  fvMesh
398  (
399  io,
400  pointField(), // points
401  faceList(), // faces
402  labelList(), // allOwner
403  labelList(), // allNeighbour
404  false // syncPar
405  ),
406  patchFaceAgglomeration_
407  (
408  IOobject
409  (
410  "patchFaceAgglomeration",
411  io.instance(),
413  *this,
414  io.readOpt(),
415  io.writeOpt()
416  ),
417  0
418  ),
419  patchFaceMap_
420  (
421  IOobject
422  (
423  "patchFaceMap",
424  io.instance(),
426  *this,
427  io.readOpt(),
428  io.writeOpt()
429  ),
430  mesh.boundaryMesh().size()
431  ),
432  reverseFaceMap_
433  (
434  IOobject
435  (
436  "reverseFaceMap",
437  io.instance(),
439  *this,
440  io.readOpt(),
441  io.writeOpt()
442  ),
443  mesh.nFaces()
444  ),
445  pointMap_
446  (
447  IOobject
448  (
449  "pointMap",
450  io.instance(),
452  *this,
453  io.readOpt(),
454  io.writeOpt()
455  ),
456  mesh.nPoints()
457  ),
458  reversePointMap_
459  (
460  IOobject
461  (
462  "reversePointMap",
463  io.instance(),
465  *this,
466  io.readOpt(),
467  io.writeOpt()
468  ),
469  mesh.nPoints()
470  )
471 {
472  const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
473 
474  labelListList agglom(oldPatches.size());
475 
476  forAll(oldPatches, patchi)
477  {
478  agglom[patchi] = identity(oldPatches[patchi].size());
479  }
480 
481  agglomerateMesh(mesh, agglom);
482 }
483 
484 
486 (
487  const IOobject& io,
488  const fvMesh& mesh,
489  const labelListList& patchFaceAgglomeration
490 )
491 :
492  fvMesh
493  (
494  io,
495  pointField(), // points
496  faceList(), // faces
497  labelList(), // allOwner
498  labelList(), // allNeighbour
499  false // syncPar
500  ),
501  patchFaceAgglomeration_
502  (
503  IOobject
504  (
505  "patchFaceAgglomeration",
506  io.instance(),
508  *this,
509  io.readOpt(),
510  io.writeOpt()
511  ),
512  patchFaceAgglomeration
513  ),
514  patchFaceMap_
515  (
516  IOobject
517  (
518  "patchFaceMap",
519  io.instance(),
521  *this,
522  io.readOpt(),
523  io.writeOpt()
524  ),
525  mesh.boundaryMesh().size()
526  ),
527  reverseFaceMap_
528  (
529  IOobject
530  (
531  "reverseFaceMap",
532  io.instance(),
534  *this,
535  io.readOpt(),
536  io.writeOpt()
537  ),
538  mesh.nFaces()
539  ),
540  pointMap_
541  (
542  IOobject
543  (
544  "pointMap",
545  io.instance(),
547  *this,
548  io.readOpt(),
549  io.writeOpt()
550  ),
551  mesh.nPoints()
552  ),
553  reversePointMap_
554  (
555  IOobject
556  (
557  "reversePointMap",
558  io.instance(),
560  *this,
561  io.readOpt(),
562  io.writeOpt()
563  ),
564  mesh.nPoints()
565  )
566 {
567  agglomerateMesh(mesh, patchFaceAgglomeration);
568 }
569 
570 
572 :
573  fvMesh(io),
574  patchFaceAgglomeration_
575  (
576  IOobject
577  (
578  "patchFaceAgglomeration",
579  io.instance(),
581  *this,
582  io.readOpt(),
583  io.writeOpt()
584  )
585  ),
586  patchFaceMap_
587  (
588  IOobject
589  (
590  "patchFaceMap",
591  io.instance(),
593  *this,
594  io.readOpt(),
595  io.writeOpt()
596  )
597  ),
598  reverseFaceMap_
599  (
600  IOobject
601  (
602  "reverseFaceMap",
603  io.instance(),
605  *this,
606  io.readOpt(),
607  io.writeOpt()
608  )
609  ),
610  pointMap_
611  (
612  IOobject
613  (
614  "pointMap",
615  io.instance(),
617  *this,
618  io.readOpt(),
619  io.writeOpt()
620  )
621  ),
622  reversePointMap_
623  (
624  IOobject
625  (
626  "reversePointMap",
627  io.instance(),
629  *this,
630  io.readOpt(),
631  io.writeOpt()
632  )
633  )
634 {}
635 
636 
637 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
638 
639 
640 
641 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
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
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const meshCellZones & cellZones() const
Return cell zones.
Definition: polyMesh.H:482
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
void addFvPatches(const List< polyPatch *> &, const bool validBoundary=true)
Add boundary patches. Constructor helper.
Definition: fvMesh.C:455
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
label nFaces() const
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:67
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:312
List< face > faceList
Definition: faceListFwd.H:43
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
HashTable< label, label, Hash< label > >::const_iterator const_iterator
Definition: Map.H:59
void clear()
Clear the zones.
Definition: MeshZones.C:401
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
const meshPointZones & pointZones() const
Return point zones.
Definition: polyMesh.H:470
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:284
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
List< label > labelList
A List of labels.
Definition: labelList.H:56
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
Foam::polyBoundaryMesh.
void resetPrimitives(pointField &&points, faceList &&faces, labelList &&owner, labelList &&neighbour, const labelList &patchSizes, const labelList &patchStarts, const bool validBoundary=true)
Reset mesh primitive data. Assumes all patch info correct.
Definition: polyMesh.C:672
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:278
fvMesh(const IOobject &io)
Construct from IOobject.
Definition: fvMesh.C:256
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
void setSize(const label)
Reset size of List.
Definition: List.C:281
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
Foam::uindirectPrimitivePatch.
writeOption writeOpt() const
Definition: IOobject.H:363
const fileName & instance() const
Definition: IOobject.H:390
label patchi
const meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:476
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
singleCellFvMesh(const IOobject &io, const fvMesh &)
Construct from fvMesh and no agglomeration.
label nPoints() const
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:430
readOption readOpt() const
Definition: IOobject.H:353
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92