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